I am simulating 1000 molecules of FLEXIBLE water with Dreiding force field (without hydrogen bonding) in the most recent version of lammps. My initial box has been prepared at 1 g/cc. First, minimization works fine (potential energy is lowered to -3000/4000 units)
After that i tried NPT (300 K, 1 atm) to anneal the system. Even though i could maintain temperature close to 300 and pressure close to 1 (barring high initial fluctuations), my volume kept increasing to ridiculously high number, which would make water density much lower than 1 even after 250000 timesteps (0.5 fmsec each).
I tried various combination of parameters, but it is the same. Then i tried simulating it in NVT, noting that initial volume is the ideal one. While the volume and temperature are okay, pressure rises to 9000 atm in 5000/10000 timesteps.
I am not able to understand, why i am not getting fluctuations in volume (for NPT) and in pressure (for NVT) around desired means such as 1 g/cc volume (for NPT) and 1 atm pressure for NVT?
It has been really surprising to me that i am having trouble with such a simple system. Does anybody has any leads?
THanks in advance
Following is the code i am using :
water with Dreiding hydrogen bonding potential
Define basics of the system
neighbor 2.0 bin
neigh_modify every 1 delay 1
bond_coeff 1 350.000000 0.980000
bond_coeff 2 350.000000 0.980000
angle_coeff 1 50.000000 104.510000
pair_style lj/cut/coul/long 12.0
pair_coeff 1 1 0.0001 2.846
pair_coeff 2 2 0.0957 3.033
pair_coeff 1 2 0.0031 2.938
kspace_style ewald 1.0e-4
pair_modify shift yes
group water molecule <> 1 1000
thermo_style custom step cpu press temp pe ke etotal enthalpy evdwl ecoul epair ebond eangle edihed emol elong vol
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all npt temp 300.0 300.0 50.0 iso 0.0 0.0. 100
restart 50000 restartf2.*
run 250000 upto
sorry for the another email again, but by mistake the pressure entered in my earlier email code is 0.0 but i simulated for 1 atm as well.
I am simulating 1000 molecules of FLEXIBLE water with Dreiding force field
(without hydrogen bonding) in the most recent version of lammps. My initial
I am not able to understand, why i am not getting fluctuations in volume
(for NPT) and in pressure (for NVT) around desired means such as 1 g/cc
volume (for NPT) and 1 atm pressure for NVT?
... and i am surprised that you expect a simple generic classical
force field to be that accurate. classical forcefields are always
compromises and compared to most other molecules, water
has a high dipole moment. the symptoms you describe correspond
to a system that is interacting too weakly. as an experiment, you
could try to raise the charges on oxygen and hydrogen until
you get the desired interactions strength. it will probably be
a bad model to simulate water, but that is how some water
potentials are derived (empirically that is).
It has been really surprising to me that i am having trouble with such a
simple system. Does anybody has any leads?
yes. have a look at the literature! there are a wealth of papers on
potentials for water interactions. people have been trying to come
up with good interactions potentials for water for over 40 years now,
and there are still new models being proposed.
even using quantum chemistry, it is not easy to get the properties
of water right (water simulated using the BLYP or PBE functionals
freezes at room temperature, if you can afford to run long enough)
bottom line: don't expect miracles from something that tries to be
I don't know why it's happening, but the 2 simulations are consistent.
If the NVT is equilibrating to a high pressure, then when you run
NPT at a target pressure of 0, that means the box will tend to
expand. I'm guessing there is something invalid about your
model that is causing the pressure to not be what you expect. Maybe
Maybe it is a stupid suggestion, but could it be an issue related to
units? If you are working in real units, then water density should be
on the order of 0.03 A^-3.
Have you first tried with a more common/simple model of water (SPC/E, TIP3P) ?
Hi Axel, Steve and Laurent
Thanks for the reply.
@Axel: Yes, you are right about charges. I was earlier using Gasteiger charges (Q-H=0.21). Now i am using HF charges (Q-H=0.41) and now the density is above 0.96 g/cc at 1 atm. It is a good point to be noted.
The reason i am using Dreiding is that i wanted to test hydrogen bonding potential of DREIDING. It is claimed in that paper that gasteiger and experimental charges work, but in my simulations with DREIDING, only HF charges (the highest possible) work.
@ Steve: My NPT simulations are at 1 atm and not zero. Anyway, I tried changing tail corrections and so on, but the results are same as above.
@Laurent: As i specifically wanted to test dreding hydrogen bonding potential, i am doing that.