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.
http://www.youtube.com/watch?v=S8sEh9RurUM
Thanks in advance
Murat Barisik
########################INPUT
units metal
dimension 3
atom_style full
boundary p p p
neighbor 2.0 bin
bond_style harmonic
angle_style harmonic
read_restart water.restart.2000000
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