NPT simulation with RATTLE holonomic constraints

Dear all,

I am trying to use lammps (nov 3, 2022 version) to do NPT simulations with holonomic constraints. As mentioned in the documentation, when velocity Verlet is used as the EoM integrator, RATTLE algorithm has to be prefered over SHAKE.

I started my investigations using the rhodo bench. The input in.rhodo file uses the two following fixes, in this order:

fix             1 all shake 0.0001 5 0 m 1.0 a 232
fix             2 all npt temp 300.0 300.0 100.0 &
                z 0.0 0.0 1000.0 mtk no pchain 0 tchain 1

Everything goes fine in this case.
Then I tried to use the RATTLE fix instead of SHAKE:

fix             2 all npt temp 300.0 300.0 100.0 &
                z 0.0 0.0 1000.0 mtk no pchain 0 tchain 1
fix             1 all rattle 0.0001 5 0 m 1.0 a 232

Note that, according to the documentation on fix rattle, it should be called after all other integration fixes. When I run this script, I obtain the following error:

ERROR: Fix rattle must come before any box changing fix (src/RIGID/fix_shake.cpp:365)

Indeed, NPT algorithm changes box dimensions…

Then, I tried to apply the rattle fix before the npt fix (analogously with the shake case):

fix             1 all rattle 0.0001 5 0 m 1.0 a 232
fix             2 all npt temp 300.0 300.0 100.0 &
               z 0.0 0.0 1000.0 mtk no pchain 0 tchain 1

I then get a WARNING and an ERROR:

rattle-npt.log:WARNING: Fix rattle should come after all other integration fixes  (src/RIGID/fix_rattle.cpp:138)
...
rattle-npt.log:ERROR on proc 6: Out of range atoms - cannot compute PPPM (src/KSPACE/pppm.cpp:1887)

The WARNING is totally in agreement with the documentation. The ERROR must be a consequence of the wrong order of the 2 fixes.

Finally, I have 2 problems:

  • if the fix rattle is placed after the fix npt, I get an error “Fix rattle must come before any box changing fix…” else if the fix rattle is placed before the fix npt, I obtain a warning “Fix rattle should come after all other integration fixes…”. Thus, the rattle fix can’t be used before nor after the npt fix.

  • more importantly, is it possible to use rattle with npt?

I (tried to) upload a pack.tgz archive with 4 different in.rhodo files and log files corresponding to most of the cases mentioned here, plus an output of my lmp -h command line (but new users can not upload).

Thank you.
Bernard Rousseau

You are between a rock and a hard place.
Either you must not use a barostat at all (i.e. fix nvt instead of fix npt) or use fall back to using fix shake (it is not that bad in most applications).

Edit: you could also consider using fix press/berendsen, but then there is the question of whether what you gain with using fix rattle over fix shake you might be losing with that choice over a nose-hoover barostat.

Dear Axel,
Thank you for your response and the time spent on it.

Either you must not use a barostat at all (i.e. fix nvt instead of fix npt) or use fall back to using fix shake (it is not that bad in most applications).

Well, I need to know density at given P and T, so… I can get it without constraints and a small timestep, but for comparison purposes I wanted to do it with rattle and npt fixes. I will give a try to the shake fix, then.
One more point here. Martyna et al. Molec. Phys 1996 proposed an iterative scheme to implement RATTLE/ROLL with explicit reversible NPT integrator. Is there any technical reason why this in not implemented in lammps (appart the obvious human and time resources reason)?

Edit: you could also consider using fix press/berendsen, but then there is the question of whether what you gain with using fix rattle over fix shake you might be losing with that choice over a nose-hoover barostat…

I tried that but also get into many troubles. The following fixes in this order

fix             1 all nvt temp 300.0 300.0 100.0
fix             2 all press/berendsen iso 0.0 0.0 1000.0
fix             3 all rattle 0.0001 5 0 m 1.0 a 232

give the error message (as expected):

ERROR: Fix rattle must come before any box changing fix (src/RIGID/fix_shake.cpp:365)

while with

fix             1 all nvt temp 300.0 300.0 100.0
fix             2 all rattle 0.0001 5 5 m 1.0 a 232
fix             3 all press/berendsen iso 0.0 0.0 1000.0

the code starts but never ends (only step 0 is done). With shake, the total duration is roughly 15 seconds while here, I’ve been about 30 mn without any change in the log (reproduced below).

LAMMPS (3 Nov 2022)
  using 2 OpenMP thread(s) per MPI task
using multi-threaded neighbor list subroutines
Reading data file ...
  orthogonal box = (-27.5 -38.5 -36.3646) to (27.5 38.5 36.3615)
  2 by 5 by 2 MPI processor grid
  reading atoms ...
  32000 atoms
  reading velocities ...
  32000 velocities
  scanning bonds ...
  4 = max bonds/atom
  scanning angles ...
  8 = max angles/atom
  scanning dihedrals ...
  18 = max dihedrals/atom
  scanning impropers ...
  2 = max impropers/atom
  reading bonds ...
  27723 bonds
  reading angles ...
  40467 angles
  reading dihedrals ...
  56829 dihedrals
  reading impropers ...
  1034 impropers
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
    12 = max # of 1-3 neighbors
    24 = max # of 1-4 neighbors
    26 = max # of special neighbors
  special bonds CPU = 0.005 seconds
  read_data CPU = 0.539 seconds
Finding RATTLE clusters ...
    1617 = # of size 2 clusters
    3633 = # of size 3 clusters
     747 = # of size 4 clusters
    4233 = # of frozen angles
  find clusters CPU = 0.003 seconds
PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
  G vector (1/distance) = 0.24883488
  grid = 25 32 32
  stencil order = 5
  estimated absolute RMS force accuracy = 0.035547797
  estimated relative force accuracy = 0.00010705113
  using double precision FFTW3
  3d grid and FFT values/proc = 4914 1600
Generated 2278 of 2278 mixed pair_coeff terms from arithmetic mixing rule
Last active /omp style is kspace_style pppm/omp
Neighbor list info ...
  update: every = 1 steps, delay = 5 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 10 13 13
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/charmm/coul/long/omp, perpetual
      attributes: half, newton on, omp
      pair build: half/bin/newton/omp
      stencil: half/bin/3d
      bin: standard
Setting up Verlet run ...
  Unit style    : real
  Current step  : 0
  Time step     : 2
RATTLE stats (type/ave/delta/count) on step 0
Bond:    12   1.07999   0.000176421      252
Bond:    22   1.1       0.000185329      396
Bond:    28   1.08      0.000136975       20
Bond:    32   1.111     0.000174163       80
Bond:    33   1.11099   0.000117405       40
Bond:    37   1.08299   0.000101864        6
Bond:    39   1.06998   0.000121966        6
Bond:    46   1.08099   0                  1
Bond:    48   1.08101   0                  1
Bond:    50   1.08101   0                  1
Bond:    54   1.08098   0                  1
Bond:    55   1.081     0                  1
Bond:    63   1.081     0                  1
Bond:    65   1.081     0                  1
Bond:    70   1.11099   0.000155271      108
Bond:    71   1.07999   0.000169132      305
Bond:    78   1.111     0.000179398      626
Bond:    79   1.07999   0.000149891       46
Bond:    88   1.11099   0.000196405      717
Bond:    91   1.111     0.00016334        99
Bond:    95   1.111     0.000245892     6194
Bond:    96   1.08      0.000162783      198
Bond:    99   1.111     0.000204859      600
Bond:   100   1.08      0.000225048      891
Bond:   102   0.996992  0.000218916      327
Bond:   103   0.999995  0.000173504       54
Bond:   104   1         0.00016159        12
Bond:   105   0.97599   3.19739e-05        5
Bond:   106   0.960002  0.000106691       60
Bond:   107   0.999991  0.000156041       35
Bond:   108   1.04001   0.000129251       33
Bond:   109   1.01301   0                  1
Bond:   110   1.32501   0.000100997        6
Bond:   112   0.957197  0.000272687     8466
Angle:  232   104.52    0.0237494       4233
Per MPI rank memory allocation (min/avg/max) = 21.88 | 22.15 | 22.51 Mbytes
------------ Step              0 ----- CPU =            0 (sec) -------------
TotEng   =    -25356.5617 KinEng   =     21444.4759 Temp     =       299.0347 
PotEng   =    -46801.0377 E_bond   =      2537.9940 E_angle  =     10921.3742 
E_dihed  =      5211.7865 E_impro  =       213.5116 E_vdwl   =     -2307.8634 
E_coul   =    207025.8927 E_long   =   -270403.7333 Press    =      -142.1259 
Volume   =    307995.0335

Any idea?
Thank you again for the time spent,
Bernard

Try a larger relaxation time constant for fix press/berendsen.

Most features in LAMMPS are implemented and contributed by users that need and care about them. The support for RATTLE is a good example. It was done by a LAMMPS user that was computing properties where the (subtle) difference between RATTLE and SHAKE mattered.
Before that, nobody needed it enough to spend the time to implement it. Similarly, the support for being able to use fix shake/rattle during minimization was added because it got old to explain to beginners how you can just replace the constraints with stiff force constants during minimization.

I tried 10 times more for the relaxation constant

fix             3 all press/berendsen iso 0.0 0.0 10000.0

and get the same issues than with rattle/npt:

-rw-r--r-- 1 rousseau users 6289909 12 déc.  20:06 data.rhodo
Bond:   105   0.97599   3.19739e-05        5
Bond:   106   0.960002  0.000106691       60
Bond:   107   0.999991  0.000156041       35
Bond:   108   1.04001   0.000129251       33
Bond:   109   1.01301   0                  1
Bond:   110   1.32501   0.000100997        6
Bond:   112   0.957197  0.000272687     8466
Angle:  232   104.52    0.0237494       4233
Per MPI rank memory allocation (min/avg/max) = 21.88 | 22.15 | 22.51 Mbytes
------------ Step              0 ----- CPU =            0 (sec) -------------
TotEng   =    -25356.5617 KinEng   =     21444.4759 Temp     =       299.0347
PotEng   =    -46801.0377 E_bond   =      2537.9940 E_angle  =     10921.3742
E_dihed  =      5211.7865 E_impro  =       213.5116 E_vdwl   =     -2307.8634
E_coul   =    207025.8927 E_long   =   -270403.7333 Press    =      -142.1259
Volume   =    307995.0335
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1798)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1798)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1798)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1798)

Which means that fix press/berendsen is changing the box too much with these settings, so try an even larger time constant.

Well, I may end up in an NVT simulation :slight_smile:

OK. I will try.

Please see the following note in the fix press/berendsen docs:

The relaxation time is actually also a function of the bulk modulus of the system (inverse of isothermal compressibility). The bulk modulus has units of pressure and is the amount of pressure that would need to be applied (isotropically) to reduce the volume of the system by a factor of 2 (assuming the bulk modulus was a constant, independent of density, which it’s not). The bulk modulus can be set via the keyword modulus. The Pdamp parameter is effectively multiplied by the bulk modulus, so if the pressure is relaxing faster than expected or desired, increasing the bulk modulus has the same effect as increasing Pdamp. The converse is also true. LAMMPS does not attempt to guess a correct value of the bulk modulus; it just uses 10.0 as a default value which gives reasonable relaxation for a Lennard-Jones liquid, but will be way off for other materials and way too small for solids. Thus you should experiment to find appropriate values of Pdamp and/or the modulus when using this fix.

Yes, I can remember this. As I was working on a liquid and that many codes are based on water isothermal compressibility, I thought the value for rhodo was fine.

A large value of

fix             3 all press/berendsen iso 0.0 0.0 10000000.0

did work for rhodo :slight_smile:

To sum up: I can either try shake/npt or nvt+rattle+berendsen barostat. And if needed, I can try to implement Martyna stuff…

Thank you for your help.
Bernard