Molecular Dynamics Simulation of Water
Molecular Dynamics Water Simulator Test (1-3-95)
Results for Water simulations using the SPC and SPC/E models and with temperature and pressure coupling using the Berendsen coupling technique.
Note: (2/2/2008) This old open source C code is a test of a molecular dynamics simulation of water using the Berendsen Water model. Updated the code (Took out some old archaic code) on 2/2/2008,,, And 13 years later it still runs!
Contents |
Source Code
- Water2a.c — open source.c file (Molecular Dynamics simulation software source)
Notes:
- The ensemble is not defined for Berendsen coupling technique. (labeled below “NPT” for convenience)
- Cv not valid during NPT ensemble, also since not a true NPT ensemble Cp cannot be calculated correctly.
- Velocity AutoCorrelation Coef. (uses Oxygen velocity, not center of mass velocity)
- Virial and kinetic energy, calculated using atomic interactions, instead of center of mass. (This will cause larger fluctuations in pressure [per Berendsen et al, 1984].)
- Pressure fluctuations became unstable when used isothermal compressibility value much above 2.0e-10.
- Property definations:
- KE,PE,E: in Kcal/mole
- P: pressure in Bar
- T: temperature in Kelvin
- V: volume (meters cubed)
- Hbonds: average # of hydrogen bonds per molecule (EPD < -2.8 kcal/mol)
- Density: density (gm/cm cubed)
- Cv: specific heat (constant volume, Joules/mole)
- D: self-diffusion coef. (10-9 metersq/sec)
- <XXXrms>: rms values (fluctuations)
- gOO,gOH: Radial Pair Distribution, positionAngstrom(height)
- EPD: Interaction Energy (Energy Pair Dist) (kcal/mol)
- VACF: Velocity AutoCorrelation function (femptoseconds,fs)
- Runlength: (in picoseconds)
- time step (dt): 5e-016 sec. (for Prevost runs = 2e-15)
- Leonard-jones sigma: 3.16556e-010
- Leonard-jones epslon: 1.07935e-021
- bondlength H-O: 1e-010
- bondangle: 109.47
- virtualbondlen H-H: 1.63298e-010
- rcut: 8.5e-010
- tt: coupling coef for temp .1ps
- pt: coupling coef for press .1ps
- isothermal compressibility: 2.0e-10 [berendsen 1984] uses 4.9e-10]
- To Convert Kcal/mole to Kjoule/mole : kcal/mole x 4.1858 = kj/mole
Results
Simulation results from literature:
Prevost Prevost Prevost
Property SPC(ewald) SPC(mei4) TIPS(ewald) Exp
---------- ---------- ---------- ----------- ---------
Ensemble NVE NVE NVE -
T target - - - -
P target - - - -
RunLength - - - -
<KE> - - - -
<PE> -9.43 -8.99 -8.11 -9.90
<E> - - - -
<P> - - - -
<PV/NKT> 1.3 2.0 .75 -
<T> 299 297 293 300
<V> 6.47e-27
<Hbonds> - - - -
<Density> 1.00 1.00 1.00 -
Cv 77.4 68.5 68.0 75.1
D 4.69 4.70 6.69 2.40
<KErms> - - - -
<PErms> - - - -
<Erms> - - - -
<Prms> - - - -
<Trms> - - - -
<Vrms> - - - -
gOO Zero1 - 2.5 - -
gOO Max1 2.75(2.8) 2.75(2.66) 2.79(2.46) 2.85(2.3)
gOO Zero2 - 3.2 - -
gOO Min1 - 3.8 (.94) - -
gOO Zero3 - 4.1 - -
gOO Max2 - 4.5 (1.04) - -
gOH Zero1 - 1.7 - -
gOH Max1 1.75(1.35) 1.80(1.30) 1.80(1.10) 1.97(1.19)
gOH Zero2 - 1.95 - -
gOH Min1 - 2.3(.22) - -
gOH Zero3 - 3.0 - -
gOH Max2 - 3.3(1.5) - -
EPDmin -2.9 -2.9 - -
EPDmax -4.9 -4.9 - -
VACFmin1 100 100 100 -
VACFmax1 130 130 130 -
VACFzero1 200 200 250 -
VACFmin2 340 340 400 -
Simulation results from simulator:
Property Water-run1 Water-run2 Water-run3 Water-run4 Water-run5 ---------- ---------- ---------- ---------- ---------- ---------- Ensemble NPT NPT NPT NPT NPT T target 300 300 300 300 300 P target 1000 100 10 4 1 RunLength 7 11 12 19.5 12 <KE> 1.78 1.78 1.78 1.78 1.78 <PE> -9.60 -9.55 -9.54 -9.54 -9.55 <E> -7.82 -7.77 -7.76 -7.76 -7.77 <P> 972 94 9.9 1.81 1.8 <PV/NKT> .68 .07 .01 .005 .005 <T> 299.8 299.76 299.64 299.67 299.82 <V> - - 6.42e-27 6.42E-27 6.44e-27 <Hbonds> 2.98 3.16 2.95 2.86 3.04 <Density> 1.04 1.01 1.01 1.01 1.00 Cv NA NA NA NA NA D 4.23 3.98 3.31 3.99 4.13 <KErms> .048 .050 .048 .049 .048 <PErms> .055 .063 .072 .061 .070 <Erms> .028 .039 .053 .036 .051 <Prms> 500 491 495 511 462 <Trms> 8.1 8.4 8.1 8.3 8.1 <Vrms> - - 7.02e-29 8.09e-29 5.45e-29 gOO Zero1 2.6 2.6 2.6 2.6 2.6 gOO Max1 2.76(2.74) 2.76(2.76) 2.76(2.77) 2.76(2.75) 2.76(2.83) gOO Zero2 3.52 3.26 3.26 3.30 3.24 gOO Min1 3.72(.95) 3.6(.918) 3.7(.948) 3.56(.934) 3.62(.933) gOO Zero3 4.84 4.16 4.12 4.22 4.14 gOO Max2 4.84(1.00) 4.44(1.03) 4.48(1.03) 4.62(1.03) 4.50(1.05) gOH Zero1 1.68 1.68 1.68 1.68 1.68 gOH Max1 1.78(1.35) 1.78(1.38) 1.78(1.37) 1.78(1.37) 1.78(1.40) gOH Zero2 1.96 1.92 1.92 1.92 1.92 gOH Min1 2.38(.26) 2.36(.24) 2.42(.25) 2.42(.25) 2.42(.24) gOH Zero3 2.94 2.96 2.96 2.96 2.96 gOH Max2 3.24(1.53) 3.24(1.51) 3.26(1.53) 3.26(1.52) 3.24(1.54) EPDmin -2.9 -2.88 -2.82 -2.84 -2.84 EPDmax -4.86 -4.94 -4.98 -4.92 -4.92 VACFmin1 - - 92.5 92.5 92.5 VACFmax1 - - 130.0 130.0 130.0 VACFzero1 - - 180 177.5 180 VACFmin2 - - 242.5 340.0 242.5
Simulation results from simulator:
Started from a previous NPT run with these targets Property Water-run6 Water-run7 Water-run8 Water-run9 Water-rn10 ---------- ---------- ---------- ---------- ---------- ---------- Ensemble NVE NVE NVE NVE NPT T target 300* 300* 300* 300* 300 P target ? 10* 4* 1* 2770 RunLength 13 25 20 22 22.5 <KE> 1.77 1.84 1.77 1.80 1.78 <PE> -9.51 -9.42 -9.57 -9.49 -9.64 <E> -7.74 -7.57 -7.80 -7.69 -7.86 <P> -101 64.59 112.4 -18.8 2748 <PV/NKT> -.07 .17 .08 -.014 1.81 <T> 297.11 309.25 297.0 302.94 299.67 <V> 6.47e-27 6.42e-27 6.42e-27 6.44e-27 - <Hbonds> 3.07 2.90 3.06 2.79 3.04 <Density> 1.00 1.01 1.01 1.00 1.10 Cv 67.7 63.7 70.1 67.4 NA D 4.07 4.31 3.57 4.48 4.09 <KErms> .057 .057 .058 .058 .049 <PErms> .057 .057 .058 .058 .075 <Erms> .002 .002 .002 .002 .056 <Prms> 775 786 801 873 513 <Trms> 9.6 9.6 9.8 9.8 8.3 <Vrms> 0 0.00 0.0 0.0 - gOO Zero1 2.6 2.6 2.6 2.6 2.6 gOO Max1 2.76(2.77) 2.76(2.68) 2.76(2.78) 2.76(2.74) 2.76(2.66) gOO Zero2 3.28 3.34 3.24 3.28 3.66 gOO Min1 3.56(.928) 3.66(.937) 3.68(.927) 3.68(.938) 4.12(.93) gOO Zero3 4.12 4.36 4.20 4.14 5.5 gOO Max2 4.26(1.02) 4.54(1.02) 4.56(1.03) 4.58(1.03) 6.12(1.03) gOH Zero1 1.7 1.70 1.68 1.68 1.7 gOH Max1 1.80(1.37) 1.78(1.32) 1.78(1.38) 1.78(1.35) 1.78(1.29) gOH Zero2 1.92 1.92 1.92 1.92 1.9 gOH Min1 2.38(.25) 2.40(.24) 2.40(.24) 2.40(.25) 2.36(.27) gOH Zero3 2.96 2.96 2.96 2.96 2.92 gOH Max2 3.26(1.52) 3.26(1.50) 3.26(1.52) 3.26(1.52) 3.24(1.53) EPDmin -3.02 -2.9 -2.82 -2.94 -3.02 EPDmax -4.86 -4.86 -5.04 -4.90 -4.86 VACFmin1 92.5 95.0 92.5 97.5 - VACFmax1 127.5 130.0 130.0 130.0 - VACFzero1 185 185 182.5 187.5 - VACFmin2 377.5 317.5 232.5 312.5 -
Simulation results from simulator:
Property Water-rn11 Water-rn12 Water-rn13 ---------- ---------- ---------- ---------- Ensemble NPT NPT NPT T target 300 300 300 P target 10 4 1 RunLength 33.5 28.5 33.5 <KE> 1.78 1.78 1.78 <PE> -10.68 -10.68 -10.68 <E> -8.90 -8.89 -8.89 <P> 3.74 -0.33 -3.1 <PV/NKT> .007 .003 .001 <T> 299.83 299.79 299.82 <V> 6.29e-27 6.30e-27 6.31e-27 <Hbonds> 3.31 3.35 3.3 <Density> 1.03 1.03 1.02 Cv NA NA NA D 2.71 2.92 2.73 <KErms> .052 .052 .050 <PErms> .076 .086 .080 <Erms> .059 .072 .067 <Prms> 511 523 500 <Trms> 8.7 8.77 8.3 <Vrms> 8.06e-29 6.5e-29 8.0e-29 gOO Zero1 2.58 2.58 2.58 gOO Max1 2.74(2.93) 2.74(2.99) 2.74(3.00) gOO Zero2 3.10 3.10 3.10 gOO Min1 3.36(.872) 3.32(.874) 3.30(.860) gOO Zero3 4.08 3.98 4.00 gOO Max2 4.32(1.05) 4.50(1.07) 4.38(1.07) gOH Zero1 1.66 1.66 1.66 gOH Max1 1.76(1.51) 1.76(1.54) 1.76(1.54) gOH Zero2 1.90 1.90 1.90 gOH Min1 2.38(.198) 2.40(.200) 2.40(.198) gOH Zero3 2.96 2.96 2.96 gOH Max2 3.24(1.55) 3.24(1.57) 3.22(1.57) EPDmin -2.92 -3.06 -2.96 EPDmax -5.46 -5.48 -5.38 VACFmin1 85.0 85.0 85.0 VACFmax1 122.5 122.5 125.0 VACFzero1 162.5 162.5 165.0 VACFmin2 297.5 292.5 292.5
Simulation results from simulator:
Started from previous a NPT run with these targets
spc/e spc/e spc/e
Property Water-rn14 Water-rn15 Water-rn16 Water-rn17 Water-rn18
---------- ---------- ---------- ---------- ---------- ----------
Ensemble NVE NVE NVE NPT NVE
T target 300* 300* 300* 300 300*
P target 10* 4* 1* 10 10*
RunLength 6 11 6 80 80
<KE> 1.97 1.80 1.79 1.78 1.99
<PE> -10.31 -10.62 -10.63 -9.53 -9.13
<E> -8.34 -8.82 -8.84 -7.74 -7.13
<P> 312 142.53 -22.47 10.7 396.7
<PV/NKT> .20 .10 -.02 .01 .26
<T> 331 302.2 300.6 299.66 334.99
<V> 6.29e-27 6.30e-27 6.31e-27 6.45e-27 6.45e-27
<Hbonds> 3.13 3.46 3.3 3.04 2.81
<Density> 1.03 1.03 1.02 1.00 1.00
Cv 70.3 72.7 71.4 NA 65.4
D 2.81 3.04 2.68 4.42 5.55
<KErms> .065 .060 .059 .050 .063
<PErms> .065 .060 .059 .070 .063
<Erms> .002 .002 .002 .050 .002
<Prms> 1027 884 919 493 865
<Trms> 10.9 10.1 10.0 8.4 10.6
<Vrms> 0 0 0 6.45e-29 0.00
gOO Zero1 2.6 2.58 2.58 2.6 2.6
gOO Max1 2.74(2.71) 2.74(2.97) 2.74(2.95) 2.76(2.80) 2.76(2.52)
gOO Zero2 3.24 3.12 3.14 3.28 3.48
gOO Min1 3.54(.912) 3.44(.883) 3.34(.891) 3.6(.953) 3.86(.950)
gOO Zero3 4.08 4.04 3.90 4.08 6.10
gOO Max2 4.40(1.04) 4.36(1.07) 4.22(1.05) 4.50(1.04) 6.48(1.00)
gOH Zero1 1.68 1.66 1.66 1.68 1.70
gOH Max1 1.76(1.36) 1.76(1.52) 1.76(1.52) 1.78(1.39) 1.80(1.20)
gOH Zero2 1.90 1.90 1.90 1.92 1.90
gOH Min1 2.38(.25) 2.38(.20) 2.36(.20) 2.40(.25) 2.40(.30)
gOH Zero3 2.94 2.96 2.96 2.96 2.96
gOH Max2 3.26(1.51) 3.24(1.56) 3.24(1.55) 3.26(1.54) 3.26(1.47)
EPDmin -3.02 -3.06 -2.92 -2.88 -3.04
EPDmax -5.24 -5.56 -5.44 -4.92 -4.72
VACFmin1 90.0 85.0 92.5 92.5 97.5
VACFmax1 130.0 122.5 127.5 130.0 127.5
VACFzero1 172.5 165.0 167.5 182.5 200.0
VACFmin2 297.5 307.5 212.5 320.0 325.0
References:
- Prevost, et al. – Mol Physics 1990, 71(3): 587-603 (MEI4, F1)
- Berendsen, et al. – J phys chem, 1987, 91: 6269-6271 (SPC/E)
- Brooks et al – J chem phys 1985, 83(11): 5897 (MEI4)
- Berendsen, et al. – J chem phys, 1984, 81(8): 3684 (Temp/Press Bath)
- Jorgensen, et al. – J phys chem, 1983, 79(2): 926 (TIPS water)
- Berendsen, et al. – Intermolecular Forces 1981, 331-342 (SPC water)
- Allen-Tildesley – computer simulation of liquids