lammps fix rigid/nvt cannot generate correct temperature

LAMMPS LAMMPS Mailing List Mirror
1

Dear All,

I have a system that only has water molecules. I modeled the water with TIP3P and TIP4P model. I am running NVT simulations using the command ‘fix 1 all nvt temp 10 10 0.01’ for tip3p and command ‘fix 1 all rigid/nvt molecule temp 10 10 0.01’ for tip4p. And I use ‘velocity all create 10 123456789 units box dist gaussian’ to generate the initial velocity after the fix command is issued. However, I noticed that the tip4p system is not generating the correct temperature at step zero while the tip3p system looks fine (see both outputs below). Additionally, in the tip4p case, there is a warning message “Cannot count rigid body degrees-of-freedom before bodies are fully initialized (…/fix_rigid.cpp:1144)” in the output file. Presumably, LAMMPS is not able to correctly account for the number of degrees of freedom in the system. This error in temperature evaluation is present not only at time zero but at all times and therefore impacts the dynamics through the thermostat. I would appreciate any advice anyone has on how to “fully initialize the bodies” so that LAMMPS properly accounts for the degrees of freedom.

Best,
Qianping

TIP4P******

LAMMPS (17 Jun 2013)
clear
atom_style full
units metal
dimension 3
boundary p p p
read_data tip4p.dat
3 = max bonds/atom
3 = max angles/atom
orthogonal box = (-1.1558 -1.1558 -1.1558) to (119.855 119.855 119.855)
1 by 1 by 1 MPI processor grid
184000 atoms
138000 bonds
138000 angles
3 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
2 = max # of 1-4 neighbors
3 = max # of special neighbors

Define interaction parameters

pair_style lj/cut/coul/long 15.0
kspace_style pppm 1.0e-6

pair_coeff * * 0 0 #all other pairs
pair_coeff 1 1 0.006722 3.154 #OW-OW
pair_coeff 1 2 0.000000 0.000 #OW-HW
pair_coeff 2 2 0.000000 0.000 #HW-HW

bond_style harmonic
bond_coeff 1 26.016 0.9572 #water oh bond(oplsrigid-updated)
bond_coeff 2 26.016 0.15 #water om bond(oplsrigid-updated)

angle_style harmonic
angle_coeff 1 3.252 104.52 #OPLS
angle_coeff 2 3.252 52.26 #OPLS

group water type 1 2 3
184000 atoms in group water

run_style verlet

fix 1 all rigid/nvt molecule temp 10 10 0.01
1 rigid bodies with 184000 atoms
neigh_modify delay 0 every 1 check yes page 200000 one 20000
velocity all create 10 123456789 units box dist gaussian
WARNING: Cannot count rigid body degrees-of-freedom before bodies are fully initialized (…/fix_rigid.cpp:1144)
compute h2o_inter water group/group water kspace yes
compute wintrape water pe/atom bond angle
compute h2o_intra water reduce sum c_wintrape
compute totpe all pe
compute h2o_ke water ke

thermo_style custom step temp c_h2o_ke c_h2o_inter c_h2o_intra c_totpe etotal press vol
thermo_modify lost warn flush yes line one
thermo 1000
dump mydump001 all xyz 1000 dump001.xyz
dump_modify mydump001 element O H Hm
run 100000
PPPM initialization …
G vector (1/distance) = 0.227608
grid = 120 120 120
stencil order = 5
estimated absolute RMS force accuracy = 1.50194e-05
estimated relative force accuracy = 1.04304e-06
using double precision FFTs
3d grid and FFT values/proc = 2048383 1728000
Memory usage per processor = 1359.37 Mbytes
Step Temp h2o_ke h2o_inte h2o_intr totpe TotEng Press Volume
0 19.193294 0.002480928 -34192.522 0.000911214 -34192.521 -34192.518 -182939.49 1772035.4

2

Dear All,

  I have a system that only has water molecules. I modeled the water
with TIP3P and TIP4P model. I am running NVT simulations using the command 'fix
1 all nvt temp 10 10 0.01' for tip3p and command 'fix 1 all rigid/nvt
molecule temp 10 10 0.01' for tip4p. And I use 'velocity all create 10
123456789 units box dist gaussian' to generate the initial velocity after
the fix command is issued. However, I noticed that the tip4p system is not
generating the correct temperature at step zero while the tip3p system
looks fine (see both outputs below). Additionally, in the tip4p case,
there is a warning message "Cannot count rigid body degrees-of-freedom
before bodies are fully initialized (../fix_rigid.cpp:1144)" in the
output file. Presumably, LAMMPS is not able to correctly account for the
number of degrees of freedom in the system. This error in temperature
evaluation is present not only at time zero but at all times and therefore
impacts the dynamics through the thermostat. I would appreciate any advice
anyone has on how to "fully initialize the bodies" so that LAMMPS properly
accounts for the degrees of freedom.

​the first thing you have to do when encountering a problem with LAMMPS
using an old(er) version is to try with the current version. LAMMPS is a
large package with many features and thus it cannot be bug free, but if
bugs are discovered and verified, they usually get fixed quickly and
patches are released immediately. so to make sure you are not encountering
a problem that has already been corrected, please check against the current
version.

in order to trigger a complete initialization of the rigid bodies, i would
suggest to assign the initial velocity distribution only after you did a
"run 0"

axel.

3

Dear Axel,

Thank you very much for your reply. I took your advise and run the simulation again with the current version (1Feb2014) and I also did a ‘run 0’ before I assign the initial velocity distribution. The good news is that the warning message is gone. However, the tip4p system is still not generating the assigned initial temperature (10 K). Based on the above facts, it seems that there are other things messed up during the temperature calculation process. We’ve been looking into this problem for weeks, and we desperately need help on how to solve this issue. Please kindly help us with this. Thank you very much as always.

Best,
Qianping

4

Dear Axel,

Thank you very much for your reply. I took your advise and run the
simulation again with the current version (1Feb2014)

​the *real* current version ​is 13 May 2014. 1Feb2014 has quite a few known
bugs that have since been fixed.

and I also did a 'run 0' before I assign the initial velocity
distribution. The good news is that the warning message is gone. However,
the tip4p system is still not generating the assigned initial temperature
(10 K). Based on the above facts, it seems that there are other things
messed up during the temperature calculation process. We've been looking
into this problem for weeks, and we desperately need help on how to solve
this issue. Please kindly help us with this. Thank you very much as
always.

Best,
Qianping

# initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

​*here* is your problem. you are simulating your entire box of water as a
single rigid body. this *cannot* work and LAMMPS behaves as to be expected.
and if i remember correctly, this was pointed out to you before in a
previous post. bad input => bad simulation.

axel.

​[...]​

5

Dear Axel,

Thank you very much for your reply. I took your advise and run the simulation again with the current version (1Feb2014)

​the real current version ​is 13 May 2014. 1Feb2014 has quite a few known bugs that have since been fixed.

and I also did a ‘run 0’ before I assign the initial velocity distribution. The good news is that the warning message is gone. However, the tip4p system is still not generating the assigned initial temperature (10 K). Based on the above facts, it seems that there are other things messed up during the temperature calculation process. We’ve been looking into this problem for weeks, and we desperately need help on how to solve this issue. Please kindly help us with this. Thank you very much as always.

Best,
Qianping

initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

​here is your problem. you are simulating your entire box of water as a single rigid body. this cannot work and LAMMPS behaves as to be expected. and if i remember correctly, this was pointed out to you before in a previous post. bad input => bad simulation.
thanks for the prompt reply. I did correct this problem and I got exactly the same result. And even if I simulate with a single rigid body, it should at least gives me the assigned temperature at step 0 whatsoever. I am truly confused by how lammps works. Also, the last email you wrote to me is on May.8th, I am sorry that I didn’t notice there’s a new version coming out. I will test the code again and see if that fixes my problem.

LAMMPS (1 Feb 2014)
clear
atom_style full
units metal
dimension 3
boundary p p p
read_data tip4p.dat
orthogonal box = (-1.1558 -1.1558 -1.1558) to (119.855 119.855 119.855)
1 by 1 by 1 MPI processor grid
reading atoms …
184000 atoms
scanning bonds …
3 = max bonds/atom
scanning angles …
3 = max angles/atom
reading bonds …
138000 bonds
reading angles …
138000 angles
3 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
2 = max # of 1-4 neighbors
3 = max # of special neighbors

Define interaction parameters

pair_style lj/cut/coul/long 15.0
kspace_style pppm 1.0e-6
pair_coeff * * 0 0 #all other pairs
pair_coeff 1 1 0.006722 3.154 #OW-OW
pair_coeff 1 2 0.000000 0.000 #OW-HW
pair_coeff 2 2 0.000000 0.000 #HW-HW
bond_style harmonic
#bond_style none
bond_coeff 1 26.016 0.9572 #water oh bond(oplsrigid-updated)
bond_coeff 2 26.016 0.15 #water om bond(oplsrigid-updated)
angle_style harmonic
#angle_style none
angle_coeff 1 3.252 104.52 #OPLS
angle_coeff 2 3.252 52.26 #OPLS
group water type 1 2 3
184000 atoms in group water
run_style verlet

initialize

fix 1 water rigid/nvt molecule temp 10 10 0.1
46000 rigid bodies with 184000 atoms
run 0
PPPM initialization …
G vector (1/distance) = 0.227608
grid = 120 120 120
stencil order = 5
estimated absolute RMS force accuracy = 1.50194e-05
estimated relative force accuracy = 1.04304e-06
using double precision FFTs
3d grid and FFT values/proc = 2048383 1728000
Memory usage per processor = 1171.01 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -34192.522 0.0009112142 -34192.521 -19027.655
Loop time of 3.09944e-06 on 1 procs for 0 steps with 184000 atoms

Pair time () = 0 (0) Bond time () = 0 (0)
Kspce time () = 0 (0) Neigh time () = 0 (0)
Comm time () = 0 (0) Outpt time () = 0 (0)
Other time (%) = 3.09944e-06 (100)

FFT time (% of Kspce) = 0 (0)
FFT Gflps 3d (1d only) = 0 0

Nlocal: 184000 ave 184000 max 184000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 202398 ave 202398 max 202398 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1.95448e+08 ave 1.95448e+08 max 1.95448e+08 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 195448084
Ave neighs/atom = 1062.22
Ave special neighs/atom = 3
Neighbor list builds = 0
Dangerous builds = 0
neigh_modify delay 0 every 1 check yes page 200000 one 20000
velocity all create 10 123456789 units box dist gaussian
compute h2o_inter water group/group water kspace yes
compute wintrape water pe/atom bond angle
compute h2o_intra water reduce sum c_wintrape
compute totpe all pe
compute h2o_ke water ke

------------- Equilibration and thermalisation ----------------

thermo_style custom step temp c_h2o_ke c_h2o_inter c_h2o_intra c_totpe etotal press vol
thermo_modify lost warn flush yes line one
thermo 1
dump mydump001 all xyz 1000 dump001.xyz
dump_modify mydump001 element O H Hm
run 100000
PPPM initialization …
G vector (1/distance) = 0.227608
grid = 120 120 120
stencil order = 5
estimated absolute RMS force accuracy = 1.50194e-05
estimated relative force accuracy = 1.04304e-06
using double precision FFTs
3d grid and FFT values/proc = 2048383 1728000
Memory usage per processor = 1367.57 Mbytes
Step Temp h2o_ke h2o_inte h2o_intr totpe TotEng Press Volume
0 0 0 -34192.522 0.0009112142 -34192.521 -34192.521 -30617.941 1772035.4
1 0.41393981 4.9224911 -34197.473 0.0009112142 -34197.472 -34192.55 -16829.017 1772035.4
2 1.6173436 19.233133 -34211.871 0.0009112142 -34211.87 -34192.637 -16848.219 1

axel.

​[…]​

6

And also, I would like to point out that my tip3p system which I follow exactly the same rule to form the input file were able to generate the correct temperature. And my colleges are having the same issue with temperature when they model the water as tip4p. I sincerely hope you could have a careful look at this issue.

Dear Axel,

Thank you very much for your reply. I took your advise and run the simulation again with the current version (1Feb2014)

​the real current version ​is 13 May 2014. 1Feb2014 has quite a few known bugs that have since been fixed.

and I also did a ‘run 0’ before I assign the initial velocity distribution. The good news is that the warning message is gone. However, the tip4p system is still not generating the assigned initial temperature (10 K). Based on the above facts, it seems that there are other things messed up during the temperature calculation process. We’ve been looking into this problem for weeks, and we desperately need help on how to solve this issue. Please kindly help us with this. Thank you very much as always.

Best,
Qianping

initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

​here is your problem. you are simulating your entire box of water as a single rigid body. this cannot work and LAMMPS behaves as to be expected. and if i remember correctly, this was pointed out to you before in a previous post. bad input => bad simulation.
thanks for the prompt reply. I did correct this problem and I got exactly the same result. And even if I simulate with a single rigid body, it should at least gives me the assigned temperature at step 0 whatsoever. I am truly confused by how lammps works. Also, the last email you wrote to me is on May.8th, I am sorry that I didn’t notice there’s a new version coming out. I will test the code again and see if that fixes my problem.

LAMMPS (1 Feb 2014)
clear
atom_style full
units metal
dimension 3
boundary p p p
read_data tip4p.dat
orthogonal box = (-1.1558 -1.1558 -1.1558) to (119.855 119.855 119.855)
1 by 1 by 1 MPI processor grid
reading atoms …
184000 atoms
scanning bonds …
3 = max bonds/atom
scanning angles …
3 = max angles/atom
reading bonds …
138000 bonds
reading angles …
138000 angles
3 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
2 = max # of 1-4 neighbors
3 = max # of special neighbors

Define interaction parameters

pair_style lj/cut/coul/long 15.0
kspace_style pppm 1.0e-6
pair_coeff * * 0 0 #all other pairs
pair_coeff 1 1 0.006722 3.154 #OW-OW
pair_coeff 1 2 0.000000 0.000 #OW-HW
pair_coeff 2 2 0.000000 0.000 #HW-HW
bond_style harmonic
#bond_style none
bond_coeff 1 26.016 0.9572 #water oh bond(oplsrigid-updated)
bond_coeff 2 26.016 0.15 #water om bond(oplsrigid-updated)
angle_style harmonic
#angle_style none
angle_coeff 1 3.252 104.52 #OPLS
angle_coeff 2 3.252 52.26 #OPLS
group water type 1 2 3
184000 atoms in group water
run_style verlet

initialize

fix 1 water rigid/nvt molecule temp 10 10 0.1
46000 rigid bodies with 184000 atoms
run 0
PPPM initialization …
G vector (1/distance) = 0.227608
grid = 120 120 120
stencil order = 5
estimated absolute RMS force accuracy = 1.50194e-05
estimated relative force accuracy = 1.04304e-06
using double precision FFTs
3d grid and FFT values/proc = 2048383 1728000
Memory usage per processor = 1171.01 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -34192.522 0.0009112142 -34192.521 -19027.655
Loop time of 3.09944e-06 on 1 procs for 0 steps with 184000 atoms

Pair time () = 0 (0) Bond time () = 0 (0)
Kspce time () = 0 (0) Neigh time () = 0 (0)
Comm time () = 0 (0) Outpt time () = 0 (0)
Other time (%) = 3.09944e-06 (100)

FFT time (% of Kspce) = 0 (0)
FFT Gflps 3d (1d only) = 0 0

Nlocal: 184000 ave 184000 max 184000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 202398 ave 202398 max 202398 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1.95448e+08 ave 1.95448e+08 max 1.95448e+08 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 195448084
Ave neighs/atom = 1062.22
Ave special neighs/atom = 3
Neighbor list builds = 0
Dangerous builds = 0
neigh_modify delay 0 every 1 check yes page 200000 one 20000
velocity all create 10 123456789 units box dist gaussian
compute h2o_inter water group/group water kspace yes
compute wintrape water pe/atom bond angle
compute h2o_intra water reduce sum c_wintrape
compute totpe all pe
compute h2o_ke water ke

------------- Equilibration and thermalisation ----------------

thermo_style custom step temp c_h2o_ke c_h2o_inter c_h2o_intra c_totpe etotal press vol
thermo_modify lost warn flush yes line one
thermo 1
dump mydump001 all xyz 1000 dump001.xyz
dump_modify mydump001 element O H Hm
run 100000
PPPM initialization …
G vector (1/distance) = 0.227608
grid = 120 120 120
stencil order = 5
estimated absolute RMS force accuracy = 1.50194e-05
estimated relative force accuracy = 1.04304e-06
using double precision FFTs
3d grid and FFT values/proc = 2048383 1728000
Memory usage per processor = 1367.57 Mbytes
Step Temp h2o_ke h2o_inte h2o_intr totpe TotEng Press Volume
0 0 0 -34192.522 0.0009112142 -34192.521 -34192.521 -30617.941 1772035.4
1 0.41393981 4.9224911 -34197.473 0.0009112142 -34197.472 -34192.55 -16829.017 1772035.4
2 1.6173436 19.233133 -34211.871 0.0009112142 -34211.87 -34192.637 -16848.219 1

axel.

​[…]​

7

​[...]​

# initialize

fix 1 all rigid/nvt molecule temp 10 10 0.1

1 rigid bodies with 184000 atoms

​*here* is your problem. you are simulating your entire box of water as a
single rigid body. this *cannot* work and LAMMPS behaves as to be expected.
and if i remember correctly, this was pointed out to you before in a
previous post. bad input => bad simulation.

thanks for the prompt reply. I did correct this problem and I got exactly
the same result. And even if I simulate with a single rigid body, it
should at least gives me the assigned temperature at step 0 whatsoever. I
am truly confused by how lammps

​you have a combination of problems. one is now fixed, the other is a bit
trickier and from running my own tests, i see now that there are still some
open issues with how LAMMPS handles counting the degrees of freedom for
rigid bodies and updates the internal rigid body data. there have been
significant improvements to the rigid body​ and related code and for some
situations, there seem to be some oversights.

works. Also, the last email you wrote to me is on May.8th, I am sorry
that I didn't notice there's a new version coming out. I will test the
code again and see if that fixes my problem.

LAMMPS follows the "frequent updates" philosophy, and thus ​there are new
versions coming out frequently; as soon as bugs are fixed or features
improved, a new version/patch is made available.​ the changes and releases
are documented here: http://lammps.sandia.gov/bug.html

8

And also, I would like to point out that my tip3p system which I follow
exactly the same rule to form the input file were able to generate the
correct temperature. And my colleges are having the same issue with
temperature when they model the water

​i think this is a misconception. for tip3p you are likely using fix shake
(i would) instead of fix rigid, and that is a cat of a different color. the
kind of problem that i am talking about is specific ​to various fix rigid
versions.

as tip4p. I sincerely hope you could have a careful look at this issue.

​please let me point out that there are several simple workarounds that
would work well​, if you would not be so fixated on one little detail and a
bit more creative in how you set up your simulations. anyway, it seems that
this is a case, where a systematic look at multiple initialization and
equilibration scenarios will be necessary and that requires a discussion
the issue with the main man.

i'd also like to remind you that with the very latest version of LAMMPS,
there is now a fix rigid/nvt/small which should result in a much better
parallel performance for your kind of system than fix rigid/nvt.

for now, i suggest you experiment with using the run 0 hack or skipping it.
or using a short run where you use plain fix rigid in combination with fix
temp/csvr and a rather short time constant to quickly crank up the
temperature to where you want (fix temp/csvr is the same as mixing in a
fraction of a velocity command repeatedly during the run and thus should
result in a cleaner initialization than using the velocity command just
once), then unfix both and switch to fix rigid/nvt or rigid/nvt/small.

axel.

9

Thank you very much. I will go ahead and test those out.

10

Thank you very much. I will go ahead and test those out.

​in addition to those suggestions, i also have been wondering if there is a
special reason why you are using fix rigid and not one of the tip4p pair
styles and fix shake. that would allow you a much larger time step and none
of the issues that you are experiencing (which doesn't mean that they don't
need to be corrected, only that you currently seem to make your life much
harder than it need to be).

axel.​