VasiliBaranov / packing-generation

Hard-sphere packing generation in C++ with the Lubachevsky–Stillinger, Jodrey–Tory, and force-biased algorithms and packing post-processing.
MIT License
107 stars 44 forks source link

loose FCC as initial configuration #14

Open hfrusinque opened 5 years ago

hfrusinque commented 5 years ago

Hi Vasili,

I want to approximate the rearrangement after relaxation of a loose FCC configuration, from a colloidal crystal whose strong electrostatic forces are "suddenly shut down". I would like to use your code (namely the LS algorithm). Should I provide the initial positions in "packing_init.xyzd" in Little-Endian format or can the code "deal" with ASCII as input?

PS: I normally change the format of the final positions via command line with: hexdump -v -e '4 "%f "' -e '"\n"' < packing.xyzd > centers_in_ascii.dat

Best regards, Hector Rusinque PhD Candidate TU Clausthal

hfrusinque commented 5 years ago

after reading your documentation in this Matlab file and building on this Matlab function, I can generate a file "packing.xyzd" with the initial positions of an FCC crystal in Little-Endian format and it works at least for the ideal FCC configuration. Now, the problem is when I want to get a loose FCC configuration with a packing density lower than the FCC one (i.e. lower than ~0.74). The code stops always after reading the packing with "'packing.xyzd' was opened". I suppose one of the conditions for a valid initial configuration is that the spheres are in contact, maybe even that there is some degree of overlapping. I will try by adding some noise to the positions of the spheres to get some degree of contact/overlapping.

Here is the function to get the initial configuration so far:

function FCCLattice = FCCMake(a,V, finalDensity)
    % "a" is the Lattice Parameter:
    a = a; # (e.g. a=1.0)
    %Set the Volume of the cube in variable V:
    V = V; % (#e.g. V=3.0) This means your cube will be of dimensions VxVxV.
    %%%----Code-----%%%
    %Create an array that stores all your coordinates:
    %First you will need to establish the number of atoms
    %that the array will have, we will do this by establishing
    %the known number of atoms, in this case, FCC 4xVxVxV.
    N = V*V*V*4;
    %This array will give a matrix of Nx3, that will store data.
    FCCLattice = zeros(N,3);
    %Create the vectors for the atoms position in the FCC Lattice:
    FCCatoms = [0 0 0;(a/2) (a/2) 0;(a/2) 0 (a/2);0 (a/2) (a/2)];
    %Set a variable to change rows in the array storing the coordinates:
     n = 0;
     k=0;
     FCC_Density = pi/(3*sqrt(2));
     R_FCC = sqrt(2)*a/4.0 #Radius for FCC
     finalDensity = 0.73;
     R = R_FCC*(finalDensity / FCC_Density)^(1/3) # final Radius
    %Create 3 For loops to mix the positions between xyz
    for x = 0:V-1
        for y = 0:V-1
            for z = 0:V-1
                for i = 1:4
                    %Create a vector that translate your location in the
                    %coordinate system into the position in the crystal:
                    coordinatestranslation = a*[x y z];
                    %Add 1 to move one row down in the array:
                    n = n+1;
                    FCCLattice(n,:) = coordinatestranslation + FCCatoms(i,:);             
                end
                 k = k + 12 + 4;
                 M_xyzd(k-15)  = FCCLattice(n-3,1);
                 M_xyzd(k-14)  = FCCLattice(n-3,2);
                 M_xyzd(k-13)  = FCCLattice(n-3,3);
                 M_xyzd(k-12)  = 2*R;   

                 M_xyzd(k-11)  = FCCLattice(n-2,1);
                 M_xyzd(k-10)  = FCCLattice(n-2,2);
                 M_xyzd(k-9)  = FCCLattice(n-2,3);
                 M_xyzd(k-8)  = 2*R;    

                 M_xyzd(k-7)  = FCCLattice(n-1,1);
                 M_xyzd(k-6)  = FCCLattice(n-1,2);
                 M_xyzd(k-5)  = FCCLattice(n-1,3);
                 M_xyzd(k-4)  = 2*R;

                 M_xyzd(k-3)  = FCCLattice(n,1);
                 M_xyzd(k-2)  = FCCLattice(n,2);
                 M_xyzd(k-1)  = FCCLattice(n,3);
                 M_xyzd(k)    = 2*R;
            end
        end
    end

    floatDataType = 'float64';
    LittleEndian = "ieee-le";

    save packing_ascii.txt M_xyzd;

    FCC_binary = fopen('packing.xyzd','wb'); # FCC_binary: File ID

    fwrite(FCC_binary, M_xyzd, floatDataType, LittleEndian);

    fclose(FCC_binary);
    end
VasiliBaranov commented 5 years ago

Hi Hector,

you are right, one has to fill in the packing.xyzd file. packing_init.xyzd is just a backup.

That's strange that the program stops so abruptly. Could you please attach a complete log? Please make sure that there is no file packing.nfo in this folder, it is a marker of a finished program (if i recall correctly).

Also (if i recall correctly), the LS algorithm implementation rescales the particles at the beginning of execution so that the closest pair of particles touches. In case of a FCC packing (even if it was rescaled before running the program), it leads to a completely jammed configuration.

To avoid rescaling, you can try to run the program with the option -md, see https://github.com/VasiliBaranov/packing-generation#3-post-processing , option 6. I think it does not rescale the particles (i can check if you need), but it computes some other things that you might not need, so it might be a little slower.

Best Regards, Vasili

hfrusinque commented 5 years ago

Thanks for your reply, Vasili. It was a rookie mistake, I forgot to delete the file "packing.nfo" before running the program. It works perfectly. For perfectly monodisperse particles nothing interesting happens, LS simply rescales the radii until reaching the FCC packing (the completely jammed configuration you mentioned with \phi~ 0.74). It becomes interesting though when a certain degree of polydispersity is added.

I tried the -md option with a small packing (108 spheres) I produced after running the -lsgd algorithm but it stopped with the following warning message:

decreasing the diameters will lead to unexpected consequences (e.g., incorrect density).

The -md log file is attached. log_md.txt

The theoretical and calculated densities were very similar, so the rescaling step must not have had a strong effect. It would be nice to use the -md option as you suggested to suppress the rescaling of the radii.

Calc. porosity is 0.349714 Theoretical porosity is 0.344944

Liebe Grüße, Hector

VasiliBaranov commented 5 years ago

Dear Hector,

it seems (and i forgot the details already) that to run with the -md option you need to rescale particle radii to the calculated density, as explained in https://github.com/VasiliBaranov/packing-generation#21-note-on-final-diameters . The difference in densities that you see can be quite significant for some applications. Here is a line in the script that does it: https://github.com/VasiliBaranov/packing-generation/blob/master/Docs/MATLAB%20scripts%20for%20interpreting%20results/ReadPackingScript.m#L44 (and please see the explanation above).

After rescaling, you can try to run with the -md option again (the packing.nfo file shall be present this time, but you have it after the -lsgd run).

Hope this helps.

Beste Grüße :-) Vasili