CHARMM TIP3P

When I calculate the temperature of water in a region, the output temperature was alway lower than the expected one.

The input file was like this:

units real
atom_style full
boundary p p p

read_data 00.data.lammps

group water type 1 2

bond_style harmonic
bond_coeff 1 450 0.9572
angle_style harmonic
angle_coeff 1 55 104.52
pair_style lj/cut/coul/long 12.0
pair_coeff 1 1 0.1521 3.15061 12.000
pair_coeff 1 2 0.0836 1.77530 12.000
pair_coeff 2 2 0.0460 0.40000 12.000

kspace_style pppm 0.0001
special_bonds lj/coul 0.0 0.0 0
neighbor 2.0 bin
neigh_modify delay 0 one 2000 check yes

timestep 1.0

region temp1 block INF INF INF INF -15 15 units box
compute Tregion1 water temp/region temp1

velocity water create 300.0

min_style cg
min_modify dmax 1.0
minimize 0.01 0.001 1000 1000000
thermo_style custom step temp etotal press c_Tregion1

thermo_modify flush yes line one

thermo 100

fix regidwater water shake 0.0001 10 10 b 1 a 1

fix 2 water npt temp 300.0 300.0 100.0 z 1.00 1.00 1000.0
run 500000
unfix 2

I found it is because the presence of ‘fix regidwater water shake 0.0001 10 10 b 1 a 1’. If this command was deleted, temperature of region ‘temp1’ was 300K.

Why? Are there something wrong with shake?

Thanks,

Hangyan

2011-09-24

When I calculate the temperature of water in a region, the output
temperature was alway lower than the expected one.

[...]

I found it is because the presence of
'fix regidwater water shake 0.0001 10 10 b 1 a 1'. If this
command was deleted, temperature of region 'temp1' was 300K.

Why? Are there something wrong with shake?

no. not at all. PEBCAC!
remember that fix shake takes away degrees of freedom.
a water molecule has 3N = 9 DOF, however you constrain
two bonds and an angle and thus have only 6 DOF.
i bet that makes exactly the difference. have you looked at
the internal compute that npt fix creates?

axel.