Issues with increasing timestep

LAMMPS LAMMPS General Discussion
1

Hello all users,

I have been simulating 1000 molecules water system. I’m using TIP3P water model and using lj/cut/coul/long potential to run NPT and then NVT simulation (the units are set to real & I’m setting periodic boundary conditions). After this I save the system to a restart file and later run it with metadynamics. For running metadynamics I’ve been using PLUMED (file attached along with the input script below). The simulation runs perfectly with 1fs timestep but when I try to run the same metadynamics simulation with 2fs timestep, I’m observing large fluctuations in pressure and temperature and can’t obtain gaussian distributed total energy as I do for 1fs.

I’ve tried changing Tdamp for the thermostat but I wasn’t able to get any good results. Would accept any suggestions/feedback on how to make the simulation work with the new timestep.

Thanks,
Nisarg

Input script

read_restart nvt.restart.5000000
kspace_style pppm 1.0e-5

#variables
variable sysvol equal vol
variable sysmass equal mass(all)/6.0221367e+23
variable sysdensity equal v_sysmass/v_sysvol/1.0e-24
variable etotal equal etotal
variable pe equal pe
variable ke equal ke
variable p equal press
variable time equal step*dt+0.000001

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

reset_timestep 0
thermo_style custom step v_time press vol v_sysdensity temp etotal pe ke
thermo 10
thermo_modify norm no flush yes

fix 3 all nve
fix 4 all temp/csvr 300.0 300.0 100.0 54324
fix 5 all ave/time 100 1 100 v_time c_thermo_temp v_p v_sysdensity v_etotal v_pe v_ke file wte_out.txt
fix 6 all plumed plumedfile plumed.dat outfile p_metad.log
timestep 2

#forces
compute 1 all property/atom fx fy fz
dump 1 all custom 100 force_wte.dump id type fx fy fz
dump_modify 1 sort id

dump 2 all custom 100 water_wte.dump id type x y z
dump_modify 2 sort id
restart 9000000 nvt.restart

run 2500000

Plumed script

UNITS ENERGY=kcal/mol

energy: ENERGY

METAD …

LABEL=metad

ARG=energy

HEIGHT=0.6

BIASFACTOR=32

TEMP=300.0

SIGMA=35

FILE=HILLS

PACE=125

… METAD

PRINT ARG=energy STRIDE=250 FILE=COLVAR

2

I don’t see you using fix shake or fix rattle anywhere.
TIP3P is a rigid water model and requires one of the two anyway and a 2fs timestep is not stable without constraining bonds involving hydrogen atoms.

Without constraints you have a different water model and even 1fs is too large a timestep (with typical values for bond and angle force constants you can use 0.25fs to 0.5fs with reasonable energy conservation).

In general, you should be running a test of the system without metadynamics and without a thermostat (after reaching equilibrium) and monitor energy conservation. The simulation may run without crashing at larger than reasonable timesteps, but your energy conservation will be broken and total energy will quickly go through the roof. A thermostat may mask that, but that doesn’t make the simulation any more correct.

3

Thanks for your help! I was under the impression that fix shake information was saved in the restart file but after reading the documentation I realized my mistake. I am able to get the results after your suggestions.

1 Like
4

It is always helpful to review the documentation :wink:

1 Like