How to create randomly sized spheres within lammps by command lines or reading a data file

[MevrouwA]

Hi, Currently I am working on creating a 3D system of X number of spheres with random radius between 2 different values. I have consulted everything about this, but I can’t solve this myself and hope someone can help me with this. I will briefly explain what I have done so far. I’ve defined an input script with all sections: I use lennard jones potential and atom_style sphere, because I’m interested in the positions and diameter.

Then I defined properties as follows:
variable rmin equal 1.0 # minimum value for the radius
variable rmax equal 3.0 # maximum value for the radius
variable rdifference equal {rmax}-{rmin}
variable radius equal {{rmin}+{rdifference}*random()}

create_box 1 simulation_box
create_atoms 1 random ${N} 12345 NULL overlap 1.0 maxtry 1

group all type 1
set group all diameter ${radius}

When running the simulation, the random generator function “random()” does not work and I get “ERROR: Variable radius: Invalid syntax in variable formula”. I don’t understand why? So then I tried to solve this with the following line: variable radius equal ceil(random(0,3,1)), where 0 is minimum value, 3 is maximum value, and 1 represents seed (which I don’t understand by the way?).
This works, but all the spheres keep the same constant radius. In addition, the value remains the same even though I run it several times.

Then I tried a loop:

Loop over the spheres and create them

Label loop_spheres
variable i loop {N} # number of iterations variable radius equal ceil(random(1,3,1)) variable x equal ceil(random(1,20,1)) variable y equal cell(random(1,20,1)) variable z equal cell(random(1,20,1)) create_atoms 1 sphere {x} {y} {z} ${radius}
next i
jump SELF loop_spheres
The “create_atoms …” line gives the following error: Unknown create_atoms command option sphere. If I delete this line nothing happens, so the radius and positions do not change!

Finally, I created a data_file.txt myself with python script and used the following format for atom_style sphere:

atom-ID atom-type diameter density x y z
1 1 5.48 1 1 9 1
2 1 4.40 1 7 7 0
3 1 2.10 1 11 19 1
4 1 4.03 1 17 4 3
5 1 2.74 1 5 2 9
6 1 5.32 1 6 17 13
7 1 2.49 1 16 12 13
8 1 2.21 1 19 10 10
9 1 3.54 1 2 5 11
10 1 2.06 1 4 13 9

After the initialization, defining of sphere properties, system definition and simulation settings, I added the following line command to read the file: read_data data_read_in.txt add append
with data_read_in.txt name of file and ‘add append’ to indicate that simulation box is defined.

Error that now follows:
ERROR: Unexpected end of data file (src/read_data.cpp:1968)
Last command: read_data data_read_in.txt add append

I don’t know what else to do and hope someone can help me with this? I would also like to know when the file is read. How do I then link those values to the variable radius and x,y,z positions to create randomly sized spheres?

[MevrouwA]

label loop_spheres
variable i loop {N}
radius equal ceil(random(1,3,1))
variable x equal ceil(random(1,20,1))
variable y equal cell(random(1,20,1))
variable z equal cell(random(1,20,1))
create_atoms 1 sphere {x} {y} {z} {radius}
next i
jump SELF loop_spheres

[akohlmey]

@MevrouwA

I use lennard jones potential and atom_style sphere, because I’m interested in the positions and diameter.

That does not make sense. The Lennard-Jones potential pays no attention to the per-atom diameter setting. The effective diameter of a Lennard-Jones particles (\sigma) is set via the pair_coeff parameter and is a per-atom-type property. If you want a Lennard-Jones potential that uses the atom radius as \sigma, you or one of the LAMMPS developers would have to write a new pair style.

@MevrouwA

variable radius equal {{rmin}+rmin+{rmin}+{rdifference}*random()}

create_box 1 simulation_box
create_atoms 1 random ${N} 12345 NULL overlap 1.0 maxtry 1

group all type 1
set group all diameter ${radius}

This also does not make sense. If you want to have the different atoms to have different diameters and you want to take it from a randomly assigned value, you need to define an atom style variable, not equal style. The documentation of the “set” command explicitly says so.

@MevrouwA

group all type 1

The group “all” is predefined and cannot be changed.

@MevrouwA

When running the simulation, the random generator function “random()” does not work and I get “ERROR: Variable radius: Invalid syntax in variable formula”. I don’t understand why?

that is because your variable expression is in violation of the syntax as it is described in the documentation for the variable command. Simple as that. You cannot just make up a syntax you like. You have to exactly follow what is documented. If you stray from that, you get an error.

Here is a simple input example to create 250 spheres with a radius ranging randomly from 1.0\sigma to 3.0\sigma:

units lj
atom_style sphere
region box block 0 25 0 25 0 25
create_box 1 box

create_atoms 1 random 250 5256326 box overlap 3.0
variable diam atom random(2.0,6.0,3425634)
set group all diameter v_diam

write_data random_size_spheres.data nocoeff

[akohlmey]

If you want to use different LJ sigma values for the different atoms without writing a custom pair style, you have to assign a new atom type to each atom and set a different value of \sigma.
This can be done with something like the following:

variable num index 250
atom_style atomic
units lj
region box block 0 25 0 25 0 25
create_box ${num} box

create_atoms 1 random ${num} 5256326 box overlap 3.0
variable newtype atom id
set group all type v_newtype

pair_style lj/cut 5.0
variable sigma equal random(1.0,3.0,3425634)
variable rmass equal 4.0/3.0*PI*v_sigma*v_sigma*v_sigma

variable r loop ${num}
label radius
    pair_coeff ${r} ${r} 1.0 ${sigma}
    mass ${r} ${rmass}
next r
jump SELF radius

write_data random_size_lj.data

[MevrouwA]

Thank you very much for your response and all the explanations!

I have never worked with lammps before, so therefore it is sometimes a bit hard to understand everything well. But now I’ve managed to make a 3D system with randomly sized spheres. However, despite adding the ‘Overlap’ keyword, the spheres still overlap. Now I’ve seen that with delete_command you can remove overlapping spheres. But I prefer that if I make 250 spheres, the number will not change either. Is there any tips on how I can fix this?

Unfortunately, I don’t understand the last part about the different LJ sigma. After running the simulation, the radius and number of spheres remain constant. I understand this is because atom style atomic is applied and therefore you can’t set the diameter attribute. But I don’t understand what you exactly want to reach with that part of code, could you maybe explain it a bit more please?

Beforehand I said I’m interested in radius and positions. By using a force field as lj potential and applying a certain velocity to the spheres, the spheres will move and interact with each other and get into a position that will match a determined global porosity. When they are in equilibrium, I want to use the x, y and z coordinates and the associated radii from the data file to draw the spheres in another program (Fluent, spaceclaim). I hope it is more clear now.

However with atomic style, I do not get any information about the diameters of the spheres in the dump file.

[akohlmey]

@MevrouwA

But now I’ve managed to make a 3D system with randomly sized spheres. However, despite adding the ‘Overlap’ keyword, the spheres still overlap. Now I’ve seen that with delete_command you can remove overlapping spheres. But I prefer that if I make 250 spheres, the number will not change either. Is there any tips on how I can fix this?

You are asking the wrong question here. Either LAMMPS is not doing what the documentation says it will do, or you are not doing in your input what is necessary to have the desired conditions applied. To be able to tell which is the case, you need to provide a (minimal, like what I posted) example that reproduces the issue and explain how exactly you have determined that what you are looking at is the way you think. It is impossible to give specific advice otherwise.

@MevrouwA

By using a force field as lj potential and applying a certain velocity to the spheres

You have to understand (and looks to me that this is what you are still not fully on board with) that when using the lj/cut pair style, the atom radius is completely ignored and irrelevant. Indeed Lennard-Jones particles are rather soft and “squishy” and thus have at best an “effective radius”. The size of that is based on the value of \sigma and that value is not a per-atom, but a per-type property. If you look at the second example, you see that there are 250 atom types (one for each atom) and then there are potential parameters created for each type, so in the dump file, you don’t have a way to output the sigma parameter, since it can only output per-atom data (and it doesn’t know that there are just as many types as atoms).

It is actually possible to store the sigma value as a custom property with the atoms by using “fix property/atom”:

fix Radius all property/atom d_radius

variable r loop ${num}
label radius
    pair_coeff ${r} ${r} 1.0 ${sigma}
    mass ${r} ${rmass}
    set atom ${r} d_radius ${sigma}
next r
jump SELF radius

The d_radius per-atom property can then be output with a dump file and used for visualization purposes.