MD simulation with TIP4P water model

Dear LAMMPS users,

Hello.

I have tried to run the simulation with TIP4P2005 water model.

After several trial and error, I can run the simulation successfully without error or warning.

The shake works well, which make the bond length and bond angle of TIP4p water molecules constant.

PPPM/tip4p also seems to work fine and temperature is well maintained.

But the problem is that the water molecules don’t move.

When I look into the log file and thermo data, those have same values for each step.

Even though temperature of system is increased to 350K, the result is same.

I can’t find where the problem exist.

Is there any fault in my input data or input script?

(please find the attached files).

Thank you very much.

Best regards,
Sehun Joo

====================================================================

Define variables

====================================================================

input data file

variable lmp_data string data.water

temperature and pressure

variable T equal 350

variable p equal 3/0.101325

interactions / rcut

variable rcut equal 12

number of steps / timestep

variable nstep equal “10000000”

variable dt equal 2

output

variable freq_thermo equal “100” # thermo output every

variable freq_frame equal “10000” # frame output every

variable freq_snap equal “v_nstep/10” # xyz output every for safety (rerun)

====================================================================

Initialization

====================================================================

units real

boundary p p p

dimension 3

atom_style full

====================================================================

Force Field

====================================================================

read_data ${lmp_data}

#replicate 5 5 5

1 WH

2 WO

mass 1 1.008

mass 2 15.9994

pair_style lj/cut/tip4p/long 2 1 1 1 0.1546 ${rcut}

pair_modify mix arithmetic

pair_coeff 1 1 0 0

pair_coeff 1 2 0 0

pair_coeff 2 2 0.1852 3.1589

1 WH-WO

bond_style harmonic

bond_coeff 1 900 0.9572

5 WH-WO-WH (fix)

angle_style harmonic

angle_coeff 1 124 104.52

Coulomb interaction

kspace_style pppm/tip4p 1.0e-4

#special_bonds lj/coul 0.0 0.0 0.0

group o type 2

group water type 1 2

====================================================================

Constraints

====================================================================

fix fixtip4p all shake 0.0001 20 10 b 1 a 1

====================================================================

1. RELAX / NVT

====================================================================

velocity all create $T 87654 dist gaussian mom yes rot yes

fix thermostat all temp/berendsen {T} {T} 100.0

#fix barostat all press/berendsen aniso {p} {p} 1000.0

thermo ${freq_thermo}

thermo_style custom step cpuremain temp press etotal pe ke epair ecoul evdwl ebond eangle edihed

thermo_modify flush yes

dump dtrj all custom 100 ${lmp_data}.lammpstrj id type x y z mol

dump dsnap all custom {freq_frame} {lmp_data}_*.lammpstrj id type x y z mol

dump dxyz all xyz 100 ${lmp_data}_*.xyz

dump_modify dxyz element H O

reset_timestep 0

timestep ${dt}

run 50000

undump dtrj

undump dsnap

write_restart restart.RELAX01

print “01.RELAX done”

data.water (30.7 KB)

in.water (2.8 KB)

Dear LAMMPS users,

Hello.

I have tried to run the simulation with TIP4P2005 water model.

After several trial and error, I can run the simulation successfully without
error or warning.

The shake works well, which make the bond length and bond angle of TIP4p
water molecules constant.

PPPM/tip4p also seems to work fine and temperature is well maintained.

But the problem is that the water molecules don’t move.

When I look into the log file and thermo data, those have same values for
each step.

Even though temperature of system is increased to 350K, the result is same.

I can’t find where the problem exist.

Is there any fault in my input data or input script?

i don't see any fix that does time integration. only fix shake and the
berendsen thermostat. without time integration atoms don't move.

please see the first warning on the fix temp/berendsen documentation:

http://lammps.sandia.gov/doc/fix_temp_berendsen.html

axel.

An additional warning:

The ensemble of conformations sampled by the Berendsen thermostat
algorithm is not the canonical ensemble. (The conformations have the
correct average total energy, but if I recall correctly, the
distribution of (total) energies is typically too narrow.) This will
give you incorrect results if you use statistical reweighting
techniques later (such as importance sampling, WHAM, etc...) In most
cases, it also prevents you from using replica exchange method (REMD).
This caused a lot of headaches for us a few years ago.

I would avoid the Berendsen algorithm.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2699284/

Andrew