CC BY
12
УДК 544.272, 539.6, 532.74
An Efficient Water Model for Large Scale Molecular Dynamics Simulations
Viktor E. Zalizniak*
Institute of Mathematics and Computer Science Siberian Federal University Svobodny, 79, Krasnoyarsk, 660041
Russia
Institute of Computational Modelling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036
Russia
Received 10.06.2015, received in revised form 24.08.2015, accepted 01.10.2015 The development of simple and efficient model that correctly represent the important features of water is essential to overcome the limitations in time scale and system size currently encountered in atomistic molecular dynamics simulations. The proposed one site model includes Lennard-Jones interaction and the angular averaged dipole-dipole interaction. Experimental data of liquid water at various temperatures are used for parametrization of the model. The values of density were chosen as primary target properties. These properties cover a temperature range from 300 to 350 K and pressures up to 10.1 MPa. The model properties are compared with those obtained from experiment and from general purpose TIP4P/2005 model. The comparison shows that all chosen properties are quite well reproduced by the proposed model. Computational scheme that is used in simulations is also presented. The proposed water model reproduces the key characteristics of liquid water while being computationally considerably more efficient than standard multi-site atomistic water models. The model is for use on large scale simulations of the fluid behavior in nanosized structures.
Keywords: water model, molecular dynamics simulation. DOI: 10.17516/1997-1397-2015-8-4-487-496
Introduction
Due to its role in biological and technological processes, there is a great interest in an accurate knowledge of the molecular interactions in water. For this reason, a large number of potential models for water have been proposed in the past. One can recommend the reviews by Guillot [1] and Finney [2] for a critical discussion of the results for various water models. One can also recommend the Chaplin website [3] containing up to date information on water structure and behavior.
Theoretical studies of water usually employ simple empirical models. Computer simulations can then predict properties of water which are to be compared with measurements. In the commonly used interaction-site models a rigid geometry is adopted and the interaction between water molecules is calculated via pair wise summation. These models typically use three to six
* vzalizniak@sfu-kras.ru © Siberian Federal University. All rights reserved
interaction sites. The SPC/E [4] and TIP3P [5] models use three sites for the electrostatic interactions. The partial positive charges on the hydrogen atoms are balanced by an appropriately negative charge located on the oxygen atom. The intermolecular interaction between two water molecules is computed using a Lennard-Jones type potential with just a single interaction point per molecule centered on the oxygen atom. No van der Waals interactions involving the hydrogen atoms are calculated. The four-site models such as OPC [6] and TIP4P/2005 [7] model shift the negative charge from the oxygen atom to a point along the bisector of the HOH angle towards the hydrogen atoms. The most commonly used five-site models are the ST2 potential [8] and TIP5P [9] model. Here, charges are placed on the hydrogen atoms and on the two lone-pair sites on the oxygen.
The computational cost of a water simulation increases with the number of interaction sites in the water model. The CPU time is approximately proportional to the number of interatomic distances that need to be computed. In addition, there is degradation in efficiency due to the longrange Coulombic calculation performed via the particle-particle, particle-mesh (PPPM) method or the related particle-mesh Ewald (PME) method. These are FFT based methods with an O (N log N) complexity in the number of atoms, so their cost grows as the system size increases.
The system size and the time scale of atomistic simulations are, however, limited to tens of nanometers and nanoseconds by computational cost. A possible remedy is coarse graining (CG) of "unimportant" degrees of freedom, in which a group of atoms is represented by a single bead or particle. This may considerably reduce the number of particle-particle interactions to be calculated [10] (see also the references cited therein). However, CG models are less accurate in comparison with multi-site models in predicting water properties.
With the development of novel micro- and nanofluidic devices, applications of flows are in increasing demand. Computer simulations of such devices require massive amounts of supercomputer time as systems such as water trapped inside nanopores are computationally very demanding. The limitation of the present water models lies mainly in the computer power. The number of molecules, which can be handled, is still not sufficient to enhance engineering design. Hence, there is the need to develop water models, which will be able efficiently simulate physical phenomena in micro- and nano devices.
The aim of this paper is to present simple, computationally efficient specialized model for liquid water for a given range of temperatures and pressures.
1. A potential model for water
We start by assuming that the total energy of interaction Uint (r) between two water molecules is divided into different terms with a physical meaning according to
Uint (r) = Urep (r) + Udis (r) + Udd (r) + Udid (r) , (1)
where Urep (r) is energy of repulsion, Udis (r) is energy of dispersion interaction, Udd (r) is energy of dipole-dipole interaction and Udid (r) is energy of dipole-induced dipole interaction.
The first two terms in (1) are represented using a Lennard-Jones type potential of the form
/ d\ 12 / d\ ^ Urep M + Udis (r)= £1^- £2^r) ,
where £1, £2 and d are adjustable parameters.
Let us consider the interaction of two dipoles, which in vector notation can be written as
U« M = - . (2)
where pn = pen and pm = pem are dipole moment vectors, r = rn — rm is the vector connecting the two water molecules, r = |r| and p is the value of a dipole moment (for water p = 0.388e • A [11]). Equation (2) can also be written in terms of three angles: and y>2, are the angles between the dipole vectors and the vector r, and <^3 is the dihedral vector formed by
the two dipole moment vectors and the vector r. Then we obtain the following expression
p2
Udd (r) = --3 (sin (^1) sin (^2 )cos(^) — 2 cos (^1) cos (^2)).
Let us now calculate the angular averaged Boltzmann weighted energy of interaction between two dipoles. That is, keep r fixed while integrating over the relevant angles
n n n
///Udd (r)exp(-¡Udd (r)) sin (^1) sin (^2) dyid^d^ (Udd (r)> = -, (3)
/ //exp(-^Udd (r)) sin (y>i) sin (^2) d<pid<p2d^3 ooo
where 3 = 1/kBT with kB being the Boltzmann constant. Let us introduce the dimensionless parameter
p2 d
4neod3kBT \r /
The integral (3) has no analytical solution so it numerically calculated for various values of parameter x in the range 0 ^ x ^ 3. The results are shown in Fig. 1. The numerical data are approximated by the piecewise third order polynomials as follows
( 0.0965x3 - 0.71x2, 0 < x < 1, (Udd (r)> = kBT 2 (4)
[ 0.0907x3 - 0.6925x2 - 0.0175x + 0.0058, 1 <x < 3.
All molecules have an electronic polarizability. This means that a molecule with a permanent dipole moment, like water, will polarize a nearby molecule. In the same way as above this interaction can be thermally averaged and the resulting interaction keeping only the leading term reads
ap2 f dN
(Udid (r)> = -2-
CO
AneodP \ r
where a is the average molecular polarizability (for water a = 1.45A3 [12]).
Finally, we arrive to the following expression for the averaged interaction energy between two water molecules (ADDI model — Averaged Dipole-Dipole Interaction model):
Uint (r)> = r) - ^d) + (Udd (r)> + (Udid (r)> .
This averaging procedure is quite justified. There are polar fluctuations due to reorientation of the dipolar molecules in liquid water. The relaxation time for the transverse fluctuations is about 10 ps [13], but molecular dynamics simulations are usually performed with the time step of about 1 fs.
-4-1-1-1-1-1-
0 0.5 1 1.5 2 2.5 3
X
Fig. 1. Averaged energy of dipole-dipole interaction versus parameter x. The energy is expressed in terms of kBT. Circles mark the calculated values, solid line is approximation (4)
In practical applications, it is desirable to employ a switching function in order to terminate the potential and forces smoothly at the cut off distance because the energy conservation is sensitive to the truncation of the force field. For this purpose, a simple polynomial switching function fc (r) can be applied to the potential in a region just below the cut off distance rc:
1, r < Tsw
fc (r) = ci (r*)3 + C2 (r*)4 + C3 (r*)5 , rsw < r < rc k 0, r > rc
r
r* =--1,
rc
where rsw is the distance at which the switching function is applied. The values of coefficients c1, c2, c3 follows from the conditions
fc (rsw ) = 1, f (rsw ) = 0, df (rsw ) = 0.
2. Simulation details
A cubic box with 1000 water molecules was used together with periodic boundary conditions. All simulations were performed under isothermic-isobaric (NPT) conditions. The equilibration time for each simulation was 0.6 ns. The production simulations were 0.6 ns long. The fitting procedure was performed using a cutoff distance rc = 10 A.
For the integration of the equations of motion kinetic energy conserving integrator [14] with
the pressure coupling algorithm is employed:
fc + l/2 _ + i vfc rn _ rn + 2 Tk Vn,
k , Tk Vk +--
n
m
F (rn+1/2)
rk+1/2 + i Tk V*
rn + 2 Tk Vn,
V
r
n _!,..., N; k _0,1,...,
1 + akn K)
vn+1 _ (1 + akn (v*J) v*n, (1 + AK (vk)/Kn (vn))1/2 , if AK (vk)/Kn (vkn) > -1
1, otherwise,
rkn+1 = rn (6fc)1/3 , bk = (1 - c (Po - P (tk)))1/2 ,
where rn is the particle position, vn is the particle velocity, m is the particle mass and v = (v1,..., vN). Here AK (vk) = (K0 — K (p^) /N is the kinetic energy difference per particle, where K (vk) is the kinetic energy of the system, K0 = 1.5NfcB T0 is the prescribed kinetic energy and Kn (v^) is the particle kinetic energy. An compressibility parameter c = 16.02 (eV/A3) 1 was used in the pressure coupling. This computational scheme admits relatively large time steps and conserves the total energy of the system over the interval of simulation for NVE conditions. The integration time step was 10.19 fs for all simulations.
3. Results of fitting
It is well known that it is not possible to fit the overall water properties with a single set of parameters. Our aim is to develop simple specialized potential for a given range of temperatures and pressures. The set of properties to be fitted consists of the densities of liquid water at various conditions and configurational energy at 300 K and 0.1 MPa.
There are three unknown parameters to determine, namely, e1, e2 and d. The values of the optimized parameters for the ADDI model are given in Tab. 1.
Table 1. Optimized parameters for the ADDI model (P = 0.1 MP)
£i, eV £2, eV d (300 K), A d (310 K), A d (330 K), A d (350 K), A
0.156 0.084 3.151 3.14 3.0925 3.084
The dependence of parameter d on pressure can be approximated as
d (T, P) = d (T) — 3 • 10"5 ^P0^ ^ , P0 = 0.1 MPa.
4. Results of simulations
In the following section we describe the results of simulations of liquid water with ADDI model in comparison with available experimental data. A summary of the properties at T = 300 K and P = 0.1 MPa is presented in Tab. 2. The rms fluctuations for density p, potential energy Epot, temperature T and pressure P are 2 kg/m3, 0.001 eV, 0.01 K and 10 MPa, respectively. Water properties obtained with the use of general purpose TIP4P/2005 model are also included in Tab. 2 for comparison. The dependence of the density with temperature is plotted in Fig. 2 together with the experimental data. The agreement is good in the temperature between 300 K and 350 K.
Table 2. Computed properties for liquid water at T = 300 K and P = 0.1 MPa
Property ADDI model TIP4P/2005 [7] Experiment
p, kg/m3 996.7 996.5 996.56 [15]
Epot, eV -0.420 - -0.415 [16]
, 10-4 K-1 2.31 2.8 2.57 [17]
kr, 10-10 Pa-1 4.4 4.65 4.58 [17]
Vrms , m/s 644 - - 640 [18]
D, 10-9 m2/s 2.49 2.08 2.3 [17]
10-3 Pa-s 0.47 - 0.85 [17]
1000
9701-i-i-i-i-
300 310 320 330 340 350
T, K
Fig. 2. Densities of the ADDI model (full line) at P = 0.1 MPa compared to experimental data (circles)
The thermal expansion coefficient aV of the ADDI model was obtained using the finite-difference expression
/d^A _ fin(p2/piA
\dT)p ~ \ T2 - Ti )v ■
v \3Tjp v t2 — t1
The isothermal compressibility kT of the ADDI model was calculated using the finite-
difference expression
* = -1 m „ Çg-M) .
v\dPjT v P2 — P1
Coefficients of self-diffusion D and viscosity n were estimated with the use of the Enskog theory [19] as
D = Do
g (deff )'
1
M = Mo ( r7d-T + 0.8b + 0.761g (deff ) b^ ,
,g (deff )
b = 2nndff /3,
where g (deff ) is the radial distribution function at contact, n is the particle number density and deff is the effective diameter of water molecule. The subscript 0 refers to the low-density coefficients. They are given by the kinetic theory [20] as
d0 = 1 V Mo = 0i9Mp
3 n%/2nd2ff nV2nd2ff
where (v) is the average particle velocity and p is the fluid density. The radial distribution function at contact can be calculated from the Carnahan-Starling equation [21]
g (deff ) = 1 053 ' n = nnd3ff/6.
(1 - n)
The effective diameter of water molecule can be estimated from the following equation
(Uint (deff )) = 1 kBT.
For T = 300 K deff = 3.07 A and deff = 2.993 A for T = 350 K. The root mean square velocity was used instead of (v) in estimation of transport coefficients.
The pair distribution function (PDF) of the ADDI model is shown in Fig. 3. The computed PDF has a broad first maximum at 3.25 A. The PDF obtained in experiment shows a first peak positioned at 2.80 A [22].
Conclusions
In this paper the results for a new simple potential model for liquid water is presented. The calculated properties include a number of thermodynamic and transport properties of liquid water. The model gives good performance for this variety of properties and thermodynamic conditions. In this model water molecules cannot form hydrogen bonds, so the model is not in a position to correctly describe structural properties of water. The proposed model is a one-site pair wise additive potential. In addition, there is no need to include rotational dynamics of molecules and the long-range Coulombic calculation. All these features result in significant increase of
4
r,A
Fig. 3. Pair distribution function from simulations with the ADDI model at T = 300 K
the computational speed in computer simulations over multi-site models. The proposed model is intended for use on large scale simulations of the fluid behavior in nanosized structures. In a next step, the model will be applied to simulation of aqueous electrolyte solutions. In particular, the model is intended for simulating electrokinetic fluid flow through active nanomembranes [23].
This research was funded by project 15-19-10017 of the Russian Science Foundation. Author would like to thank O. A. Zolotov for stimulating discussions and technical help.
References
[1] B.Guillot, A reappraisal of what we have learnt during three decades of computer simulations on water, J. Mol. Liq, 101(2002), 219-260.
[2] J.L.Finney, The water molecule and its interactions: the interaction between theory, modelling, and experiment, J. Mol. Liq., 90(2001), 303-312.
[3] M.Chaplin, Water structure and science (http://www.lsbu.ac.uk/water).
[4] H.J.C.Berendsen, J.R.Grigera, T.P.Straatsma, The missing term in effective pair potentials, J. Phys. Chem., 91(1987), 6269-6271.
[5] W.L.Jorgensen, J.Chandrasekhar, J.D.Madura, R.W.Impey, M.L.Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 79(1983), 926-935.
[6] S.Izadi, R.Anandakrishnan, A.V.Onufriev, Building water models: A different approach, J. Phys. Chem. Lett, 5(2014), 3863-3871.
[7] J.L.F.Abascal, C.Vega, A general purpose model for the condensed phases of water: TIP4P/2005, J. Chem. Phys., 123(2005), 234505.
[8] F.H.Stillinger, A.Rahman, Improved simulation of liquid water by molecular dynamics, J. Chem. Phys., 60(1974), 1545-1557.
[9] M.W.Mahoney, W.L.Jorgensen, A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions, J. Chem. Phys., 112(2000), 8910-8922.
10] S.Riniker, W.F.van Gunsteren, A simple, efficient polarizable coarse-grained water model for molecular dynamics simulations, J. Chem. Phys., 134(2011), 084110.
11] D.R.Lide, Ed., CRC Handbook of Chemistry and Physics, 75th Ed., CRC Press, Boca Raton, FL, 1994.
12] K.J.Miller, Additivity methods in molecular polarizability, J. Am. Chem. Soc., 112(1990), 8533.
13] D.P.Shelton, Slow polarization relaxation in water observed by hyper-Rayleigh scattering, Phys. Rev. B, 72(2005), 020201(R).
14] O.A.Zolotov, V.E.Zalizniak, Accurate energy conservation in molecular dynamics simulation, Nanosystems: Physics, Chemistry, Mathematics, 4(2013), 657-669.
15] NIST Chemistry WebBook. NIST Standard Reference Database Number 69 (http://webbook.nist.gov/chemistry/).
16] M.W.Mahoney, W.L.Jorgensen, A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions, J. Chem. Phys., 112(2000), 8910-8922.
17] A.P.Kunz, W.F.van Gunsteren, Development of a Nonlinear Classical Polarization Model for Liquid Water and Aqueous Solutions: COS/D, J. Phys. Chem. A, 113(2009), 11570.
18] M.Kinoshita, Roles of translational motion of water molecules in sustaining life, Fronteirs Biosci, 14(2009), 3419-3454.
19] D.Enskog, Der warmeleitung, reibung und selbstdiffusion in gewissen verdichteten gasen und flussigkeiten, Kungliga Svenska Vetenskapsakademiens Handlingar, Band 63(1922).
20] I.E.Irodov, Physics of Macro-systems. The Fundamental Laws, 4th ed., BKL Publishers, Moscow, 2010.
21] N.F.Carnahan, K.E.Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys., 51(1969), 635-636.
22] V.Petkov, Y.Ren, M.Suchomel, Molecular arrangement in water: random but not quite, J. Phys.: Condens. Matter, 24(2012), 155102.
23] M.M.Simunin, S.V.Hartov, A.V.Shiverskiy, V.Ya.Zyryanov, Iu.V.Fadeev, A.S.Voronin, Structures on the basis of graphitic nanotubulenes a common electrode in a matrix of porous anodic aluminum oxide for the task of forming switchable electric field membrane, Pisma v Zhurnal Tehnchieskoy Fiziki, 41(2015), 52-59 (in Russian).
Эффективная модель воды для моделирования больших молекулярных систем
Виктор Е. Зализняк
В статье предлагается простая и эффективная, с вычислительной точки зрения, модель воды. Эта одноточечная модель включает в себя потенциал Леннард-Джонса и усреднённое по углам диполь-дипольное взаимодействие. Для параметризации модели используются экспериментальные данные для воды при различных температурах. Основным целевым параметром была выбрана плотность воды. Параметризация проводилась в области температур от 300 до 350 К и давлений до 10.1 МПа. Свойства воды, вычисленные с помощью предлагаемой модели,, сравниваются с экспериментальными данными и данными, полученными с помощью модели воды Т1Р4Р/2005. Сравнение показывает, что предлагаемая модель достаточно хорошо воспроизводит выбранные свойства. В статье также приводится вычислительная схема метода молекулярной динамики для ИРТ ансамбля. Предлагаемая модель не только воспроизводит основные свойства воды, но и требует значительно меньше вычислений по сравнению со стандартными многоточечными моделями воды. Данная модель предназначена для моделирования больших молекулярных систем, возникающих в различных задачах наногидродинамики.
Ключевые слова: модель воды, молекулярная динамика.
