# Energy drift in SPC/E

Dear all

I model water using SPC/E with periodic BC in all directions to obtain RDFs. Using NVT, I obtained equilibrium at 300K (2ns) and then measure RDFs by NVE ensemble for 7ns. However, water temperature was not constant and increased each time step. Being specific, in 7ns temperature increased 300K to 316K. My input and output are simply like following. My theory is there is an energy drift due to the shake algorithm. I believe that during the process of keeping the constrains constant, there is a force exerted by shake algorithm which appears as a heat dissipation. Am I doing something wrong or shake does that?

You may also find my water video on youtube with following link.

Murat Barisik

########################INPUT

units metal
dimension 3
atom_style full

boundary p p p
neighbor 2.0 bin
bond_style harmonic
angle_style harmonic

pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 1 1 0.00674 3.16556
pair_coeff 1*2 2 0.0 0.0

kspace_style pppm 1.0e-4

timestep 0.001

fix 1 all shake 0.0001 20 0 b 1 a 1

fix 22 all nve

compute ke all ke/atom
variable wtemp atom c_ke/1.51.6/1.380650310000*3/2
fix 88 all ave/spatial 1 500000 500000 z lower 0.01 v_wtemp file tmp100.profile units reduced
fix 99 all ave/spatial 1 500000 500000 z lower 0.1 v_wtemp file tmp10.profile units reduced
fix 999 all ave/spatial 1 500000 500000 z lower 0.5 v_wtemp file tmp2.profile units reduced

compute myRDF all rdf 200 1 1 1 2 2 1 2 2
fix 220 all ave/time 1 500000 500000 c_myRDF[1] c_myRDF[2] c_myRDF[3] c_myRDF[4] c_myRDF[5] c_myRDF[6] c_myRDF[7] c_myRDF[8] file RDF2.rdf mode vector

compute myMSD all msd
fix 110 all ave/time 1 500000 500000 c_myMSD file MSD.rdf mode vector

thermo 1000
thermo_style custom step temp etotal vol
restart 100000 water.restart

run 4000000

###################OUTPUT for temperature

Dear all

I model water using SPC/E with periodic BC in all directions to obtain RDFs. Using NVT, I obtained equilibrium at 300K (2ns) and then measure RDFs by NVE ensemble for 7ns. However, water temperature was not constant and increased each time step. Being specific, in 7ns temperature increased 300K to 316K. My input and output are simply like following. My theory is there is an energy drift due to the shake algorithm. I believe that during the process of keeping the constrains constant, there is a force exerted by shake algorithm which appears as a heat dissipation. Am I doing something wrong or shake does that?

expecting perfect energy conservation with fix nve is
ignoring the realities of doing floating point math.

but you also have some rather sloppy settings for
PPPM convergence and particularly for SHAKE convergence,
that are very likely to have an impact on energy
conservation.

axel.

Hi,

you also should be aware that the ik-differentiation employed in the
PPPM does not conserve the energy. You will observe a drift even when
specifying a very small accuracy.

Best,

Rolf