Difficult LAMMPS Problem

Morroa2 · March 31, 2022, 3:14am

I’m having problems with LAMMPS.

The simulation I am working with starts from a restart file. I divide all atoms in the system into a group for liquid atoms, and one for gaseous atoms, based on a Potential Energy criterion.

compute         PotEng all pe/atom

compute         spec all reduce sum c_PotEng
thermo_style    custom step temp pe c_spec
run             0

variable         ifFind atom "c_PotEng < -1.3"
variable        elseFind atom "c_PotEng >= -1.3"

group           liquid13 dynamic all var ifFind every 10
group           gas dynamic all var elseFind every 10

Then I find the center of mass of the liquid.

compute        liquidCenter liquid13 com
variable        liquidCenterX equal c_liquidCenter[1]
variable        liquidCenterY equal c_liquidCenter[2]
variable        liquidCenterZ equal c_liquidCenter[3]

And I apply a force to the liquid and an equal and opposite force to the gas.

variable        liqForceValue equal 0.125
variable        gasForceValue equal "-1.0*v_liqForceValue *(v_liqCount/v_gasCount)"
fix             liqForceadd liquid13 addforce v_liqForceValue 0.0 0.0 every 1000
fix             gasForceAdd gas addforce v_gasForceValue 0.0 0.0 every 1000

Finally, based on the center of mass I calculated and a pre-defined radius, I create a series of “shells” and divide the atoms into groups based on the shells, then record their velocities.

variable        radius  equal 11.23489 #for 0.125 force sphere radius
variable        twiceRad equal "2*v_radius"
variable        threeRad equal "3*v_radius"
variable        fourRad equal "4*v_radius"
variable        fiveRad equal "5*v_radius"
variable        mid1Rad equal "2.01*v_radius"
variable        mid2Rad equal "3.01*v_radius"
variable        mid3Rad equal "4.01*v_radius"

region          deltaTwo sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_twiceRad side out 
region          deltaThree sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_threeRad side out
region          deltaFour sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fourRad side out
region          deltaFive sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fiveRad side out

region          mid1Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid1Rad side in
region          mid2Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid2Rad side in
region          mid3Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid3Rad side in
region          innerSphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_radius

region          intersectSphere1 union 2 deltaTwo innerSphere side out
region          intersectSphere2 union 2 mid1Sphere deltaThree side out
region          intersectSphere3 union 2 mid2Sphere deltaFour side out
region          intersectSphere4 union 2 mid3Sphere deltaFive side out

group           allTwo dynamic all region intersectSphere1 every 10
group           allThree dynamic all region intersectSphere2 every 10
group           allFour dynamic all region intersectSphere3 every 10
group           allFive dynamic all region intersectSphere4 every 10

dump            26 allTwo custom 10000 dumpAllTwo.data id x y z vx vy vz 
dump            326 allThree custom 10000 dumpAllThree.data id x y z vx vy vz
dump            426 allFour custom 10000 dumpAllFour.data id x y z vx vy vz
dump            526 allFive custom 10000 dumpAllFive.data id x y z vx vy vz

I am getting very strange and inconsistent results when trying to graph the shell versus the average x-velocity of atoms in it. Based on the physical conditions of the simulation, I would expect the velocity to grow more negative as the shell increases, but instead it will vary randomly.

Anyone have any ideas how to stop this?

akohlmey · March 31, 2022, 3:41am

First and foremost, you are describing a complex simulation, not necessarily a difficult problem. If you want people to take a closer look, you have to make it simpler and so that you can provide a data file and complete input, so people can look at the whole thing and make experiments/evaluations on their own.

Second, it is not clear from your description why you do not expect large fluctuations at this size scale from instantaneous snapshots.

Morroa2 · March 31, 2022, 5:15pm

#read_restart   liqDrop.01.restart

#units           lj  # lj implies quantities are in LJ internal units
#atom_style      atomic # Defining attributes associated with the atoms
#lattice         sc 0.82 # Defining a simple cubic lattice with atom number density of 0.1
#region          simbox block 0 20 0 20 0 20 # Creating a region called simbox within the limits 0 and 30 in each direction. The format of the limits is as follows: xlo xhi ylo yhi zlo zhi. This is the size of the latticei
#create_box      1 simbox  # Create 1 type of atoms in the simbox
#create_atoms    1 region simbox  # Create atoms of type 1 in the region simbox
#mass            1 1.0     # Defining the mass of atoms of type
#velocity        all create 0.1 87287 # Defining the velocity of the atoms to set the temperature. The final number is a random seed required for the probability distribution while the quantity before that is the temperature

#pair_style      lj/sf 3.0 # Defining the lj potential with a interaction cutoff of 2.5

#pair_coeff      1 1 1.0 1.0 3.0 # Defining the coefficients of the interaction between types of atomsi

#neighbor       0.3 bin # Building neighbour lists using the linked cell method with a skin layer of 0.3
#neigh_modify   every 20 delay 0 check no # Frequency of building neighbour lists

#compute         mxsdcal all msd #Compute msd
#compute        myRDF all rdf 50 # Compute rdf
#fix            4 all ave/time 100 100 100000 c_myRDF file RDFv2.data mode vector # output average rdf at the last timestep (50000) as an average of the last 100 frames in steps of 100
#fix             5 all ave/time 1 1 1 c_msdcal[4] file MSDv2.data mode scalar #4th column of the compute is the total msd


#thermo          1000  # Frequency of outputting thermodynamic data
#thermo_style    custom step temp ke pe etotal vol press pxx pyy pzz  lx ly lz enthalpy # Thermodynamic data to output
timestep        0.0025 # Timestep

compute         PotEng all pe/atom

compute         spec all reduce sum c_PotEng
thermo_style    custom step temp pe c_spec
run             0

variable         ifFind atom "c_PotEng < -1.3"

variable        elseFind atom "c_PotEng >= -1.3"

#group          liquid variable ifFind

group           liquid13 dynamic all var ifFind every 10
group           gas dynamic all var elseFind every 10

variable        cenX equal vcm(liquid13,x)
variable        cenY equal vcm(liquid13,y)
variable        cenZ equal vcm(liquid13,z)


variable        GcenX equal vcm(gas,x)
variable        GcenY equal vcm(gas,y)
variable        GcenZ equal vcm(gas,z)


compute         liqForceComp liquid13 reduce sum fx fy fz
compute         gasForceComp gas reduce sum fx fy fz
compute         allForceComp all reduce sum fx fy fz

set             group all image 0 0 0
compute         liquidCenter liquid13 com

thermo_style    custom step c_liqForceComp[1] c_gasForceComp[1] c_allForceComp[1] c_liquidCenter[1]
run             0

variable        liqForceVar equal c_liqForceComp[1]
variable        gasForceVar equal c_gasForceComp[1]
variable        allForceVar equal c_allForceComp[1]

variable        liqCount equal count(liquid13)
variable        gasCount equal count(gas)

compute         liqTem liquid13 temp
compute         gasTem gas temp

set             group all vx 0 vy 0 vz 0
#set            group all nx 0 ny 0 nz 0


#compute        liquidCenter liquid13 com
variable        liquidCenterX equal c_liquidCenter[1]
variable        liquidCenterY equal c_liquidCenter[2]
variable        liquidCenterZ equal c_liquidCenter[3]
variable        radius  equal 11.23489 #for 0.125 force sphere radius
variable        twiceRad equal "2*v_radius"
variable        mid1Rad equal "2.01*v_radius"
variable        mid2Rad equal "3.01*v_radius"
variable        mid3Rad equal "4.01*v_radius"
variable        threeRad equal "3*v_radius"

variable        fourRad equal "4*v_radius"

variable        fiveRad equal "5*v_radius"

variable        liqForceValue equal 0.125
variable        gasForceValue equal "-1.0*v_liqForceValue *(v_liqCount/v_gasCount)"


fix             liqForceadd liquid13 addforce v_liqForceValue 0.0 0.0 every 1000
fix             gasForceAdd gas addforce v_gasForceValue 0.0 0.0 every 1000

fix             liq liquid13 ave/time 2 5 10 v_cenX v_cenY v_cenZ file liquidVel.data
fix             gas gas ave/time 2 5 10 v_GcenX v_GcenY v_GcenZ file gasVel.data
fix             forceLiq all print 1000 ${liqForceVar} file forceLiq.data
fix             forceGas all print 1000 ${gasForceVar} file forceGas.data
fix             forceAll all print 1000 ${allForceVar} file forceAll.data

region          deltaTwo sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_twiceRad side out #v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_Halfrad
region          deltaThree sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_threeRad side out
#region          deltaThree sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_threeRad side out
region          deltaFour sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fourRad side out
region          deltaFive sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fiveRad side out

region          mid1Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid1Rad side in
region          mid2Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid2Rad side in
region          mid3Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid3Rad side in
region          innerSphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_radius

region          intersectSphere1 union 2 deltaTwo innerSphere side out
region          intersectSphere2 union 2 mid1Sphere deltaThree side out
region          intersectSphere3 union 2 mid2Sphere deltaFour side out
region          intersectSphere4 union 2 mid3Sphere deltaFive side out

group           allTwo dynamic all region intersectSphere1 every 10
group           allThree dynamic all region intersectSphere2 every 10
group           allFour dynamic all region intersectSphere3 every 10
group           allFive dynamic all region intersectSphere4 every 10
#fix             gasMath all print 1000 ${gasCount} file gasMath.data


#compute       PotEng all pe/atom
thermo          100000 # Frequency of outputting thermodynamic data
#thermo_style    custom step c_  #Thermodynamic data to output


variable        atom_pot equal pe #gStore the potential energy of each atom

dump            24 liquid13 custom 10000 dumpLiq.data id x y z
dump            temp liquid13 custom 10000 dumpLiqUnwrapped.data id xu yu zu ix iy iz
dump            25 all custom 10000 dumpAll.data id x y z vx # dump the Coordinates data file in custom format.
dump            26 allTwo custom 10000 dumpAllTwo.data id x y z vx vy vz # dump the Coordinates data file in custom format.
dump            326 allThree custom 10000 dumpAllThree.data id x y z vx vy vz
dump            426 allFour custom 10000 dumpAllFour.data id x y z vx vy vz
dump            526 allFive custom 10000 dumpAllFive.data id x y z vx vy vz

dump            777 all custom 10000 dumpShellCode.data id x y z vx vy vz
#dump           777 all dcd 100000 dumpTest.dcd
#dump_modify    777 unwrap yes
#dump           27 gasTwo custom 10000 dumpGas.data id type x


compute         shell2 allTwo reduce sum vx
compute         shell3 allThree reduce sum vx
compute         shell4 allFour reduce sum vx
compute         shell5 allFive reduce sum vx
#compute         shell2 allTwo reduce mode vx vy vz

# comPos liquid13 ave/time 100 10 1000 c_liquidCenter[1] c_liquidCenter[2] c_liquidCenter[3] file liquidCOM.data

thermo_style     custom step c_liquidCenter[1] c_liquidCenter[2] c_liquidCenter[3]
#fix            cmPrint all print 100000 "$liquidCenterX $liquidCenterY $liquidCenterZ $newX $newY $newZ" file comData.data


dump            28 liquid13 custom 10000 dumpLiqForce.data id fx fy fz
dump            29 gas custom 10000 dumpGasForce.data id fx fy fz
dump            30 all custom 10000 dumpAllForce.data id fx fy fz

restart         100000 liqDrop.*.restart
write_restart   liqDrop.0.restart

#fix             1 all npt temp 0.1 0.1 1 iso 0 0 100
#fix             1 all nve
#fix              1 all nvt temp 0.75 0.75 1 # nvt ensemble. The first number after temp is the initial temperature, the second temperature after temp is the final ensemble temperature and the third number is the temperature damping factor Q.



#change_box       all z final -20 40 y final -20 40 x final -20 40 units lattice #To increase the size of f the boxo
fix              1 all nvt temp 0.75 0.75 1
run              150000

That’s the entire file.

And here is the restart file: https://drive.google.com/file/d/1cGEIlEqVnn0WUnRa614rjG35WMGJnwL_/view?usp=sharing

And to clarify, I am not using instantaneous snapshots, but averages of all velocities from multiple points in time for each shell.

akohlmey · March 31, 2022, 5:32pm

I didn’t ask for the entire file, but for a simplified input deck. Something that demonstrates the issue with the least amount of commands and simulation time. Also, since restart files are in general not portable between versions and architectures, it would need to be an example using a data file anyway.

Morroa2 · March 31, 2022, 8:38pm


timestep        0.0025 # Timestep

compute         PotEng all pe/atom

compute         spec all reduce sum c_PotEng
thermo_style    custom step temp pe c_spec
run             0

variable         ifFind atom "c_PotEng < -1.3"

variable        elseFind atom "c_PotEng >= -1.3"

group           liquid13 dynamic all var ifFind every 10
group           gas dynamic all var elseFind every 10


compute         liquidCenter liquid13 com

thermo_style    custom step c_liqForceComp[1] c_gasForceComp[1] c_allForceComp[1] c_liquidCenter[1]
run             0

variable        liquidCenterX equal c_liquidCenter[1]
variable        liquidCenterY equal c_liquidCenter[2]
variable        liquidCenterZ equal c_liquidCenter[3]

variable        radius  equal 11.23489 #for 0.125 force sphere radius
variable        twiceRad equal "2*v_radius"
variable        mid1Rad equal "2.01*v_radius"
variable        mid2Rad equal "3.01*v_radius"
variable        mid3Rad equal "4.01*v_radius"
variable        threeRad equal "3*v_radius"
variable        fourRad equal "4*v_radius"
variable        fiveRad equal "5*v_radius"

variable        liqForceValue equal 0.125
variable        gasForceValue equal "-1.0*v_liqForceValue *(v_liqCount/v_gasCount)"


fix             liqForceadd liquid13 addforce v_liqForceValue 0.0 0.0 every 1000
fix             gasForceAdd gas addforce v_gasForceValue 0.0 0.0 every 1000

region          deltaTwo sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_twiceRad side out 
region          deltaThree sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_threeRad side out
region          deltaFour sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fourRad side out
region          deltaFive sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_fiveRad side out

region          mid1Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid1Rad side in
region          mid2Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid2Rad side in
region          mid3Sphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_mid3Rad side in
region          innerSphere sphere v_liquidCenterX v_liquidCenterY v_liquidCenterZ v_radius

region          intersectSphere1 union 2 deltaTwo innerSphere side out
region          intersectSphere2 union 2 mid1Sphere deltaThree side out
region          intersectSphere3 union 2 mid2Sphere deltaFour side out
region          intersectSphere4 union 2 mid3Sphere deltaFive side out

group           allTwo dynamic all region intersectSphere1 every 10
group           allThree dynamic all region intersectSphere2 every 10
group           allFour dynamic all region intersectSphere3 every 10
group           allFive dynamic all region intersectSphere4 every 10

variable        atom_pot equal pe #Store the potential energy of each atom

dump            24 liquid13 custom 10000 dumpLiq.data id x y z
dump            25 all custom 10000 dumpAll.data id x y z vx 
dump            26 allTwo custom 10000 dumpAllTwo.data id x y z vx vy vz 
dump            326 allThree custom 10000 dumpAllThree.data id x y z vx vy vz
dump            426 allFour custom 10000 dumpAllFour.data id x y z vx vy vz
dump            526 allFive custom 10000 dumpAllFive.data id x y z vx vy vz



thermo          100000 # Frequency of outputting thermodynamic data
thermo_style     custom step c_liquidCenter[1] c_liquidCenter[2] c_liquidCenter[3]

fix              1 all nvt temp 0.75 0.75 1
run              150000

Here’s a more slimmed down version of the code. I will add in the data shortly.

akohlmey · March 31, 2022, 8:43pm

This doesn’t look much reduced. there are still many regions and lots of commented out lines that are distracting.

Morroa2 · March 31, 2022, 8:47pm

The issue is only really visible over several regions.

akohlmey · March 31, 2022, 8:57pm

Why? Why not just for one of those shells? and how should I get that information?
I only see four dump commands and dynamic groups, but no compute reduce/region which I would have expected as a means to get the average velocity of atoms in a region/shell.