Reverse Monte Carlo modeling of liquid water with the explicit use of the SPC/E interatomic potential

0.83 MB·

by  Ildikó Pethes

Reverse Monte Carlo modeling of liquid water with the explicit use of the
SPC/E interatomic potential
Ildikó Pethes1,a) and László Pusztai1
Institute for Solid State Physics and Optics, Wigner Research Centre for Physics,
Hungarian Academy of Sciences, Konkoly-Thege M. út 29-33, 1121, Budapest,
Hungary
(Dated: 19 January 2017)
Reverse Monte Carlo modeling of liquid water, based on one neutron and one X-ray diffraction data set,
applying also the most popular interatomic potential for water, SPC/E, has been performed. The strictly
7
rigid geometry of SPC/E water molecules had to be loosened somewhat, in order to be able to produce
1
0 a good fit to both sets of experimental data. In the final particle configurations, regularly shaped water
2 molecules and straight hydrogen bonding angles were found to be consistent with diffraction results. It has
been demonstrated that explicit use of interatomic potentials in RMC has a role to play in future structural
n
a modeling of water and aqueous solutions.
J

8
1 I. INTRODUCTION the story started with Monte Carlo simulations nearly
fifty yearsago20 andcontinued with molecular dynamics

] Water, the most common liquid on Earth, has been (MD) (for a review, see e.g. Ref. 21), ab initio MD (see
t
f (one of) the most frequently investigated substance for e.g. Ref. 22) and Reverse Monte Carlo (RMC) (see e.g.
o
thousands of years. The microscopic structure of liquid Ref. 5) simulations. Monte Carloand molecular dynam-
s
. water, which is a fundamental piece of information for ics simulations are based on the intra- and intermolecu-
t
a researchers in many fields of research, is one of the old- lar interactions between particles, thus their outcome is
m est not fully resolved problems. Despite the dozens of determined by the chosen force fields. During the past
- publications appearing year after year, our knowledge is decades, several water models with different force field
d still uncertain in this respect1. The unique behavior of parametershavebeen developedthathavebeenfittedto
n water stems from the hydrogen bonded network of the some chosen experimental data (for reviews, see7,21,23).
o
c molecules. However,crucialfeatures of this network,the In Reverse Monte Carlo structural modeling24 large
[ averagenumber of the hydrogenbonds per molecules,or 3-dimensional particle configurations are generated that

the intermolecularO-Hbonddistancearestillsomewhat areconsistentwithallthesuppliedinputdatasetswithin
1
controversial2–5. their uncertainty. Any experimental (and/or theoret-
v
5 The microscopic structure of liquid water has been ical quantities) that can be expressed in terms of the
3 studied by different spectroscopic and scattering tech- atomic coordinates may be fitted simultaneously. Con-
0 niques, for example small and wide angle X-ray ventionalRMCalgorithms24–27 arenotabletotakeener-
5 scattering6–10, neutron diffraction11–13, and X-ray ab- getic considerations into account, although the so-called
0 sorption and emission spectroscopies2,14–17. X-ray ’hybrid RMC’ scheme by Opletal et al.28,29 does operate
.
1 diffraction is suitable for the determination of oxygen- with specific potential parameters. Another method of
0 related correlations, but it is less sensitive to hydro- structural modeling, the ’Empirical Potential Structure
7 gen. Neutron diffraction with H/D isotopic substitu- Refinement’ (EPSR)30,31, starts with known interatomic
1
tioncanbeusefulforthedetectionofhydrogen-hydrogen potential parameters that are varied during the calcula-
:
v and hydrogen-oxygen correlations, since H has nega- tions. EPSR has been applied to liquid water over the
i tive coherent scattering length, bH = -3.7406 fm, while past 20 years, starting with the original publication in
X c
bD=6.671 fm. However, the determination of the co- 199630, to a very recent extensive paper11.
r c
a herent structure factor from the measured neutron scat- The RMC_POT algorithm (and the software that
tering intensities is difficult due to the large incoher- makes use of it)32 combines traditional RMC modeling
ent scattering cross section of H, and the strong inelas- with some features of standard molecular simulations.
ticity effects caused by the similar masses of H nuclei Instead of ’Fixed Neighbour Constraints’ (FNC) used
and incident neutrons. There are numerous attempts previously (see, e.g., Ref. 26), the RMC_POT algo-
known for resolving these issues, e.g. via oxygen isotope rithm keeps molecules together via (more or less) flex-
substitution18 andpolarizedneutrondiffractionwithpo- ible intramolecular forces: bond stretching, angle bend-
larization analysis19. ing and dihedralstretching potential functions (see, e.g.,
Various computer simulation methods have also been Ref. 33). Besides these, RMC_POT can handle inter-
applied in order to gain real-space correlation functions: molecular potentials of arbitrary complexity: Coulomb
and Lennard-Jones energies can be easily calculated. It
isworthstressingthatwhileinEPSR30 potentialparam-
etersarebeingmodifiedcontinuously,RMC_POTkeeps
a)[email protected] all intermolecular terms intact.
2
Total scattering structure factors (TSSF) obtained B. Molecular dynamics simulations
from neutron and X-ray diffraction experiments were
studiedbytheconventionalRMC technique,usingFNC, Molecular dynamics simulations have been performed
in an early publication5. The RMC technique has also by the GROMACS software package (version 5.1.1)33.
been used to investigate the compatibility of structure The initial particle configuration was obtained by plac-
factorsandpartialpaircorrelationfunctions(PPCF)ob- ing 3333 water molecules in the cubic simulation box by
tainedbydifferentmethodsinRefs. 34and35. Thecon- the ’gmx solvate’ program of the GROMACS package.
ventional RMC technique was used to compare different The edge length of the simulation box was 46.5701 Å,
water potential models via PPCFs obtained by MD sim- according to an atomic number density of 0.099 Å−3.
ulations and structure factors from neutron and X-ray According to the SPC/E water potential39, pairwise-
diffraction measurements in Refs. 36 and 37.
additive non-bonded interactions have been used for the
representation of dispersion and repulsion effects (in the
In this study the suitability of the RMC_POT algo-
form of the Lennard-Jones (LJ) potential), as well as
rithm for the determination of the structure of water is
for the electrostatic interactions. In the SPC/E force
explored. X-ray and neutron diffraction structure fac-
field39, which was chosen for this study, the charges are
tors from Refs. 6, 8, and 38 were applied here as in-
-0.8476eand 0.4238e for O and H atoms, respectively (e
put data. Of the many possibilities, the most frequently
is the elementary charge). The LJ σ and ǫ parameters
usedSPC/E(’ExtendedSimplePointCharge’)waterpo-
tential model39 was selected for the purpose. We are are3.16557Åand0.650194kJ/molfortheoxygenatoms
and0forhydrogens. Intramoleculardistancesare1Åfor
awarethattheintramolecularO-HdistanceoftheSPC/E
O-H and 1.633 Å for H-H pairs in the rigid molecules.
model is very slightly (by about 0.02 Å) longer than
Canonical NVT (constant number, volume and tem-
suggested by most neutron diffraction experiments (see,
perature) ensemble was applied at T = 295 K. The tem-
e.g., Ref. 12); however, the overwhelming popularity of
perature was controlled by the Berendsen thermostat40,
SPC/E over the past decades justifies the choice of this
with temperature coupling time τ = 0.01 ps. The cut-
potential for the first RMC_POT study on liquid wa- T
offsfortheCoulombandvanderWaalsinteractionswere
ter. For generating starting particle arrangements for
10 Å. The steepest-descent gradient method was used
RMC_POT,aswellasfor referencepurposes,molecular
to reach the energy equilibrium. The total energy has
dynamics computer simulations with the SPC/E water
potential39 have been carried out as an initial step. reached its minimum value in less than 100 ps. The to-
tal simulation time was 4000 ps, the time step was 2 fs;
particleconfigurationsforcalculatingaverageshavebeen
collected in every 20 ps between 2000-4000 ps. MD re-
sults reported here were averaged over 100 time frames.
The ’gmx rdf’ programwas used to calculate the partial
II. SIMULATION DETAILS pair correlation functions from the collected MD config-
urations.
MD simulations were conducted with the flexible
A. Experimental data sets SPC/Fw water model41, as well. In this model the LJ
σ and ǫ parameters are the same as in the SPC/E force
field. The charges are -0.82e and 0.41e for O and H
Totalscattering structure factors that hadbeen inves-
atoms, respectively. The flexibility is realized by har-
tigated by conventional RMC earlier5 were tested here:
monicbondstretchingandanglebendingpotentials. The
theX-raystructurefactorofH2OfromNartenandLevy8 equilibrium O-H bond length is 1.012 Å, the k force
b
(this dataset will be denoted throughout this paper as constant is 443153 kJ mol−1 nm−2, the equilibrium H-
XRD1) and the neutron structure factor of D2O from O-H angle is 113.24 °, the k force constant is 317.56 kJ
a
Soper et al.38 (indicated as ND). Additionally, a recent mol−1 rad−2. Details of the simulation were similar as
X-ray diffraction result, the SOO partial structure factor before, but the time step was smaller (0.2 fs). Results
of Skinner and co-workers6 (marked as XRD2) was con-
from this simulation will be referred to as ’MD-flexible’,
sidered. TheNDstructurefunctionwasmodeledoverthe
whilethosefromtheprecedingcalculations(withSPC/E
1.1Å−1≤Q≤15Å−1regime. Theweightsofthepartial
water model) will be denoted as ’MD-rigid’.
functionsintheNDtotalscatteringfunctionwere0.0919
(O-O),0.4225(O-D)and0.4857(D-D). XRD1S(Q)was
considered for Q values 1 Å−1 ≤ Q ≤ 16 Å−1, whereas
C. Reverse Monte Carlo modeling
theSOO partialstructurefactor(XRD2)wasfittedinthe
0.975 Å−1 ≤Q≤ 25.85 Å−1 region.
RMC modeling is described in detail in Refs. 26, 27,
In the first case (denoted as Case 1) the ND and the 32,and42. During a conventionalRMC calculationpar-
XRD1 datasets werefitted together,while inthe second ticles are moved randomly in the simulation box and
set of RMC simulations (marked as Case 2) the ND and differences between experimental and model structural
XRD2structurefactorswereapproachedsimultaneously. quantities are minimized. RMC may be used for any
3
quantity that can be expressed from the atomic coor- j is in the middle) is
dinates (e.g. structure factors from diffraction experi-
ments, EXAFS signals, or model pair correlation func- V (θ )= 1ka θ −θ0 2 (2)
tions). If the squared differences between experimental a ijk 2 ijk ijk ijk
and calculated data sets decrease by the move of a par- (cid:0) (cid:1)
whereka istheforceconstantandθ0 istheequilibrium
ticle then the move is accepted, otherwise it is only ac- ij ijk
angle.
cepted with some probability.
The non-bonded potential energy terms in these sim-
Inthepresentinvestigation,severaldifferentRMCcal-
ulations are the Coulomb and the Lennard-Jonescontri-
culations have been performed; they are summarized in
butions. The Coulomb potential is
Table I.
AlongwithRMCcalculationsthatuseinteratomicpo- 1 qiqj
V (r )= (3)
C ij
tential functions, traditional RMC modeling with Fixed 4πǫ0 rij
Neighbour Constraints (FNC) has also been carried out
for comparison. For these computations the RMC++ whereqi andqj arethe partialchargesplacedonatomsi
computer programme27 has been applied (for an early andj,ǫ0 isthevacuumpermittivity. TheLennard-Jones
application of such an approach for liquid water, see, potential is
e.g., Ref. 5). In RMC++ molecules are kept together
12 6
via FNC; no interatomic potentials are involved. The V (r )=4ǫ σij − σij (4)
LJ ij ij
FNC method connects two hydrogenatoms and the cen- r r
(cid:18) ij(cid:19) (cid:18) ij(cid:19) !
traloxygenatompermanentlyviaasimpleneighborlist.
Intramoleculardistancesarekeptbetweenminimumand where ǫ and σ are the Lennard-Jones parameters ap-
ij ij
maximum values, namely 0.95 to 1.03 Å for (covalently plied for the ij atom pair.
bonded) O-H and 1.55 to 1.70 Å for H-H (non-bonded During RMC_POT calculations differences between
intramolecular)pairs. Closestintermolecularapproaches experimental and model S(Q) functions are minimized
were 2.2 Å for O-O, 1.5 Å for O-H and 1.7 Å for H-H together with the total potential energy. A relative
pairs. The bin size was 0.05 Å, and the maximum move weight (σ parameter) is assigned to every potential re-
of atoms was set to 0.05 Å. The final MD-rigid configu- lated term. From the potential-related terms χ2i =Vi/σi
ration was taken as starting configuration for these sim- values are calculated and the sum of them is minimized
ulations. Control parameters for experimental data sets (i refersto the individualpotentialrelatedcomponents).
are presented in Table II. These runs will be referred to After moving a particle, first the potential-related terms
as ’RMC-FNC’ throughout this work. are investigated and the move is accepted if the sum of
In an additional reference series of RMC calculations the χ2i terms has decreased. If it has increased then the
the FNC method has been used again and in addition, move is accepted with some probability. Calculations
partial pair correlation functions (g (r)) from MD sim- related to conventional data sets (and geometrical con-
ij
ulations have been appliedas ’(quasi-)experimentaldata straints, if present) would be executed only if the move
sets to fit’; these RMC runs will be denoted as ’RMC- can be accepted based on potential energy.
FNC+g(r)’. A similar approach had been used and in- Properchoiceoftheσ parameterscanwarrantthatin-
vestigated previously in conjunction with several water formation related to experimental data and interatomic
modelsinPusztaietal.36andSteinczingeretal.37. g (r) potentialsaretakenintoaccountinabalancedway. Rela-
ij
curveswereobtainedherefrommoleculardynamics sim- tiveweightsofthebondingpotentialtermsshouldbecho-
ulationwith the SPC/Eforcefield(MD-rigid). The con- sen so that molecules are kept together but they can be
trolparametersforthesedatasetsarealsoshowninTable asflexibleasdiffractiondatamightrequire. Weightsthat
II; all other parameters were the same as before. are too strict do not allow the system to move around;
on the other hand, too loose relative weights result in
In the RMC_POT method, exploitation of which is
that the molecules break up. The ratio of the potential
the genuinely novel element of this study, molecules
andexperimentrelatedweights shouldbe chosenso that
are kept together via intramolecular (’bonded’) poten-
the fit to the experimental data sets is as good as possi-
tialsthatarecalculatedsimilarlytothatimplementedin
ble, i.e., at least of similar quality as it may be achieved
GROMACS33. In general, the bond stretching interac-
in an RMC++ calculation. As a result of RMC_POT,
tion (between atoms i and j) is taken into account as a
final particle configurations will be compatible with the
harmonic potential:
experimentaldatasetsandtheappliedforcefieldmodels
V (r )= 1kb (r −b )2 (1) simultaneously – provided that these two things can be
b ij 2 ij ij ij made compatible at all for a given system.
It should be noted here that as a consequence of the
where kb is the force constant, b is the equilibrium way molecules are handled in the RMC_POT method,
ij ij
distance of the bonded pair, r is the actual distance and of that movements are overwhelmingly atomic in
ij
of the atoms in the bonded pair. The harmonic angle RMC modeling in general, molecules in practice are al-
bendingpotential(betweenatomsi,j,andk,whereatom ways flexible (or perhaps better to say, never strictly
4
TABLE I.Reverse Monte Carlo calculations performed.
RMC calculation Short name Data sets Starting configuration Potential/g(r) data sets
RMC-FNC Case 1 FNC_X1 ND+XRD1 MD-rigid -
RMC-FNC Case 2 FNC_X2 ND+XRD2 MD-rigid -
RMC-FNC+g(r) Case 1 FNC_g_X1 ND+XRD1 MD-rigid g(r) sets
RMC-FNC+g(r) Case 2 FNC_g_X2 ND+XRD2 MD-rigid g(r) sets
RMC_POT Case 1 POT_X1 ND+XRD1 MD-rigid SPC/Ef
RMC_POT Case 2 POT_X2 ND+XRD2 MD-rigid SPC/Ef
RMC_POT random POT_r ND+XRD2 random SPC/Ef
RMC_POT SPC/Fw POT_Fw ND+XRD2 MD-flexible SPC/Fw
notbedefinedseparately,sinceintra-andintermolecular
TABLEII.Final controlparameters (σ values) of inputdata
regionsmayoverlap;consideringthispossibility,theH-H
sets and potential terms.
minimal interatomic distance was set at 1.4 Å. Bin size,
Data Set/Potential σ maximum moves of the particles and relative weights of
ND 0.003 the experimental data sets were the same as in the tra-
XRD1 0.003 ditional RMC runs.
XRD2 0.005 As a cross-check with a genuine flexible force field,
LJ 1.0a RMC_POT calculations have been performed applying
Coulomb 1.0a theflexibleSPC/Fw41,alongwiththe’flexible’versionof
Bond 0.5 (the originally rigid) SPC/E, a.k.a. SPC/Ef (see above)
Angle 0.5 potential. The final MD-flexible configurationwas taken
gO−O 0.04 as starting configuration for the RMC_POT runs per-
gO−H,intra 0.1 formed with the SPC/Fw potential (this calculation is
gO−H,inter 0.04 denotedas’RMC_POTSPC/Fw’). RMC_POTsimula-
gH−H,intra 0.01 tions with SPC/Ef were started from the final MD-rigid
gH−H,inter 0.04 configuration. The (possible) influence of the starting
configurationwas checkedby a special simulation run: a
a Thisvaluewas0.6forPOT_Fw randomconfigurationofwatermoleculeshadbeengener-
atedby placingwatermoleculesrandomlyto the simula-
tion box, and this random initial configuration was used
rigid). Flexible water models are treated the same in this ’RMC_POT random’ run. (We note in passing
way as in an MD simulation. However, instead of the thataccordingtothebestofourknowledge,nosuchtest
rigid SPC/E water model a modified variant of this hasbeenconductedforanyotherpotential-relatedRMC-
model (referred to as SPC/Ef throughout this work) like methodbefore.) This ’RMC_POT random’calcula-
was effective during the present RMC_POT simula- tionrequiredmorecomputationaltime(highernumberof
tions. It was realized in a similar way as in the work acceptedmoves)toreachanequilibrium,thereforeithas
of Teleman and co-workers43: the basic idea is that a onlybeenperformedforonecombinationofpotentialand
rigid model can be treated as a flexible model with in- experimental data sets. During the ’RMC_POT ran-
finitely strong force constants (ka = kb = ∞). In the dom’ simulation the SPC/Ef potential was applied and
SPC/Ef model every single parameter of the SPC/E the ND and XRD2 experimental data sets were fitted.
model are maintained, but finite force constants are in-
Good values of the σ parameters have resulted from
troduced. The force constant for bond stretching (k )
b severaltestruns,inwhichthenumberofacceptedmoves
was 463700 kJmol−1nm−2 and the angle bending force was about 106. In the final runs the number of the ac-
−1 −2
constant (ka) was 383 kJmol rad , in alignment with cepted moveswas around 2−8×107. The final σ values
Ref. 43. The exact values of the force constants are
are shown in Table II.
not critical, since the relative weights of the potential-
related terms in the RMC_POT scheme determine the
significanceofthese’intramolecularbonding’energycon-
tributions to the total χ2 and thus, to the final particle III. RESULTS AND DISCUSSION
configuration.
In the present study version 1.4 of RMC_POT32 has RMC and MD simulated totalstructure factors,along
been made use of. Intermolecular O-O and O-H mini- with their experimental counterparts, are shown in Fig-
mum distances werethe same asfor the RMC-FNC sim- ure 1. Agreement between any RMC calculation and
ulation runs (2.2 Å and 1.5 Å, respectively). For H-H experimental data is very good, whereas differences be-
pairs, inter- and intramolecular closest approaches can- tween MD simulations and experiment are well visible
5
TABLEIII.Goodness-of-fitvalues(R-factors)foreachcalcu-
lation.
(a)
R-factors (in %)
0.5
Q) ND XRD1 XRD2
D ( 0.0 MD-rigid 15.4 19.7 13.5
N
F MD-flexible 13.6 24.4 18.5
FNC_X1 3.25 2.81 -
-0.5
FNC_X2 3.27 - 3.15

FNC_g_X1 4.20 4.14 -
0.4
(b) FNC_g_X2 4.14 - 4.20
POT_X1 3.52 3.54 -
0.0 POT_X2 3.34 - 4.34
Q)
XRD1 ( -0.4 PPOOTT__rFw 33..7205 -- 44..4272
F
-0.8
allRMC_POT calculations(regardlessofthe actualpo-
0.4
(c) tential and starting configuration) appear to be identi-
cal, although the distribution of intramolecular H-O-H
0.0 angles is significantly wider when the flexible SPC/Fw
Q)
XRD2 ( -0.4 e x p e r i m e n t a l PPOOTT__XX21 faoprpcreofiaeclhdpisroavpidpeliesdig.nCifiaclacnutlalytiobnrsoaudseinrgdtishteribRuMtiCon-FsNfoCr
F MD-rigid POT_r these intramolecular parameters. As each RMC calcula-
MD-flexible POT_Fw
tion produced agreement with experimental data at the
FNC_X1 FNC_g_X1
-0.8 FNC_X2 FNC_g_X2 same(veryhigh)level,itisestablishedthatthedataused
here allow for a diversity of the H···H non-bonding in-
0 5 10 15
tramoleculardistance andof the bond angle as shownin
Q [¯]
Figure 2.
Partial pair correlationfunctions obtained for most of
FIG. 1. Experimental, MD simulated and RMC total struc- the (MD and RMC) simulations are presented in Figure
ture factors. (a) Neutron diffraction data of Soper38; (b) X- 3. The two XRD sets bring about very different O-O
ray diffraction data of Narten8; (c) the recent O-O partial PPCF-s for the RMC-FNC calculations; also, none of
structurefactor derived from X-ray data by Skinner6.
the RMC-FNC simulations produced a separated O-H
peak around the (assumed) hydrogen bonding distance
(between cca. 1.75 and 1.95 Å). Otherwise, the curves
reflect, again, the diversity that one had to get used
(although it is worth noting that both the rigid SPC/E to in the literature of the structure of water (see, e.g.,
and the flexible SPC/Fw potentials work remarkably Refs. 5, 11, 35, and 38). O-O and O-H peak positions
well). Goodness-of-fit values (R-factors) are provided in are gathered in table IV: for the former, values between
TableIII:inaccordancewithvisualinspection,R-factors 2.75 and 2.9 Å, while for the latter, those between 1.77
alsoindicatesimilarfitqualitiesforalltheRMC models, and 1.82 Å have been found (if we disregard the most
whereas differences between MD and experiment appear certainly unphysical O-H maxima at 2.6 Å detected for
as ’magnified’. the RMC-FNC runs). Again, the variously shaped func-
Since in an earlier communication35 some dispute ap- tions in Figure 3 and 4 and the different maximum po-
peared concerning intramolecular parameters that may sition values in Table IV that are related to RMC cal-
be derived from diffraction experiments, here we show culations are all equally consistent with the two sets of
distributionsofO-H(bonding)andH···H(non-bonding) diffraction data: separation between them is only possi-
intramolecular distances, as well as of H-O-H bond an- ble on the basis of external information. Note that the
gles, in Figure 2. Clearly, RMC_POT with the SPC/Ef ’RMC_POT random’ calculation produced gij(r)-s that
potential parameters leads to a molecular geometry that wereindistinguishablefromthosethathadresultedfrom
israthersimilartothatobtainedfrommoleculardynam- the other RMC_POT simulations. The intermolecular
ics simulations (with the SPC/E potential), although parts of each PPCF are identical for the 4 RMC_POT
there are slight differences: the average H···H non- calculations reported here (independently of the flexibil-
bondingdistanceissomewhat(bycca. 0.02Å)shorteras ity of the actual potential function used.).
aresultofRMC_POT,andasaresult,theaveragebond WhatisworthemphasizingisthatRMC_POTresults
angleisveryslightlysmaller. Themeanbondanglefrom from modeling the ND and XRD1 data sets have re-
6

50 15
(b)
(a) 3 MD-rigid g(r) Soper
MD-flexible gOO Skinner et al.
40 FNC_X1
10 (r)O2 FFNNCC__Xg_2X 1
(r)OH30 g(r)HH gO1 FNC_g_X2
g
20 (a)
5 0

10
0 0.9 1.0 1.1 01.5 1.6 1.7 1.8 (r)H1
O
g
r [¯] r [¯]
20
(c) H-O-H POT_X1 (b)
0 POT_X2
15 MD-rigid POT_r
) MD-flexible POT_X1 POT_Fw
os FNC_X1 POT_X2
P (c1H-O-H 0 FFFNNNCCC___Xgg__2XX 12 PPOOTT__rFw g(r)HH1
5
(c)
0
1 2 3 4 5
0
-0.6 -0.4 -0.2 0.0 r [¯]
cos
FIG.2. IntramolecularO-Hbonding(a),H···Hnon-bonding FIG.3. Partialpaircorrelationfunctionsforeachcalculation,
(b)distancesandthedistributionofthecosinesofintramolec- together with literature data from Soper11 and from Skinner
ular H-O-Hangles (c). et al.6.
produced the ’experimental’ gOO(r) that belongs to the
TABLE IV. Positions of intramolecular O-H and H-H, and
XRD2 data6. This is yet another sign of that data pre-
of the first intermolecular O-O and O-H maxima, as deter-
sentedbySkinneretal.6areofhighqualityandthatthey
mined from thePPCF-s, for the different simulations (in Å).
may, indeed, be called as ’consensual’(although perhaps
Note that ’RMC-FNC’ calculations have not provided any
not quite as ’benchmark’). well distinguishable maximum for the hydrogen bonding dis-
Distributions of intermolecular O···O···O, H-O···O tance (around 1.8 Å).
and (c) O-H···O angles are shown in Figure 5: these
Simulation O-H H-H O···O O···H
are all characteristic to the local environment (includ-
intra intra inter inter
inghydrogenbonding)ofwatermolecules. Althoughthe
MD-rigid 1.00 1.63 2.76 1.77
mainfeatures,thetetrahedrallocationofneighboringwa-
MD-flexible 1.03 1.66 2.73 1.72
ter molecules and the approximately straight hydrogen
FNC_X1 1.00 1.58 2.9 2.6
bonds, may be detected for each RMC calculations, the
FNC_X2 1.00 1.58 2.8 2.6
extents of these vary significantly. RMC_POT results
FNC_g_X1 1.00 1.64 2.75 1.79
follow characteristics of the original molecular dynamics
FNC_g_X2 1.00 1.64 2.76 1.79
simulationnearlywithinthe linethicknessinthefigures,
POT_X1 1.00 1.63 2.81 1.82
whereasitisratherhardtodetectthefeaturesinquestion
POT_X2 1.00 1.63 2.81 1.82
oncurvesresultingfromRMC-FNCcalculations. (Again,
POT_r 1.00 1.63 2.80 1.82
’RMC_POTrandom’resultswereindistinguishablefrom
POT_Fw 1.02 1.66 2.79 1.78
thoseoftheotherRMC_POTsimulations;also,theorig-
Soper11 2.77 1.84
inally rigid SPC/Ef and the flexible SPC/Fw potentials
Skinner6 2.80
produce identical curves.)
This behaviorof RMC_POT results is very encourag-
7

POT_X1
1.0
POT_X2 (a)
3
POT_r
(a) )
POT_Fw s
co
) 2 P (0.5 O...O...O
r
(O
O
g

0.0
1
(b)
g(r) Soper 3 H-O...O
)
gOO Skinner et al. s
o
02.5 2.6 2.7 2.8 2.9 3.0 P (c 2
MD-rigid
1
MD-flexible
FNC_X1 (b)
0
FNC_X2
10 O-H...O
FNC_g_X1
(c)
FNC_g_X2 MD-rigid MD-flexible
()grOH1 P (cos ) 5 FFNNCC__XX12 PPOOTT__XX12
FNC_g_X1 POT_r
FNC_g_X2 POT_Fw
0
-1 0 1
0
1.4 1.5 1.6 1.7 1.8 1.9 2.0 cos
r [¯]
FIG. 4. O-O and O-H partial pair correlation functions for FIG. 5. Distributions of the cosines of (a) O···O···O, (b)
each simulation, together with literature data from Soper11 H-O···O, and (c) O-H···O angles.
and from Skinner et al.6; the focus is on the regions of the
first intermolecular maxima.
IV. CONCLUSIONS
WehavepresentedresultsfromthefirstReverseMonte
Carlo study of liquid water in which a popular inter-
atomic potential model, SPC/E39, was applied explic-
ing: there seems to be no need for inventing extensive itly. After some ’learning period’, the calculations could
coordinationconstraints(c.f. Ref. 5)inRMCanylonger be tuned so that the RMC_POT algorithm was run-
if we wish to obtain physically meaningful particle ar- ningwiththe sameefficiencyasusualformodelingstud-
rangements for liquid water. Interestingly, RMC_POT ies without potentials; we therefore expect the exten-
provided slightly more regular tetrahedral local environ- sive exploitation of this approach for water (with var-
ment than even the originalMD (cf. Figure 5, part (a)). ious input diffraction data sets and/or thermodynamic
conditions) and aqueous solutions. Despite all of its at-
tractive features, it has to be made clear that although
A few words may be appropriate for a brief compari- RMC_POT applies interatomic potentials explicitly, it
son with EPSR results on liquid water11,31: there is a still is a method of structural modeling and cannot be
general agreement between RMC_POT and EPSR in used to reliably calculate properties beyond structural
terms of the main characteristics, although some of the ones.
details may appear quite differently. Perhaps the most The set of diffraction data applied here allowed the
apparent of these differences concern the intramolecular formation of local structural motifs in the RMC_POT
structure, for which RMC_POT seems to produce con- particle configurations that are hardly distinguishable
siderably narrower distributions for the O-H and H-H from those produced by molecular dynamics simulations
distances. The main reason behind may be the different that use the same (SPC/E) water potential (cf. Fig-
handling of intramolecular contributions, cf. Ref.31. ure 5). On the other hand, clear differences between
8
MD and RMC_POT findings were detected in terms of 12A. Zeidler, P. S. Salmon, H. E. Fischer, J. C. Neue-
total structure factors and partial pair correlation func- feind, J. Mike Simonson, and T. E. Markland,
tions. The very regular ’V-shape’ of water molecules, J.Phys.: Condens.Matter24,284126(2012).
13A. H. Narten, W. E. Thiessen, and L. Blum,
the tetrahedralenvironmentandstraighthydrogenbond
Science217,1033(1982).
angles reflect the present ’collective wisdom’ concerning 14S. Myneni, Y. Luo, L. Å. Näslund, M. Cavalleri, L. Ojamäe,
the structure of liquid water: it is therefore rather com- H. Ogasawara, A. Pelmenschikov, P. Wernet, P. Väterlein,
forting that the present study is able to offer large 3D C. Heske, Z. Hussain, L. G. M. Pettersson, and A. Nilsson,
particle configurationsthat, in addition, arefully consis- J.Phys.: Condens.Matter14,L213(2002).
tentwithdiffractiondata. Thesestatementshaveproven 15J. D. Smith, C. D. Cappa, K. R. Wilson, B. M. Messer, R. C.
Cohen, andR.J.Saykally,Science306,851(2004).
to be valid independently of the initial particle configu-
16T. Tokushima, Y. Harada, O. Takahashi, Y. Senba,
ration, as well as of the actual variant (rigid or flexible) H. Ohashi, L. G. M. Pettersson, A. Nilsson, and S. Shin,
of the SPC/E force field. Chem.Phys.Lett.460,387(2008).
Eventhoughreference(FNC-based)RMCcalculations 17O. Fuchs, M. Zharnikov, L. Weinhardt, M. Blum,
M. Weigand, Y. Zubavichus, M. Bär, F. Maier, J. D.
reproduced experimental diffraction data at the same
Denlinger, C. Heske, M. Grunze, and E. Umbach,
(very high) level as RMC_POT, and much better than Phys.Rev.Lett.100,027801 (2008).
moleculardynamicssimulations,the outcomefromthese 18A. Zeidler, P. S. Salmon, H. E. Fischer, J. C. Neuefeind,
runs is less attractive: the molecular shapes are less J. M. Simonson, H. Lemmel, H. Rauch, and T. E. Markland,
uniform and hydrogen bond angles take values that are Phys.Rev.Lett.107,145501 (2011).
19L. Temleitner, A. Stunault, G. J. Cuello, and L. Pusztai,
sometimes far from the ideal 180 degrees. Still, these ar- Phys.Rev.B92,014201(2015).
rangementsarefullyconsistentwithdiffractiondatacon- 20J.A.BarkerandR.O.Watts,Chem.Phys.Lett.3,144(1969).
sidered here – they just tend to be the most disordered 21B.Guillot,J.Mol.Liq.101,219(2002).
ones of those that may be called ’realistic’ (cf. Ref. 35). 22K. Szalewicz, C. Leforestier, and A. van der Avoird,
Chem.Phys.Lett.482,1(2009).
23W. L. Jorgensen and J. Tirado-Rives,
Proc.Natl.Acad.Sci.U.S.A.102,6665(2005).
ACKNOWLEDGMENTS 24R.L.McGreevyandL.Pusztai,Mol.Simul.1,359(1988).
25M. G. Tucker, M. T. Dove, and D. A. Keen,
J.Appl.Crystallogr. 34,630(2001).
The authors acknowledge financial support from the 26G.EvrardandL.Pusztai,J.Phys.: Condens.Matter17,S1(2005).
National Research, Development and Innovation Office 27O.Gereben, P. Jóvári, L. Temleitner, and L. Pusztai, J. Opto-
electron. Adv.Mater. 9,3021(2007).
of Hungary (NKFIH), under grant No. SNN 116198.
28G.Opletal,T.Petersen,B.O’Malley,I.Snook,D.G.McCulloch,
N.A.Marks, andI.Yarovsky, Mol.Simul.28,927(2002).
1Seee.g.http://www1.lsbu.ac.uk/water/. 29G. Opletal, T. Petersen, B. O’Malley, I. Snook, D. McCulloch,
2P. Wernet, D. Nordlund, U. Bergmann, M. Cavalleri, andI.Yarovsky,Comput.Phys.Commun.178,777(2008).
M. Odelius, H. Ogasawara, L. Å. Näslund, T. K. Hirsch, 30A.K.Soper,Chem.Phys.202,295(1996).
L. Ojamäe, P. Glatzel, L. G. M. Pettersson, and A. Nilsson, 31A.K.Soper,Phys.Rev.B72,104204(2005).
Science304,995(2004). 32O.Gereben andL.Pusztai,J.Comput.Chem.33,2285(2012).
3T. Head-Gordon and M. E. Johnson, 33M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith,
Proc.Natl.Acad.Sci.U.S.A.103,7973(2006). B.Hess, andE.Lindahl,SoftwareX 1-2,19(2015).
4M. Leetmaa, K. T. Wikfeldt, M. P. Ljungberg, M. Odelius, 34Z. Steinczinger and L. Pusztai,
J. Swenson, A. Nilsson, and L. G. M. Pettersson, Condens.MatterPhys.15,23606(2012).
J.Chem.Phys.129,084502(2008). 35I.Pethes andL.Pusztai,J.Mol.Liq.212,111(2015).
5L.Pusztai,Phys.Rev.B60,11851(1999). 36L. Pusztai, O. Pizio, and S. Sokolowski,
6L. B. Skinner, C. Huang, D. Schlesinger, L. G. M. J.Chem.Phys.129,184103(2008).
Pettersson, A. Nilsson, and C. J. Benmore, 37Z. Steinczinger and L. Pusztai,
J.Chem.Phys.138,074506(2013). Condens.MatterPhys.16,43604(2013).
7G. Hura, D. Russo, R. M.Glaeser, T. Head-Gordon, M. Krack, 38A. K. Soper, F. Bruni, and M. A. Ricci,
andM.Parrinello,Phys.Chem.Chem.Phys.5,1981(2003). J.Chem.Phys.106,247(1997).
8A.H.NartenandH.A.Levy,J.Chem.Phys.55,2263(1971). 39H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma,
9C. Huang, K. T. Wikfeldt, T. Tokushima, D. Nordlund, J.Phys.Chem.91,6269(1987).
Y.Harada,U.Bergmann,M.Niebuhr,T.M.Weiss,Y.Horikawa, 40H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren,
M.Leetmaa,M.P.Ljungberg,O.Takahashi,A.Lenz,L.Ojamae, A.DiNola, andJ.R.Haak,J.Chem.Phys.81,3684(1984).
A.P.Lyubartsev, S.Shin,L.G.M.Pettersson, andA.Nilsson, 41Y. Wu, H. L. Tepper, and G. A. Voth,
Proc.Natl.Acad.Sci.U.S.A.106,15214(2009). J.Chem.Phys.124,024503(2006).
10G.N.I.Clark,G.L.Hura,J.Teixeira,A.K.Soper, andT.Head- 42R.L.McGreevy, J.Phys.: Condens.Matter 13,R877(2001).
Gordon, Proc.Natl.Acad.Sci.U.S.A.107,14003(2010). 43O. Teleman, B. Jönsson, and S. Engström,
11A.K.Soper,J.Phys.Chem.B119,9244(2015). Mol.Phys.60,193(1987).