[lammps-users] thermostat/barostat problems with water

LAMMPS LAMMPS Mailing List Mirror
1

LAMMPS folks,
I am trying to judge the limits of the range of applicability of water potentials found in the literature. However, I am having trouble using certain combinations of thermostats and barostats cited in these papers. I am using water that was first equilibrated by the NPT thermostat. Could someone please tell me why they are failing?

The error is reported in the thermo output:

Step Temp peperh2o Press density
190000 249.77017 -12.826013 130.91496 0.89062524
190001 298 -12.310568 1100.8691 0.89062524
190002 298 nan nan 0.89062524
ERROR: Computed temperature for fix temp/rescale cannot be 0.0

BEGIN INPUT FILE:
units real
dimension 3
boundary p p p
atom_style full
read_restart 298K.1Bar.restart.190000
#read_data Ice40.dat

## alternate way to set charges - beside manually in file ###
group ox type 2
group hy type 1
set group ox charge -0.8476
set group hy charge 0.4238

### Flexible SPC/E Potential Parameters ###
### Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216 ###
pair_style lj/cut/coul/cut 15.5
pair_modify tail yes #long range correction to lj
pair_coeff 2 2 0.1502629 3.1169
pair_coeff 1 2 0.0341116368 2.04845
pair_coeff 1 1 0.00774378 0.98
bond_style harmonic
bond_coeff 1 176.864 0.9611
angle_style harmonic
angle_coeff 1 42.1845 109.4712

#the order of the ensemble and thermostat/barostat commands might be important!
fix ensemble all nve
fix ptherm all press/berendsen iso 1.0 1.0 1000.0
fix ttherm all temp/rescale 1 298.0 298.0 0.01 1.0

#Compute
compute totalpe all pe

#Variables
variable peperh2o equal c_totalpe/count(ox) #total energy per H2O molecule
variable tempv equal vol
variable v equal ${tempv}
variable density equal count(ox)/v_v*29.91526 #density of water
#Computes
v_keperO file ave.spc.dat
compute myRDF all rdf 500 1 1 1 2 2 2
fix 3 all ave/time 2 100 5000 c_myRDF file zhang.rdf mode vector

neighbor 2.0 bin
neigh_modify delay 0 every 10 check yes
thermo_modify norm no flush yes press ptherm_press temp ttherm_temp #order of press & temp call is very important
thermo_style custom step temp v_peperh2o press v_density
thermo 1

#dumps
dump waterdump all atom 5000 298K.1Bar.dump
restart 10000 298K.1Bar.restart

#run variables
timestep 0.75
run 200000

2

christopher,

LAMMPS folks,
I am trying to judge the limits of the range of applicability of water potentials found in the literature. However, I am having trouble using certain combinations of thermostats and barostats cited in these papers. I am using water that was first equilibrated by the NPT thermostat. Could someone please tell me why they are failing?

to call temp/rescale a "thermostat" is a bit of a stretch.
in any case, have you tried running with a shorter time step?
0.75 fs seems pretty large for a flexible water potential.

axel.

3

If I add/delete these lines to examples/peptide/in.peptide

#fix 1 all nvt temp 275.0 275.0 100.0 tchain 1
#fix 2 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
fix ttherm all temp/rescale 1 298.0 298.0 0.01 1.0
fix nve all nve

it runs fine, with a 2.0 fmsec timestep. That system is mostly
water (small peptide). So there must be something perverse
about your system. From your output it looks like you
ran the first 190K steps with something else, then it croaked
within 2 timesteps of starting to use fix temp/rescale every
step. To rescale the T every step is very dramatic. Try
doing it less often and more gradually.

Steve