Journal of Computational Physics 326 (2016) 222233

Contents lists available at ScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

Electronic coarse graining enhances the predictive power of

molecular simulation allowing challenges in water physics to
be addressed
Flaviu S. Cipcigan a,b, , Vlad P. Sokhan b , Jason Crain a,b , Glenn J. Martyna c
a
b
c

School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom
National Physical Laboratory, Hampton Road, Teddington, Middlesex TW11 0LW, United Kingdom
IBM T. J. Watson Research Center, Yorktown Heights, NY 10598, USA

a r t i c l e

i n f o

Article history:
Received 20 March 2016
Received in revised form 10 August 2016
Accepted 24 August 2016
Available online 29 August 2016
Keywords:
Molecular dynamics
Force eld
Electronic coarse graining
Dispersion
Polarisation
Path integral

a b s t r a c t
One key factor that limits the predictive power of molecular dynamics simulations is
the accuracy and transferability of the input force eld. Force elds are challenged by
heterogeneous environments, where electronic responses give rise to biologically important
forces such as many-body polarisation and dispersion. The importance of polarisation
in the condensed phase was recognised early on, as described by Cochran in 1959
[Philosophical Magazine 4 (1959) 10821086] [32]. Currently in molecular simulation,
dispersion forces are treated at the two-body level and in the dipole limit, although the
importance of three-body terms in the condensed phase was demonstrated by Barker
in the 1980s [Phys. Rev. Lett. 57 (1986) 230233] [72]. One approach for treating both
polarisation and dispersion on an equal basis is to coarse grain the electrons surrounding a
molecular moiety to a single quantum harmonic oscillator (cf. Hirschfelder, Curtiss and Bird
1954 [The Molecular Theory of Gases and Liquids (1954)] [37]). The approach, when solved
in strong coupling beyond the dipole limit, gives a description of long-range forces that
includes two- and many-body terms to all orders. In the last decade, the tools necessary
to implement the strong coupling limit have been developed, culminating in a transferable
model of water with excellent predictive power across the phase diagram. Transferability
arises since the environment automatically identies the important long range interactions,
rather than the modeller through a limited set of expressions. Here, we discuss the role of
electronic coarse-graining in predictive multiscale materials modelling and describe the
rst implementation of the method in a general purpose molecular dynamics software:
QDO_MD.
2016 Elsevier Inc. All rights reserved.

1. Introduction
Molecular modelling has become an integral part of research in all areas of condensed matter, forming a new part of the
scientic method that elucidates scientic principles and accelerates the discovery of new materials [1]. This computational
revolution has been driven by the exponential scaling of computer hardware [2,3] along with the implementation of scalable
software on these platforms [47]. An important part of the equation remains to improve the accuracy and transferability of

Corresponding author.
E-mail address: aviu.cipcigan@ed.ac.uk (F.S. Cipcigan).

http://dx.doi.org/10.1016/j.jcp.2016.08.030
0021-9991/ 2016 Elsevier Inc. All rights reserved.

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

223

Fig. 1. Examples of environments studied using electronically coarse-grained methods and shown to match experiment: ice II, ambient temperature liquid
water and the surface of liquid water. The images illustrate the electronic responses of individual molecules, with red and blue isosurfaces corresponding to
regions of enhancement and depletion of electronic density, respectively. (For interpretation of the references to colour in this gure, the reader is referred
to the web version of this article.)

the employed force elds, which ultimately dene the accuracy of results, thereby laying the basis for predictive modelling,
given sucient sampling of phase space.
Standard force elds, designed to be numerically ecient, are created by tting a subset of thermodynamic and structural
properties to an analytic expression for the potential energy, typically consisting of sums of pairwise long-range terms [8].
In many systems, this approach has proved successful even if the accuracy of such force elds is at best around 1 kcal/mol.
For example, Schames et al. [9] used molecular dynamics simulations to discover a hidden binding site in HIV integrase,
a molecule that enables the integration of HIVs genetic material into the host cell DNA. This discovery led to experimental
conrmation [10] and subsequent development of raltegravir, a successful medicine that halts the progression of HIV into
AIDS. Another success was the simultaneous discovery by simulation [11] and experiment [12] of the molecular basis of
mutations in a gene responsible for Gauchers disease, a genetic disease causing the build-up of fatty substances in the
organism.
Nonetheless, strategies based on tting to a xed functional are not guaranteed to be transferable to heterogeneous
environments far from the region of parametrisation and validity of terms selected. One such environment is, for example,
a bacterial membrane penetrated by a peptide. The shape of the peptide inside the membrane depends strongly on the
force elds used [13], with only experiments currently being able to distinguish between possibilities.
Here we present our recent development of a class of force elds for materials modelling [1420] that are transferable,
physically motivated, and implemented in the general purpose molecular dynamics software QDO_MD. The technique treats
many-body polarisation and many-body dispersion on the same footing by representing electronic distributions of individual
atoms and molecular moieties using a single coarse grained particle. This particle, known as a Quantum Drude Oscillator
[19,20] (QDO), consists of a negative charge bound to a positive centre by a quantum harmonic oscillator. QDOs generate
many-body polarisation, many-body dispersion interactions and cross terms beyond the dipole limit [21] and represent a
non-perturbative approach that generates all long-range forces within Gaussian statistics when sampled in strong coupling.
The QDO approach has been employed to construct a general purpose water model (QDO-water) that is transferable
without tting to any condensed phase data (at present neglecting bond breaking and making, which could however be
included using a TersoffPettifor approach [17,22,23]). QDO-water consists of a rigid molecular geometry (decorated with
xed charges), short-range empirical repulsion and a QDO. The model is parametrised using the static multiple moments of
an isolated molecule, its polarisability, the dipoledipole dispersion coecient (C 6 ) and one cut through the dimer energy
surface [14] used to set the short-range repulsion. Its properties have been studied in a diverse set of environments (depicted in Fig. 1), ranging from high-pressure ice (ice II), liquid water, steam, the liquidvapour interface, and supercritical
water. The model predicts density, surface tension, enthalpy of vaporisation and dielectric constant across these wide range
of environments, demonstrating a formerly unrealised level of transferability.
2. Predictive materials modelling using electronically coarse grained methods
Fig. 2 presents the spectrum of modelling techniques, with the QDO technique occupying a previously vacant place
between ab initio methods and all-atom molecular dynamics. QDOs guiding philosophy is instead of treating an exact
system in an approximate way, treat an approximate system in an exact way.
In a physical sense, the QDO model is solvable in strong coupling, allowing for an approach that contains all long-range,
many-body interactions to all orders. The absence of truncation means that the modeller does not pick the symmetry by
choice of truncation but rather the environment picks the key terms. This type of simplied but rich model is at the heart
of modern theories because its lack of bias allows the essential physics to emerge naturally without the modeller choosing
the important interactions a priori.

224

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

Fig. 2. Various simulation techniques arranged according to the timescales and lengthscales which they sample. The left-most image, originally by Chaplin
[24], represents electron densities of the water pentamer obtained via ab initio calculations. The next is the electronically coarse grained representation
of the air-water surface, where red and blue isosurfaces depict regions of enhancement and depletion of electronic density. This is followed by a classical
molecular dynamics simulation and an atomically coarse grained study of the same interface. The right-most image represents a particle-based realtime
simulation of water ow over complex landscape using a method developed by Chentanez and Mller [25]. The left and right-most images are used with
permission from their respective authors while the rest are original work by the authors of this paper. (For interpretation of the references to colour in this
gure, the reader is referred to the web version of this article.)

The price paid for the richness of the QDO model at low computational cost (order N log N [26,27] or N [28]) is the
approximation made in coarse graining: assuming localised electrons in an insulator respond according to Gaussian statistics
and the neglect of exchange effects, which are taken to be short range. This is consistent with the general spirit of coarse
grained approaches where exactness is sacriced while maintaining sucient accuracy for the problem at hand. Basically,
QDOs are a coarse graining of a full high level electronic structure without truncating the long range interaction, while
standard force elds can be thought of as a coarse graining of QDOs, with a truncation to dipole limit pair-wise dispersion
and a neglect of polarisation or its inclusion at the many-body dipole level. As demonstrated in the next section, systems
of single QDOs per atom or molecular moiety solved in strong coupling are sucient to model molecules that are isotropic
(such as noble gases) or close to isotropic (such as simple hydrides like methane and water). To model anisotropic responses,
multiple interacting QDOs can be pinned to individual points in the molecular frame [29].
The quantum Drude oscillator model has a rich history dating back to London [30], who used it to derive from quantum
mechanical principles the 1/ R 6 functional form of the dispersion interaction between atoms and molecules and of course
Drude himself, who used the classical limit of the model to understand the dispersion of light passing through insulators
[31]. In more details, the classical limit describes dipole many-body polarisation and has been successfully employed in theory and simulation for many years following the work of Cochran [3236]. In the quantum dipole limit, the Drude oscillator
has been employed as an exemplar model for dispersion in the 1950s [37,38], following Londons initial calculations. The
quantum dipole limit model has been used to simulate xenon clusters [35] via an N 3 scaling method and electron attachment to water clusters using a Conguration Integral approach with high order polynomial scaling [39]. The full model was
examined as an approximate treatment for dispersion of one valance electron systems by Fontana in the 1960s [40]. In a
series of papers beginning in 2006, the tools to simulate the full model with order N to N log N complexity (depending
on the method used to treat long-range interactions) via non-perturbative path integral techniques [20,17,21] were developed. Applications begun with rare gases [19,21] and moved forward to liquid water [18,1416]. In related work, the
dipole-limit quantum model has been embedded within density functional theory (within the local density approximation)
by Tkachenko et al. [41] to treat many-body dipole dispersion beyond empirical potentials using the adiabatic connection
formula [42]. This merger of embedded dipole-limit quantum oscillators and density functional theory lead to applications
to heterogeneous systems such as benzene adsorbed onto metal surfaces [43] or DNAgraphene interactions [44].
Water is an ideal system in which to illustrate the importance of a complete treatment of electronic responses. Waters
structure arises due to a competition between directional hydrogen bonds, which favour an open ice-like local structure and
van der Waals forces, which favour a close-packed local structure. Many-body polarisation and dispersion interaction are
key to distinguishing the dominant motifs under the conditions of interest.
Hydrogen bonds are known to be cooperative, meaning that their interaction strength changes depending on environment
thus leading to the emergence of a wide variety of motifs. A simple reporter of this cooperation is the molecular dipole
moment, which changes from a value of 1.85 D in the gas phase to an estimated value of 2.53.6 D in the condensed phase
[18,45], where four hydrogen bonded motifs are dominant. This sensitivity to local structure can be used as a ngerprint of
various local motifs, as we have reported in Ref. [15]. There, we showed that a cusp in the dipole moment as a function
of temperature coincides with experimental maxima in heat capacity, separating the supercritical region in gas-like and
liquid-like regions [46] and revealing the Widom line [47] at a molecular level. We also have shown that such cooperation
makes hydrogen bonds more asymmetric than previously thought [16], with a molecule favouring the loss of an acceptor
bond (on the side of the oxygen) over that of a donor bond (on the side of the hydrogen) in both the liquid phase and at
the liquidgas interface.
Dispersion interactions are also important in water. Including these responses is essential to generate even the basic
structure of water at room temperature. An illustrative example of this balance is the overstructuring of room-temperature

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

225

water by local approximations to DFT [48], where including the electron correlations that lead to accurate van der Waals
interactions is challenging (although, there have been promising results using the similar technique of embedded quantum
harmonic oscillators at the dipole level [41]). In our water model, it is possible to increase the dispersion coecients linearly
while keeping the polarisabilities constant and vice versa, allowing the study to the effect of both contributions to water
structure independently. In Ref. [18], we showed that by tuning only the strength of van der Waals interactions, the density
of liquid water can be decreased by 15% while still maintaining the desired hydrogen bonded local structure.
In addition, three-body effects account for as much as 25% of the binding energy of water, according to the estimates
of Keutsch et al. [49]. Thus, monoatomic two-body potentials for the water molecule can only reproduce waters condensed
phase properties if their parameters are allowed to vary with state point [50,51]. Even a full-atom description of water
cannot be transferable if its electrostatics is xed. Vega and Abascal [52] showed that a nonpolarisable model of water
cannot simultaneously reproduce the melting temperature and temperature of maximum density.
Therefore, in order to create a truly transferable energy surface for molecular simulations, which is therefore predictive
outside the regime of parametrisation, the electronic effects need to be taken into account since they give rise to non-trivial
many-body polarisation and dispersion. QDOs do so while keeping eciency close to that of all-atom molecular dynamics
and scaling as N log N with the number of atoms N (assuming the use of particle mesh type Ewald summation techniques;
the scaling can be further reduced to order N via the use of cell multipole methods [28]). This strategy can be used to
create models that are parametrised without any condensed phase input, yet generate the properties of condensed phases.
3. Electronic coarse graining using Quantum Drude Oscillators
3.1. Dispersion interactions between quantum harmonic oscillators
In order to gain intuition into how quantum harmonic oscillators reproduce electronic responses, consider a simple
model proposed rst by London [30,37]: a negative charge q of mass m, free to move in a dimension x. The charge is
localised by connecting it by a spring of frequency to a positive charge +q xed at the origin.
Consider adding a positive, xed test charge Q at a distance R  x from the origin. The system is in equilibrium when
the Coulomb force Eq on the negative charge equals the spring force m2 x. Balancing the two forces results in the oscillator
acquiring a dipole moment = qx = q2 E /(m2 ) and thus having a dipole polarisability dened as:

1 :=

q2
m 2

(1)

In order to account for dispersion effects, the system has to be treated quantum mechanically. Consider two identical
oscillators separated by a distance R interacting via a term c ( R ). The leading order interaction is dipoledipole, meaning
that the term c ( R ) is proportional to R 3 . Thus, the Hamiltonian is:

=
H

p 21
2m

p 22
2m



1
+ m2 x21 + x22 + c ( R ) x1 x2 .
2

Changing coordinates to = 1
monic oscillators of frequency


=
H

p 2+
2m

p 1 p 2 and = 1 (x1 x2 ) decouples the Hamiltonian into two independent har-

1 c.
 

(2)

+ m (1 + c ) + +
2

p 2
2m


2

+ m (1 c )
2

(3)

The ground state energy of the coupled system is hence given by, to leading order in c:

E0 =

1
2

1+c+

1c



1
h 1 c 2 + O (c 3 ) .

(4)

Fig. 3 shows the ground state probability distribution 2 (x1 , x2 ) in the independent (c = 0) and correlated (c > 0) cases.
When c increases, it becomes more probable for

x1 and x2 to have the opposite sign and less probable for them to have the
same sign. Since the energy is proportional to dx1 dx2 2 (x1 , x2 )x1 x2 c ( R ), this correlation leads to the attractive force between two neutral molecules proportional to c ( R )2 R 6 , just as the leading order term in the van der Waals potential.
This analysis shows that even a one dimensional quantum harmonic oscillator can capture two basic effects of long range
forces: polarisation and dispersion. In this case, these interactions are limited to dipole-limit forces. To capture the full set
of interactions one has to consider the full three dimensional system.

226

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

Fig. 3. The ground state probability density of two one dimensional quantum harmonic oscillators interacting via a c ( R )x1 x2 term. The left panel illustrates
the non-interacting case, with c ( R ) = 0. The right panel illustrates the interacting case, with c ( R ) > 0. The spread of the wavefunction across the diagonal
represents the electron correlation the gives rise to attractive van der Waals forces in real systems.

3.2. Full quantum model

A Quantum Drude Oscillator (QDO) is made out of a light, negative particle of charge q connected to a heavy positive
centre by a harmonic spring of frequency [19]. The bound state of this system is a drudon, with reduced mass m and
Hamiltonian

=
H

p 2

2m

+ m 2 x2 .

(5)

Perturbing this Hamiltonian via a point charge Q at a distance R gives a rst order correction to the energy of zero and
a second order correction of [21]:

E (2 ) ( R ) =


Q 2 l
l =0

2R 2l+2

(6)

l are given by


l1
(2l 1)!!
h

where the polarisabilities

l =

q2
m 2

2m

(7)

Considering two QDOs separated by a distance R gives a rst order correction to the energy of zero and a second order
correction of [21,40]:

E (2 ) ( R ) =

C 2n R 2n

(8)

n =3

with

C6 =

1 1 h ,
4
C 8 = 51 2 h ,

(9)

...
Therefore, the full model has a complex and rich set of long range responses, which are replicated with only three
parameters. This means the responses are correlated and thus one needs to check whether these correlations are satised
in real atoms and molecules.
3.3. Invariant relationships between response coecients
As the responses of QDOs depend only on three parameters, these parameters can be eliminated, resulting in the following invariants that can be used to verify the accuracy of the QDO approximation [21]:

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

227

Fig. 4. Three types of invariant ratios between polarisation and dispersion coecients, predicted by Quantum Drude Oscillators. Polarisation ratios involve
only polarisabilities and analogously for dispersion ratios. Mixed ratios involve both polarisation and dispersion coecients. Deviation from theory is shown
for three types of atoms and molecules: noble gases, alkali metals and small hydrides.




20
9

49
40

2
= 1,
1 3

C8
C 6 C 10
C 6 1
4C 9

= 1,

(10)

= 1.

Comparing these ratios with the responses of real molecules gives the results shown in Fig. 4 (adapted from Ref. [21]).
The agreement is good for the noble gases, with ratios within 10% of experiment, as their shells are fully lled and thus
the distribution of electrons is nearly spherical. The agreement is also within 10% for most alkali metals, for which the last
electron contributes mostly to polarisation. Small hydrides such as water and methane agree within 15%. This agreement is
due to the electronegativity of atoms such as oxygen, which centres most of the electronic charge on them thus resulting in
charge distributions and responses that are close to spherical. For example, the dipole polarisability tensor of water exhibits
a 4% deviation from spherical symmetry.
The agreement shows that, in the case of the materials illustrated in Fig. 4, QDOs produce a dispersion and induction
series that matches experiments and ab initio calculations. This agreement also means that, in real materials, the individual
terms in these series are not independent, but related by symmetry constraints, with spherically symmetric materials having
responses approximated well by a Gaussian model. However, a spherical model has its limitations. For instance, the dipole
hyperpolarisability is identically zero due to the use of the on-site harmonic interactions. In order to move beyond the
response of (approximately) spherically symmetric systems, more QDOs can be placed the molecular frame to break the
symmetry. Enriching the model by changing the on-site interaction or introducing multiple QDOs per molecular moiety,
forms the basis for future work.
3.4. Parametrisation of the full quantum model
Invariants aside, the responses of a QDO can be inverted to give its properties as a function of the responses. The relevant
relations, which are used to parametrise a QDO are as follows:

1 4C 6
h 312

h 31
h 5C 6
or m =
,
42
C8

q = m 2 1 .

m=

(11)
(12)
(13)

4. Simulating Quantum Drude Oscillators via two temperature path integral molecular dynamics
Path integral molecular dynamics, shortened to PIMD, is a method of sampling a quantum mechanical partition function
such as that of a QDO, using classical methods [53,54]. This section highlights its basic derivation, with the full, QDO-based
method being described in Refs. [20,17].

228

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

Consider a quantum mechanical system (such as electrons sampling the BornOppenheimer surface of a nuclear conguration), with degrees of freedom x. The state of this system is described by a density matrix (x, x ; D ). The partition
function at inverse temperature D = 1/kT D is thus [55]:

Z ( D ) = tr e D H =

dx (x, x; D ),

(14)

where dx represents a small volume element in the dim(x)-dimensional space spanned by the systems degrees of freedom.
Since we are interested in the ground state electronic surface of insulators, the inverse temperature D associated with the
electronic degrees of freedom is taken as large as possible while ensuring the ground state dominates the statistics. For
harmonic systems, D h > 8.
Since the density matrix is an exponential of the Hamiltonian, it can be factorised into higher temperature (smaller
= D / P ) slices to give the following path integral representation,

Z ( D ) =

P 


dxi (xi , xi +1 ; D / P ).

(15)

i =1

The reason for this factorisation is that increasing the temperature, thus decreasing D makes density matrices more
classical and thus easier to approximate. Hence the use of the smallest D (or largest electronic temperature) that still
keeps the system close to the ground state serves to reduce P . The approximation used comes from a Trotter factorisation

is split into a reference H 0 , with known density matrix

[20,56], where the Hamiltonian H
V (x),

0 = e H 0 , and a perturbation

(x, x ; ) e V (x)/2 0 (x, x ; )e V (x )/2 + O ( 3 ).

(16)

In the case of QDOs it is convenient to set H 0 as the Hamiltonian of the isolated oscillator, with the following reference
density matrix:

0 (x, x ; ) =

P ( ) 2
P ( )
exp P ( )(x x )2 +
(x + (x )2 ) .

(17)

The coecients P ( ) and P ( ) are dened in Appendix A.

Note that this form introduces strong nearest-neighbour coupling between the coordinates x and x . To remove this cou i via a staging transformation with unit Jacobian
pling, the density matrix is diagonalised to the independent coordinates u
[20]. The coecients of the transformation are given in Appendix A, leading to the following partition function:

Z ( D ) =

P 


dxi 0 (xi , xi +1 ; ) exp

i =1

P 



i
du

3N /2

i =1


exp

2i2

i =1

 2i
u

V (xi )


exp

2i2


V (xi (
u i ))

(18)

i =1

staging coordinates

external potential

One step remains in order to make the partition function isomorphic to a ctitious classical system: the addition of
 i with corresponding faux masses m.
This transforms the partition function to:
conjugate momenta p

Z ( D ) =

P 

i =1

i
du

1
2i2

3N /2

exp

 2i
u
2i2

exp

P

i =1

 


 2i
p
 i exp
V (xi (
u i )) C
dp
.

2m

(19)

This nal transformation leads to a partition function with the effective classical Hamiltonian H (faux) , which can be sampled
via existing methods to sample classical systems [54,19,20].

H (faux) =

P

 2i
 2i
p
u
V (xi (
u i ))
+
+
.
2

2m
P
2 i
i =1

(20)

of the
When adding the motion of the nuclei, the adiabatic principle needs to be invoked and the ctitious mass m
electronic degrees of freedom selected such that they evolve on a time scale a factor of smaller than that of the nuclear
degrees of freedom. This selection of masses results in two effects: it allows the electronic degrees of freedom to be thermostatted at a higher temperature than the nuclei, and reduces energy transfer between the hot electrons and cold nuclei,
allowing the dynamics of these two subsystems to be evolved simultaneously [54,19,20,17]. Further energy transfer between
the two subsystems can be reduced by using an adaptive thermostat that dissipates excess heat induced by Brownian motion [57]. It is important to recognise that the canonical ensemble must be generated for the electronic degrees of freedom
in order for the model to have physical meaning, since the path integral method is explicitly derived to sample a canonical
ensemble.

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

229

Table 1
The free parameters of QDO-water. E h 27.211 eV is the Hartree energy,
a0 0.5292 is the Bohr radius and e 1.60 1019 C is the electron
charge.
Parameter

Value

Molecular geometry
R OH
0.9572

HOH
104.52
Ground state electrostatics
qH
0.605|e |
R OM
0.2667
Quantum Drude Oscillator
mD
0.3656 me
D
0.6287 E h /h
qD
1.1973|e |

Parameter

D
C
1
1

2
2

Value

Coulomb damping
= H = M
0.1 a 0
1.2 a0
Short range repulsion
613.3 E h
2.3244 a01
10.5693 E h
1.5145 a01

5. QDO_MD: a molecular dynamics software implementing Quantum Drude Oscillators

The QDO methodology has been implemented in QDO_MD, a new software based on the work of the Tuckerman and
Martyna groups that includes a classical molecular dynamics engine, CarParrinello ab initio molecular dynamics (CPAPIMD)
[58] using plane wave basis set, as well as a QM/MM capability with N log( N ) scaling. The force eld based part of the software has the standard force elds built in, targeting complex chemical and biological applications. QDO_MD is based on
a number of algorithms developed by its principal authors, which are now widely accepted, including reversible reference
system propagator (RESPA) [59], NosHoover chains [60] and isothermalisobaric ensemble [61] and allows simulation in
various equilibrium statistical ensembles, perform energy and structure minimisation using several techniques including
steepest descent, conjugate gradient or direct inversion in the iterative subspace. It is parallelised [62] at several levels including classical force parallelisation, parallelisation over beads in path integral simulation, state and hybrid state/reciprocal
space parallelisation of the electrostatic interactions and is aimed at distributed memory architectures using the Message
Passing Interface (MPI) protocol [63]. QDO_MD is undergoing an extensive refactoring and is currently available by contacting the authors. In the near future, it will be placed on the National Physical Laboratory website and open sourced on
GitHub.
QDO_MD has modular structure and easily allows extensions conforming to the standard structure of atomistic MD
software. In order to perform adiabatic path integral molecular dynamics with QDO (APIMD-QDO), a new simulation type
keyword, qdo_pimd, has been added. Two new atom types, drude_bead and drude_center specify the drudon. For
generality, the Drude oscillator centres are treated as ghost atoms. Electrostatic interactions in QDO_MD are calculated using
smooth particle mesh Ewald (PME) method [27], and in addition to standard 3D-periodic systems the code can handle the
systems of reduced dimensionality using modied reciprocal space summation [64,65]. Although drudons do not introduce
new types of interactions, special terms have been added between the bead and the centre in drudon to correct for terms
necessary included in the reciprocal space part of the Ewald scheme. In addition to that, Coulomb regularisation has been
added at short range as well as triexponential repulsion terms to correct for the missing exchange repulsion. The pressure
tensor, which is calculated in QDO_MD from the atomic virial expression is currently being converted to a virial form based
on heavy atom centres. Following the code convention, all interactions are splined.
APIMD-QDO integration is done in the canonical ensemble using a velocity Verlet propagator with RESPA and due to
separation of the nuclear and QDO degrees of freedom this implementation requires separate thermostats for the Drude and
atomic degrees of freedom. The separation is controlled by the adiabaticity parameter , with QDO degrees of freedom kept
at elevated temperature. In order to remove large vacuum energy of the Drude subsystem, low variance harmonic staging
virial estimator [20] acting on QDO sites only has been implemented to sample the instantaneous energy.
APIMD-QDO can handle mixed systems, where only a subsystem contains QDO with the rest treated classically. Since in
the APIMD-QDO algorithm all ingredients for the classical Drude oscillator are included, QDO_MD could also simulate the
classical polarisable molecules. These, however, cannot be mixed with QDOs at present.
6. QDO-water: an electronically-coarse grained model of water
The rst application of QDOs to complex systems was water [18,14] after proof of concept studies on the noble gases
[19,21] (where the importance of many-body dispersion was recognised early on by Barker [72]). Water was chosen due to
its fundamental importance for life. Its properties have been studied for decades, yet the essential physics underlying many
of its live-giving anomalies remains unresolved [66].
To construct a non-dissociative water model out of QDOs, one needs three elements: a rigid molecular frame with
embedded point charges to replicate the lowest order electrostatic moments, a QDO to replicate electronic responses and
a short range, pairwise repulsion potential. These elements are shown in Fig. 5 with their respective parameters given in
Table 1.

230

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

Fig. 5. Schematic of QDO-water.

In more detail, the frame is xed in the experimental molecular geometry, with charges and distances chosen to give
an exact dipole moment for the isolated molecule (1.85 D) and a best t to the quadrupole moment. The parameters of
the QDO are t using the relations in equation (13) to the dipole and quadrupole polarisabilities and the C 6 induceddipoleinduced-dipole dispersion coecient. The repulsive potential is t by calculating the interaction energy between two
molecules using ab initio and the repulsion-free model and tting the difference to a double exponential 1 e 1 r + 2 e 2 r
[18,14]. The Coulomb interaction between charges is regularised by replacing each point charge X by a Gaussian distribution
of width X .
Since QDO-water is parameterised from the properties of a single molecule and the dimer, condensed phase properties are a prediction rather than a tting target. The model was studied in the following environments: high pressure ice (ice II) [14], liquid water [14], steam [14], the liquidvapour interface [16] and supercritical water [15]. The
agreement with experiment was good, with densities within few percent of experiment. The model predicts the temperature of maximum density of 5.5(2) C (compared to the experimental value of 3.98 C [67]) and the critical point of
{ T C = 649(2) K, C = 0.317(5) g/cm3 } (compared to the experimental value of { T C = 647.096 K, C = 0.322 g/cm3 } [67]). It
also accurately tracks the dielectric constant, enthalpy of vaporisation and surface tension between T = 300 K and T = 600 K
[14].
QDO water also provided insights into the local structure at the liquidvapour interface [16] and in supercritical water
[15]. At the liquidvapour interface the orientation of the water molecules is inuenced by an intrinsic asymmetry in the
way water molecules hydrogen bond. When losing a hydrogen bond, molecules prefer to lose it on the side of the oxygen
(an acceptor bond) rather than on the side of the hydrogen (a donor bond) [16]. This leads to a surface with more oxygen
than hydrogen atoms, which is negatively charged.
In supercritical water [15], the dipole moment was discovered to be a reliable ngerprint of the Widom line extending
the liquidgas transition [15] into the supercritical region. The phase diagram locus where the heat capacity is maximum
coincides with a cusp in the dipole moment as a function of density. This cusp separates liquid-like from gas-like scaling,
tracking the Widom line. This conclusion is general to any polar uid, showing the importance of polarisation in characterising highly heterogeneous materials.
At this stage of development, we have studied a classical molecular model using a rigid geometry for the water
molecules. It has been shown that the rigid frame forms a much better approximation to a fully exible model including nuclear quantum effects than a classical exible model [6870]. Also as noted in [15], the critical constants, melting
point and other thermodynamic quantities do not show strong isotope effects. It is feasible to combine QDO with PIMD to
examine a fully exible model with quantum nuclei, a subject of current interest [71].
7. Conclusion
This paper summarises our recent progress in signicantly enhancing the transferability of force elds describing complex systems of insulators by coarse graining the electrons of a molecular moiety to a single quantum harmonic oscillator, an
approach we call the Quantum Drude Oscillator method. This technique generates, naturally and to all orders, many-body
forces such as dispersion, polarisation and mixed interactions beyond the dipole limit. These interactions are important
in modelling heterogeneous environments such as those occurring in biological systems, where the symmetry of key interactions depends strongly on environment. To demonstrate the application of the model, we highlight the successful
construction of a transferable model of a water molecule. The model is parametrised using only the properties of a single
molecule and those of a dimer yet has excellent predictive power across the phase diagram, from high pressure ice to

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

231

supercritical water. To encourage the adoption of the method, we present its rst implementation in the general purpose
molecular dynamics software QDO_MD.
Acknowledgements
This work was supported by the NPL Strategic Research programme. FSC acknowledges the Scottish Doctoral Training
Centre in Condensed Matter Physics (grant number EP/G03673X/1), NPL and EPSRC for funding under an Industrial CASE
studentship. We acknowledge use of Hartree Centre resources in this work, including those from the Intel Parallel Computing Centre.
Appendix A. Coecients
For a quantum Drude oscillator of mass m, charge q and frequency
following values,

P () =
P () =
f =

2h sinh( f )
2m tanh( f /2)

h
P

, the coecients dened in the text have the

(A.1)

The staging transformation uses the following coecients:

 1 = x1 ,
u
 i = xi xi ,
u
xi =

sinh( h )
sinh[i h ]

x1 +

sinh((i 1) h )
sinh(i h )

(A.2)

x P +1 = x1

12 =
i2 =

xi +1 ,

2m tanh( h /2)
h sinh[(i 1) h ] sinh( h )
m sinh[i h ]

References
[1] E. Winsberg, Science in the Age of Computer Simulation, University of Chicago Press, 2010.
[2] D. Frank, R. Dennard, E. Nowak, P. Solomon, Y. Taur, H.-S.P. Wong, Device scaling limits of Si MOSFETs and their application dependencies, Proc. IEEE
89 (3) (2001) 259288, http://dx.doi.org/10.1109/5.915374.
[3] R. Dennard, F. Gaensslen, H.-N. Yu, V. Rideout, E. Bassous, A. Leblanc, Design of ion-implanted MOSFETs with very small physical dimensions, Proc.
IEEE 87 (4) (1999) 668678, http://dx.doi.org/10.1109/jproc.1999.752522.
[4] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R.D. Skeel, L. Kal, K. Schulten, Scalable molecular dynamics with NAMD,
J. Comput. Chem. 26 (16) (2005) 17811802, http://dx.doi.org/10.1002/jcc.20289.
[5] M.E. Tuckerman, D. Yarne, S.O. Samuelson, A.L. Hughes, G.J. Martyna, Exploiting multiple levels of parallelism in molecular dynamics based
calculations via modern techniques and software paradigms on distributed memory computers, Comput. Phys. Commun. 128 (2000) 333376,
http://dx.doi.org/10.1016/S0010-4655(00)00077-1.
[6] E. Bohm, A. Bhatele, L.V. Kale, M.E. Tuckerman, S. Kumar, J.A. Gunnels, G.J. Martyna, Fine-grained parallelization of the CarParrinello ab initio molecular
dynamics method on the IBM Blue Gene/L supercomputer, IBM J. Res. Dev. 52 (1.2) (2008) 159175, http://dx.doi.org/10.1147/rd.521.0159.
[7] R.V. Vadali, Y. Shi, S. Kumar, L.V. Kale, M.E. Tuckerman, G.J. Martyna, Scalable ne-grained parallelization of plane-wave-based ab initio molecular
dynamics for large supercomputers, J. Comput. Chem. 25 (16) (2004) 20062022, http://dx.doi.org/10.1002/jcc.20113.
[8] D.C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2004.
[9] J.R. Schames, R.H. Henchman, J.S. Siegel, C.A. Sotriffer, H. Ni, J.A. McCammon, Discovery of a novel binding trench in HIV integrase, J. Med. Chem. 47 (8)
(2004) 18791881, http://dx.doi.org/10.1021/jm0341913.
[10] D.J. Hazuda, N.J. Anthony, R.P. Gomez, S.M. Jolly, J.S. Wai, L. Zhuang, T.E. Fisher, M. Embrey, J.P. Guare, M.S. Egbertson, J.P. Vacca, J.R. Huff, P.J. Felock,
M.V. Witmer, K.A. Stillmock, R. Danovich, J. Grobler, M.D. Miller, A.S. Espeseth, L. Jin, I.-W. Chen, J.H. Lin, K. Kassahun, J.D. Ellis, B.K. Wong, W. Xu,
P.G. Pearson, W.A. Schleif, R. Cortese, E. Emini, V. Summa, M.K. Holloway, S.D. Young, J.M. Con, A naphthyridine carboxamide provides evidence for
discordant resistance between mechanistically identical inhibitors of hiv-1 integrase, Proc. Natl. Acad. Sci. USA 101 (31) (2004) 1123311238, http://
dx.doi.org/10.1073/pnas.0402357101.
[11] M.N. Offman, M. Krol, I. Silman, J.L. Sussman, A.H. Futerman, Molecular basis of reduced glucosylceramidase activity in the most common gaucher
disease mutant, n370S, J. Biol. Chem. 285 (53) (2010) 4210542114, http://dx.doi.org/10.1074/jbc.m110.172098.
[12] R.R. Wei, H. Hughes, S. Boucher, J.J. Bird, N. Guziewicz, S.M.V. Patten, H. Qiu, C.Q. Pan, T. Edmunds, X-ray and biochemical analysis of N370S mutant
human acid-glucosidase, J. Biol. Chem. 286 (1) (2010) 299308, http://dx.doi.org/10.1074/jbc.m110.150433.
[13] Y. Wang, T. Zhao, D. Wei, E. Strandberg, A.S. Ulrich, J.P. Ulmschneider, How reliable are molecular dynamics simulations of membrane active antimicrobial peptides?, Biochim. Biophys. Acta, Biomembr. 1838 (9) (2014) 22802288, http://dx.doi.org/10.1016/j.bbamem.2014.04.009.

232

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

[14] V.P. Sokhan, A.P. Jones, F.S. Cipcigan, J. Crain, G.J. Martyna, Signature properties of water: their molecular electronic origins, Proc. Natl. Acad. Sci. USA
112 (20) (2015) 63416346, http://dx.doi.org/10.1073/pnas.1418982112.
[15] V.P. Sokhan, A. Jones, F.S. Cipcigan, J. Crain, G.J. Martyna, Molecular-scale remnants of the liquidgas transition in supercritical polar uids, Phys. Rev.
Lett. 115 (11) (2015) 117801, http://dx.doi.org/10.1103/physrevlett.115.117801.
[16] F.S. Cipcigan, V.P. Sokhan, A.P. Jones, J. Crain, G.J. Martyna, Hydrogen bonding and molecular orientation at the liquidvapour interface of water, Phys.
Chem. Chem. Phys. 17 (14) (2015) 86608669, http://dx.doi.org/10.1039/c4cp05506c.
[17] A. Jones, J. Crain, F. Cipcigan, V. Sokhan, M. Modani, G. Martyna, Electronically coarse-grained molecular dynamics using quantum Drude oscillators,
Mol. Phys. 111 (2223) (2013) 34653477, http://dx.doi.org/10.1080/00268976.2013.843032.
[18] A. Jones, F. Cipcigan, V.P. Sokhan, J. Crain, G.J. Martyna, Electronically coarse-grained model for water, Phys. Rev. Lett. 110 (22) (2013) 22781, http://
dx.doi.org/10.1103/physrevlett.110.227801.
[19] T.W. Whiteld, G.J. Martyna, A unied formalism for many-body polarization and dispersion: the quantum Drude model applied to uid xenon, Chem.
Phys. Lett. 424 (46) (2006) 409413, http://dx.doi.org/10.1016/j.cplett.2006.04.035.
[20] T.W. Whiteld, G.J. Martyna, Low variance energy estimators for systems of quantum Drude oscillators: treating harmonic path integrals with large
separations of time scales, J. Chem. Phys. 126 (2007) 074104, http://dx.doi.org/10.1063/1.2424708.
[21] A.P. Jones, J. Crain, V.P. Sokhan, T.W. Whiteld, G.J. Martyna, Quantum drude oscillator model of atoms and molecules: many-body polarization and
dispersion interactions for atomistic simulation, Phys. Rev. B 87 (14) (2013) 144103, http://dx.doi.org/10.1103/physrevb.87.144103.
[22] J. Tersoff, New empirical approach for the structure and energy of covalent systems, Phys. Rev. B 37 (12) (1988) 69917000, http://dx.doi.org/
10.1103/physrevb.37.6991.
[23] D.G. Pettifor, I.I. Oleinik, Analytic bond-order potentials beyond TersoffBrenner. I. Theory, Phys. Rev. B 59 (13) (1999) 84878499, http://dx.doi.org/
10.1103/PhysRevB.59.8487.
[24] M. Chaplin, Water structure and science, www1.lsbu.ac.uk/water/ (Accessed on 15 March 2016).
[25] N. Chentanez, M. Mller, Real-time eulerian water simulation using a restricted tall cell grid, ACM Trans. Graph. 30 (2011) 8218210, http://dx.doi.org/
10.1145/2010324.1964977.
[26] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: an nlog(n) method for Ewald sums in large systems, J. Chem. Phys. 98 (12) (1993) 1008910092,
http://dx.doi.org/10.1063/1.464397.
[27] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys. 103 (19) (1995)
85778593, http://dx.doi.org/10.1063/1.470117.
[28] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (2) (1987) 325348, http://dx.doi.org/10.1016/
0021-9991(87)90140-9.
[29] A.J. Stone, A.J. Misquitta, Atomatom potentials from ab initio calculations, Int. Rev. Phys. Chem. 26 (1) (2007) 193222, http://dx.doi.org/10.1080/
01442350601081931.
[30] F. London, The general theory of molecular forces, Trans. Faraday Soc. 33 (8) (1937) 8b, http://dx.doi.org/10.1039/tf937330008b.
[31] P. Drude, The Theory of Optics, Longmans, Green & Co., 1902.
[32] W. Cochran, Lattice dynamics of alkali halides, Philos. Mag. 4 (45) (1959) 10821086, http://dx.doi.org/10.1080/14786435908238288.
[33] M. Sangster, M. Dixon, Interionic potentials in alkali halides and their use in simulations of the molten salts, Adv. Phys. 25 (3) (1976) 247342, http://
dx.doi.org/10.1080/00018737600101392.
[34] M. Sprik, M.L. Klein, A polarizable model for water using distributed charge sites, J. Chem. Phys. 89 (12) (1988) 7556, http://dx.doi.org/10.1063/
1.455722.
[35] G.J. Martyna, B.J. Berne, Structure and energetics of Xen : many-body polarization effects, J. Chem. Phys. 88 (7) (1988) 4516, http://dx.doi.org/10.1063/
1.453759.
[36] Z.H. Deng, G.J. Martyna, M.L. Klein, Quantum simulation studies of metal-ammonia solutions, J. Chem. Phys. 100 (10) (1994) 75907601, http://
dx.doi.org/10.1063/1.466852.
[37] J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, The Molecular Theory of Gases and Liquids, Wiley-Interscience, 1954.
[38] M. Sparnaay, On the additivity of Londonvan der Waals forces, Physica 25 (16) (1959) 217231, http://dx.doi.org/10.1016/S0031-8914(59)92714-4.
[39] F. Wang, K.D. Jordan, Application of a Drude model to the binding of excess electrons to water clusters, J. Chem. Phys. 116 (16) (2002) 69736981,
http://dx.doi.org/10.1063/1.1461811.
[40] P.R. Fontana, Theory of long-range interatomic forces. I. Dispersion energies between unexcited atoms, Phys. Rev. 123 (5) (1961) 18651870, http://
dx.doi.org/10.1103/physrev.123.1865.
[41] A. Tkatchenko, R.A. DiStasio, R. Car, M. Scheer, Accurate and ecient method for many-body van der Waals interactions, Phys. Rev. Lett. 108 (23)
(2012) 236402, http://dx.doi.org/10.1103/physrevlett.108.236402.
[42] A. Tkatchenko, A. Ambrosetti, R.A. DiStasio, Interatomic methods for the dispersion energy derived from the adiabatic connection uctuation
dissipation theorem, J. Chem. Phys. 138 (7) (2013) 074106, http://dx.doi.org/10.1063/1.4789814.
[43] W. Liu, J. Carrasco, B. Santra, A. Michaelides, M. Scheer, A. Tkatchenko, Benzene adsorbed on metals: concerted effect of covalency and van der Waals
bonding, Phys. Rev. B 86 (24) (2012) 245405, http://dx.doi.org/10.1103/physrevb.86.245405.
[44] Y. Cho, S.K. Min, J. Yun, W.Y. Kim, A. Tkatchenko, K.S. Kim, Noncovalent interactions of DNA bases with naphthalene and graphene, J. Chem. Theory
Comput. 9 (4) (2013) 20902096, http://dx.doi.org/10.1021/ct301097u.
[45] D.D. Kemp, M.S. Gordon, An interpretation of the enhancement of the water dipole moment due to the presence of other water molecules, J. Phys.
Chem. A 112 (22) (2008) 48854894, http://dx.doi.org/10.1021/jp801921f.
[46] G.G. Simeoni, T. Bryk, F.A. Gorelli, M. Krisch, G. Ruocco, M. Santoro, T. Scopigno, The Widom line as the crossover between liquid-like and gas-like
behaviour in supercritical uids, Nat. Phys. 6 (7) (2010) 503507, http://dx.doi.org/10.1038/nphys1683.
[47] B. Widom, Surface tension of uids, in: C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, Academic Press, 1972, pp. 7999.
[48] I.-C. Lin, A.P. Seitsonen, I. Tavernelli, U. Rothlisberger, Structure and dynamics of liquid water from ab initio molecular dynamicscomparison of BLYP,
PBE, and revPBE density functionals with and without van der Waals corrections, J. Chem. Theory Comput. 8 (10) (2012) 39023910, http://dx.doi.org/
10.1021/ct3001848.
[49] F.N. Keutsch, J.D. Cruzan, R.J. Saykally, The water trimer, Chem. Rev. 103 (7) (2003) 25332578, http://dx.doi.org/10.1021/cr980125a.
[50] M.E. Johnson, T. Head-Gordon, A.A. Louis, Representability problems for coarse-grained water potentials, J. Chem. Phys. 126 (14) (2007) 144509, http://
dx.doi.org/10.1063/1.2715953.
[51] A. Chaimovich, M.S. Shell, Anomalous waterlike behavior in spherically-symmetric water models optimized with the relative entropy, Phys. Chem.
Chem. Phys. 11 (12) (2009) 1901, http://dx.doi.org/10.1039/b818512c.
[52] C. Vega, J.L.F. Abascal, Relation between the melting temperature and the temperature of maximum density for the most common models of water, J.
Chem. Phys. 123 (14) (2005) 144504, http://dx.doi.org/10.1063/1.2056539.
[53] M. Parrinello, A. Rahman, Study of an f center in molten KCl, J. Chem. Phys. 80 (2) (1984) 860, http://dx.doi.org/10.1063/1.446740.
[54] M.E. Tuckerman, B.J. Berne, G.J. Martyna, M.L. Klein, Ecient molecular dynamics and hybrid Monte Carlo algorithms for path integrals, J. Chem. Phys.
99 (4) (1993) 2796, http://dx.doi.org/10.1063/1.465188.

F.S. Cipcigan et al. / Journal of Computational Physics 326 (2016) 222233

233

[55] R.P. Feynman, Statistical Mechanics: A Set of Lectures, (Advanced Books Classics), Westview Press, 1998.
[56] H.F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10 (4) (1959) 545, http://dx.doi.org/10.1090/s0002-9939-1959-0108732-6.
[57] A. Jones, B. Leimkuhler, Adaptive stochastic methods for sampling driven molecular systems, J. Chem. Phys. 135 (8) (2011) 084125, http://dx.doi.org/
10.1063/1.3626941.
[58] L.V. Kal, K. Schulten, R.D. Skeel, G. Martyna, M. Tuckerman, J.C. Phillips, S. Kumar, G. Zheng, Biomolecular modeling using parallel supercomputers, in:
Chapman & Hall/CRC Computer & Information Science Series, Chapman and Hall/CRC, 2005, pp. 34-134-43.
[59] M. Tuckerman, B.J. Berne, G.J. Martyna, Reversible multiple time scale molecular dynamics, J. Chem. Phys. 97 (3) (1992) 19902001, http://dx.doi.org/
10.1063/1.463137.
[60] G.J. Martyna, M.L. Klein, M. Tuckerman, NosHoover chains: the canonical ensemble via continuous dynamics, J. Chem. Phys. 97 (4) (1992) 26352643,
http://dx.doi.org/10.1063/1.463940.
[61] G.J. Martyna, D.J. Tobias, M.L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys. 101 (5) (1994) 41774189, http://dx.doi.org/
10.1063/1.467468.
[62] A.B.K. Collis, A performance analysis of two ab initio molecular dynamics simulation codes, MSc report, The University of Edinburgh, 2009.
[63] M. Forum, MPI a Message Passing Interface Standard V3.1. Ean 1114444410030, High Performance Computing Center, Stuttgart (HLRS), 2015 (Accessed
on 15 March 2016).
[64] G.J. Martyna, M.E. Tuckerman, A reciprocal space based method for treating long range interactions in ab initio and force-eld-based calculations in
clusters, J. Chem. Phys. 110 (6) (1999) 28102821, http://dx.doi.org/10.1063/1.477923.
[65] P. Minry, M.E. Tuckerman, K.A. Pihakari, G.J. Martyna, A new reciprocal space based treatment of long range interactions on surfaces, J. Chem. Phys.
116 (13) (2002) 53515362, http://dx.doi.org/10.1063/1.1453397.
[66] C.A. Angell, Supercooled water: two phases?, Nat. Mater. 13 (7) (2014) 673675, http://dx.doi.org/10.1038/nmat4022.
[67] W. Wagner, A. Pru, The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientic use, J. Phys.
Chem. Ref. Data 31 (2) (1999) 387, http://dx.doi.org/10.1063/1.1461829.
[68] J. Lobaugh, G.A. Voth, A quantum model for water: equilibrium and dynamical properties, J. Chem. Phys. 106 (6) (1997) 24002410, http://dx.doi.org/
10.1063/1.473151.
[69] Y.A. Mantz, B. Chen, G.J. Martyna, Temperature-dependent water structure: ab initio and empirical model predictions, Chem. Phys. Lett. 405 (46)
(2005) 294299, http://dx.doi.org/10.1016/j.cplett.2005.02.050.
[70] S. Habershon, T.E. Markland, D.E. Manolopoulos, Competing quantum effects in the dynamics of a exible water model, J. Chem. Phys. 131 (2) (2009)
024501, http://dx.doi.org/10.1063/1.3167790.
[71] F. Paesani, S. Iuchi, G.A. Voth, Quantum effects in liquid water from an ab initio-based polarizable force eld, J. Chem. Phys. 127 (7) (2007) 074506,
http://dx.doi.org/10.1063/1.2759484.
[72] J.A. Barker, Many-body interactions in rare gases: krypton and xenon, Phys. Rev. Lett. 57 (1986) 230233.