Liquid water MD simulation

A lot of stories on water:
Please refer to the excellent web-site
Martin Chaplin’s “Water Structure and Science”
http://www.lsbu.ac.uk/water/index2.html

ModelTypeσ/Åε / kJ
mol-1

l1 Ål2 Å

q1 (e)

SPC/Ea

3.166

0.650

1.0000-+0.423 -0.847
109.47
86-

TIP3P

a

3.1506
0.6364 0.9572
1-+0.417 -0.834
104.52
00-

TIP4P

c

3.15360.6480 0.9572
50.15+0.520 -1.040
104.52 52.26
00TIP4PEwc

3.1643 0.6809
0.95725460.125+0.524 -1.048
104.52 52.26
2244

TIP5PEwd3.097

0.70+0.241 -0.241
109.4
104.520070.7448 0.9572

q2 (e)

θ°φ°

Water model potential

SPC/E: H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The missing term in effective pair potentials, J. Phys. Chem. 91 (1987) 6269-6271.
TIP3P: W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J.
Chem. Phys. 79 (1983) 926-935
TIP4P: W. L. Jorgensen and J. D. Madura, Temperature and size dependence for monte carlo simulations of TIP4P water, Mol. Phys. 56 (1985) 1381-1392.
TIP4P-Ew: H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura and T. Head-Gordon, Development of an improved four-site water model
for biomolecular simulations: TIP4P-Ew, J. Chem. Phys. 120 (2004) 9665-9678.
TIP5P-Ew: S. W. Rick, A reoptimization of the five-site water potential (TIP5P) for use with Ewald sums, J. Chem. Phys. 120 (2004) 6085-6093.
TIP5P: M. W. Mahoney and W. L. Jorgensen, A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potentialfunctions, J. Chem. Phys. 112 (2000) 8910-8922.

Calculated results：classical MD
Average
Self diffusion, configurational
10-5 cm2/s
energy, kJ
mol-1

Densitymaximum,°C

Expansioncoefficient,10-4 °C-1-41.5+3.984

2.53

2.49

-41.5

-38

5.14

82

5.19

-41.1

-91

9.2

2.18

53

3.29-41.8-25

4.4

TIP4P-Ew

2.32

62.9

2.4-46.5+1

3.1

TIP5P-Ew

2.29

92

2.8-+84.9ModelDipole
momentDielectricconstant

Exp.

2.95

78.4

2.30

SPC/E

2.35

71

TIP3P

2.35

TIP4PAll the data is at 25°C and 1 atm==>The popular models SPC, SPC/E, TIP3P and TIP4P produce poor agreement with water's
melting point (giving melting points of 190 K, 215 K, 146 K and 232 K respectively)==>The popular TIP4P model underestimates the tetrahedrality of the water molecule's
environment, which explains its poor estimate of the dielectric constant. It is, however, remarkably
good at qualitatively describing water's phase diagram. The SPC/E and TIP4P models are reported asfailing to properly describe the experimental O···O radial distribution function. The TIP3P and SPC
show particularly poor agreement, the TIP4P, SPC/E and PPC show improved agreement.

modelpotential
SCORE==> We will use
SPC/EandTIP4P/Ewmodels.

flexible model: time consuming and not improved from rigid model
Polarizable model: so far not significantly improved, but maybe important forion solvation, for interface, and for ILs

ALS2000 and Soper 2000

ALS:X-ray 2000
Narten: X-ray 1972
Soper: Neutron 97
Soper 2000:reanalysis of1997 data

good agreement !!
SPC/E and ALS2000

Radial distributionfunction(RDF)from exp.
and MD.

Classical model
potentials canexplain
the Exp. resultsvery well.

Exp. results arealso improving.

Radial distributionfunction from exp.
and MD.

first-pinciples
MD is getting
better.

Here we use one of the water model SPC/E or TIP4P-Ew

1) L-J interaction is oxygen-oxygen only.
2) rigid molecule is assumed. ==> Cal. become fasterO-H distance is const.
angle H-O-H is const.
==> Cal. become faster==> but neglect the vibrational degree of freedom.
==> we should select how we impose the constraints onbond length and bond angle.
3) no electrostatic interaction between the intramolecularatoms.

Z

initial configuration: ice?
Cubic Ice

fcc + half the tetrahedral sites are filled.

ρ / g cm-3 = 0.996516
H2O all site(fcc+t-site)= 512
liquid water all atom number= 1536
ax,ay,az,volume= 24.86, 24.86, 24.86 Å 15370.0X

Z

SPC/E model-0.8476 e
1.0000 Å
109.47°+0.4238 er(OH)=1.0000 Åθ(HOH)=109.47
are fixed.

intermolecular interactionz(O)=-0.8476 ez(H)=0.4238 eε(O) = 0.650 kJ mol-1
σ(O) = 3.166 ÅCONFIGH2O liquid SPCE water at 300 K
2
1
24.8632886676.0000000000.0000000000.0000000000
24.8632886676.0000000000.0000000000.0000000000
24.8632886676
OW 1-12.4316443338-12.4316443338-12.4316443338-.1157118156-.1006883860.0252069353.0000000000.0000000000.0000000000
HW 2-11.8542984133-11.8542984133-11.8542853673-.1157118156-.1006883860.0252069353.0000000000.0000000000.0000000000
HW 3-13.0089902543-13.0089902543-11.8542853673-.1157118156-.1006883860.0252069353.0000000000.0000000000.0000000000
OW 4-10.8776887921-10.8776887921-10.8776887921-.0944489040.0834090526-.0070384049.0000000000.0000000000.0000000000
HW 5-10.3003428716-11.4550347126-10.3003298256-.0944489040.0834090526-.0070384049.0000000000.0000000000.0000000000
HW 6-11.4550347126-10.3003428716-10.3003298256-.0944489040.0834090526-.0070384049.0000000000.0000000000.0000000000...
OW 1534
7.7697777086
10.8776887921
10.8776887921.0179460404.0359122702.0141856681.0000000000.0000000000.0000000000
HW 1535
8.3471236291
10.3003428716
11.4550477586.0179460404.0359122702.0141856681.0000000000.0000000000.0000000000
HW 15367.1924317881
11.4550347126
11.4550477586.0179460404.0359122702.0141856681.0000000000.0000000000.0000000000

FIELD
H2O SPC/Eunits kJ
molecules 1
SPCE waternummols 512
atoms 3OW
15.9996 -0.8476
HW
1.00080.4238HW
1.00080.4238constraints 31 2 1.0000
1 3 1.0000
2 3 1.63298
finishvdW 3OWOW lj 0.6507346 3.165492
OW
HW lj 0.0
3.165492
HW
HW lj 0.0
0.0
close

SPC/E modelCONTROL NVE

CONTROL NPT 300 K 1 atom
SPC/E model of waterSPC/E model of water
integrator velocity verlet
temperature
300.00ensemble nve
steps
150000equilibration 100000
multiple step
1scale10
print50
stack100stats
10
rdf
1
traj 1 100 0
timestep0.001
cutoff9.50
delr width
0.5000ewald precision 1d-6
shake tolerance 1.0E-5
quaternion tolerance 1.0E-5
print rdf
job timeclose timefinish150000.00
500.00integrator velocity verlet
temperature
300.00ensemble npt hoover 0.1 0.2
steps
150000equilibration 100000
multiple step
1
scale
10
print50
stack100stats
10
rdf
1
traj 1 100 0timestep
0.001
cutoff9.50
delr width0.5000
ewald precision 1d-6
shake tolerance 1.0E-5
quaternion tolerance 1.0E-5
print rdf
pres 0.001
job timeclose time

150000.00
500.00finish

SPC/E modelNVE : microcanonical ensemble==> total energy conserved==> directly linked to mechanics==> kBT = kinetic energy
Temperature fluctuates
but the termal average should be the sameNVT : canonical ensemble==> T keeps constant andthe heat exchange (or momentum exchange)
with the virtual reservoir (thermostat)===> some tricks are needed.
Nose-Hoover algorithm is only well defined!!
Text
NVT

NPT ensemble==> T and P keeps constant andthe heat exchange and cell size change
dynamically.
===> Hoover barostat method is only well defined!!
Text

Pext
NPT

SPC/E model:
NVE ensemble
results

time elapsed since job start =temp 2.9618E+02 K r.m.s. fluctuation 6.3217E+00 K

27144.723 seconds

3.0
SPC/E model of liquid water: NVE ensemble

Approximate 3D Diffusion coefficients (10^-9 m^2 / s)D2.6341E+00
2.6341E+00

Average pressure tensor-8.2421E-02 2.9661E-02 5.1028E-03
2.9661E-02 -3.2231E-02 2.5959E-02
5.1028E-03 2.5959E-02 -8.3213E-02
trace/3. -6.5955E-02 katom

2.5radial distribution function

atom
OW
HW

g(O-O)
g(O-H)
g(H-H)
2.0

1.5

1.0

0.5

0.0
24

6
distance / Angstrom

8

SPC/E model:
NPT ensemble
results
T = 3.0011E+02 K r.m.s. fluc. 7.1254E+00 K
Approximate 3D Diffusion coefficients (10^-9m^2 / s)
atom
OW
HWD
2.4579E+00
2.4818E+00

movie

Average pressure 1.2925E-03 katom
Average pressure tensor
r.m.s. fluctuations-3.3618E-02 1.3665E-02 3.6790E-02
6.2210E-01 5.0643E-01 5.1488E-01
1.3665E-02 -1.0684E-02 5.6640E-03
5.0643E-01 6.1777E-01 5.0854E-01
3.6790E-02 5.6640E-03 4.8179E-02
5.1488E-01 5.0854E-01 6.2829E-01
trace/3. 1.2925E-03
Average cell vectors
24.8372539181.0000000000.0000000000
24.8372539181.0000000000.0000000000.0000000000.0000000000
24.8372539181

r.m.s. fluctuations
7.1583E-02 0.0000E+00 0.0000E+00
0.0000E+00 7.1583E-02 0.0000E+00
0.0000E+00 0.0000E+00 7.1583E-02

time elapsed since job start =24121.500 seconds

SPC/E model:
NPT ensemble
results

cf. exp. results

3.0

radial distribution function

2.5

SPC/E liquid water NPT ensembleg( O-O )g (O-H )g( H-H )

2.0

1.5

1.0

0.5

0.0
24
6
distance / Angstrom8

TIP4P/Ew modelr(OH)=0.9572 År(OM)=0.125 Åθ(HOH)=104.52
θ(MOH)=52.26
are fixed.

OActiveHMHM: dummy charge sitez(M)=-1.04844 ez(H)=0.52422 eε(O) = 0.680946 kJ mol-1
σ(O) = 3.16435 ÅYmass of M is zero.

ZXCONFIG
TIP4P H2O liquid water 300 K
2124.8632889271.0000000000.0000000000.0000000000
24.8632889271.0000000000.0000000000.0000000000
24.8632889271
OW 1-12.4316444636-12.4316444636-12.4316444636-.1157118180-.1006883881.0252069358.0000000000.0000000000.0000000000
HW 2-11.8963997541-11.8963997541-11.8457621869-.1157118180-.1006883881.0252069358.0000000000.0000000000.0000000000
HW 3-12.9668891730-12.9668891730-11.8457621869-.1157118180-.1006883881.0252069358.0000000000.0000000000.0000000000
MW 4-12.4316444636-12.4316444636-12.3066444636-.1157118180-.1006883881.0252069358.0000000000.0000000000.0000000000
OW 5-10.8776889056-10.8776889056-10.8776889056-.0944489059.0834090543-.0070384050.0000000000.0000000000.0000000000
HW 6-10.3424441962-11.4129336150-10.2918066290-.0944489059.0834090543-.0070384050.0000000000.0000000000.0000000000
HW 7-11.4129336150-10.3424441962-10.2918066290-.0944489059.0834090543-.0070384050.0000000000.0000000000.0000000000
MW 8-10.8776889056-10.8776889056-10.7526889056-.0944489059.0834090543-.0070384050.0000000000.0000000000.0000000000

FIELD
H2O TIP4P
units kJ
molecules 1
TIP4P water
nummols 512
atoms 4OW
15.9996 0.0
HW1.00080.52422
HW1.0008
0.52422
MW0.0-1.04844
constraints 3
1 2 0.95720000
1 3 0.95720000
2 3 1.51390065
rigid 12 1 4finish
vdW 6OWOW lj 0.680946 3.16435
OW
HW lj 0.0
0.0
HW
HW lj 0.0
0.0
OW
MW lj 0.0
0.0
HW
MW lj 0.00.0
MW
MW lj 0.0
0.0
close

CONTROL NVE

CONTROL NPT 300 K 1 atom

TIP4P model of water

TIP4P model of water

integrator velocity verlet
temperature
300.00ensemble nve

integrator velocity verlettemperature
300.00ensemble npt hoover 0.1 0.2

steps
150000equilibration 100000
multiple step
1scale
10
print
50
stack
100
stats
10
rdf
1
traj 1 100 0
timestep0.001
cutoff9.50
delr width
0.5000ewald precision 1d-6
shake tolerance 1.0E-5
quaternion tolerance 1.0E-5
print rdf
job timeclose timefinish300000.00
500.00steps
150000equilibration 100000
multiple step
1
scale
10
print
50
stack
100
stats
10
rdf
1
traj 1 100 0
timestep0.001
cutoff9.50
delr width
0.5000ewald precision 1d-6
shake tolerance 1.0E-5
quaternion tolerance 1.0E-5
print rdf
pres 0.001
job timeclose timefinish300000.00
500.00

TIP4P/Ew model: NVE ensemble results
temp_tot 3.0270E+02 K r.m.s. fluctn.

movie

6.4702E+00 K

Approximate 3D Diffusion coefficients (10^-9 m^2 / s)
OW
HWMW

atom
D
2.6537E+00
2.6627E+00
2.6538E+00

Average pressure tensor

r.m.s. fluctuations

4.8443E-02 6.8694E-03 1.4413E-02
6.8694E-03 3.8647E-02 2.9359E-02
1.4413E-02 2.9359E-02 4.2101E-02

7.6492E-01 5.0866E-01 5.1297E-01
5.0866E-01 7.9529E-01 5.0658E-01
5.1297E-01 5.0658E-01 8.0716E-01

trace/3. 4.3064E-02
3.0
TIP4P/Ew NVE ensembleg(O-O)
g(O-H)
g(H-H)

2.5

RDF

2.0

1.5

1.0

0.5

0.0
24

6

distance / Angstrom

8time elapsed since job start =34789.530 seconds

TIP4P/Ew model: NPT ensemble results
Temp. 2.9996E+02 K +- 7.7193E+00

Approximate 3D Diffusion coefficients(10^-9 m^2 / s)
atomOW
HWMWmovieD
2.5677E+00
2.5421E+00
2.5601E+00

Average pressure tensor
2.3043E-02 -3.2620E-02 1.5465E-02-3.2620E-02 4.4084E-03 -1.0378E-02
1.5465E-02 -1.0378E-02 -1.3775E-02

r.m.s. fluctuations
8.9626E-01 5.0277E-01 4.8920E-01
5.0277E-01 8.5233E-01 4.9858E-01
4.8920E-01 4.9858E-01 8.6417E-01

trace/3. 4.5587E-03

Average pressure
4.5587E-03 katom
Average cell vectors
24.8869811588.0000000000.0000000000.0000000000
24.8869811588.0000000000

r.m.s. fluctuations.0000000000.0000000000
24.8869811588

1.0249E-01 0.0000E+00 0.0000E+00
0.0000E+00 1.0249E-01 0.0000E+00
0.0000E+00 0.0000E+00 1.0249E-01

time elapsed since job start =34773.057 seconds

TIP4P/Ew model: NPT ensemble results
cf. exp. results
3.0

2.5

TIP4PEw NPT ensembleg(O-O)
g(O-H)
g(H-H)

RDF

2.01.51.00.50.0
246

distance / Angstrom8