A classical reactive potential for molecular clusters of sulphuric acid and water

3.6 MB·

by  Jake L. Stinson

A classical reactive potential for molecular clusters of sulphuric acid and water
Jake L. Stinsona, Shawn M. Kathmannb and Ian J. Forda∗
aDepartment of Physics and Astronomy and London Centre for Nanotechnology,
University College London, Gower Street, London, WC1E 6BT, United Kingdom and
bPhysical Sciences Division, Pacific Northwest National Laboratory, Richland, Washington 99352, United States
We present a two-state empirical valence bond (EVB) potential describing interactions between
sulphuric acid and water molecules and designed to model proton transfer between them within a
classical dynamical framework. The potential has been developed in order to study the properties
of molecular clusters of these species, which are thought to be relevant to atmospheric aerosol
nucleation. The particle swarm optimisation method has been used to fit the parameters of the
6
EVB model to density functional theory (DFT) calculations. Features of the parametrised model
1
and DFT data are compared and found to be in satisfactory agreement. In particular, it is found
0
thatasinglesulphuricacidmoleculewilldonateaprotonwhenclusteredwithfourwatermolecules
2
at 300 K and that this threshold is temperature dependent.

n
a
J I. INTRODUCTION consider in this paper, for example, Kakizaki et al. [13]

1 concluded that nuclear zero point motion in clusters of
1 sulphuricacidandwaterat250Kgivesrisetonoticeably
It has long been suspected that sulphuric acid plays increased fluctuations and liquid-like behaviour. Sug-
] animportantroleinatmosphericparticlenucleationasa
h awaraet al. [14]studiedthedegreeofhydrationrequired
consequence of its affinity to water and its low volatility
p forthedissociationofthesulphuricacidmolecule,atthe
[1], thoughammoniaandorganicspecies, aswellasions,
- same level of theory [15]. Further evidence of small but
m are also likely to participate [2]. Considerable experi-
significant effects of zero point motion in similar clusters
mentaladvancesincharacterisingnucleationphenomena
e was provided by Stinson et al. [16].
h in atmospherically relevant conditions, with particular
However, whencomputingfreeenergiesusingab initio
c emphasisonsulphuricacid, havebeenreportedinrecent
methods, the harmonic oscillator approximation is com-
s. years[3,4]. Theinterpretationofsuchexperimentaldata
monly employed, based on identifying the lowest energy
c inordertounderstandbehaviouroverwiderconditionsis
i clusterconfigurationatzerotemperatureandestimating
s critically dependent on calculations of the free energies
the vibrational entropic contributions to the free energy
y of clusters of the nucleating species. Such thermody-
h from a characterisation of the low temperature normal
namic information is employed within a well-established
p modespectrum. Suchanapproachislikelytobeaccurate
theoreticalframeworkofclustergrowthanddecaytopre-
[ attemperatureswheretheclusterbehavesasavibrating
dict rates of formation of stable clusters [5, 6]. Although
solid-like structural network, but would seem to be less
1 further processes such as cluster coalescence or removal
appropriate for liquid-like systems. A more general ap-
v also play a role in particle formation [7], this kinetic and
1 proachisthentypicallyrequired,suchasthermodynamic
thermodynamicframeworkliesattheheartofourunder-
0 integration[17]wherethefreeenergyofasystemiscom-
standing of the phenomenon.
4 pared with that of a better understood reference system
2 Cluster free energies, however, are not straightforward at the same temperature, through performing a series of
0 to calculate. Simple models of cluster thermodynamic canonical averages with interpolating Hamiltonians, of-
.
1 properties have repeatedly been sought, but often such tenusingMonteCarlo(MC)methods. Approachesbased
0 approaches rely on an extrapolation of the properties of on nonequilibrium molecular dynamics (MD) simulation
6 larger droplets down to clusters consisting of only a few have also received attention [18, 19]. Methods of free en-
1 molecules [8–10]. The prime example of such a model is ergy computation are indeed rather numerous [20]. In
:
v classical nucleation theory (CNT) [5, 11], based on the many systems of interest, however, ab initio energy cal-
i capillarityapproximation. Cautionisnecessarywhenus- culations are too expensive, and modelling necessarily
X
ing such models if they suggest that the nucleation rate proceeds on a more coarse grained, classical level.
r is sensitive to the properties of very small clusters, since
a A number of classical schemes have been employed to
accuracy is not to be expected [12], though surprisingly
study clusters of sulphuric acid and water molecules. In
it is often observed.
an early study Kusaka et al. [21] developed a grand
More advanced models of molecular clusters can be canonical MC model based on rigid molecules and con-
used, of course, though at greater cost in computational cludedthattheclustersarehighlynon-sphericalandthat
effort. At the most fundamental level, molecular inter- bulk-like behaviour only emerges when there are at least
actionmodelsbasedonab initio quantummechanicsare 240 water molecules and 1−3 molecules of sulphuric
available. Indeed, very detailed studies have been con- acid,oritsdissociationproductbisulphate,inthecluster.
ducted, even taking account of the quantum nature of Later, Kathmann and Hale [22] presented a model based
the lighter nuclei present [13, 14] in addition to that of on rigid water and sulphate molecules and a free Hδ+
theelectrons. Inastructuralstudyofthesystemthatwe ion, treatedwithinanMCapproach. Dinget al. [23]de-
2
velopedaflexiblebondingmodelforsulphuricacid,bisul- for water with excess protons. A recent summary of the
phate,hydroniumandwaterspecies,whichwasemployed EVB model in the context of modelling proton transfer
byToivolaetal. [24]tostudyaquasiplanarliquid-vapour in a water network is given by Knight and Voth [33].
interfaceusingclustersof2000molecules. Amongstother We describe the basis of an EVB model in Section II.
matters, they found that when the sulphuric acid mole Section III is an account of the specific elements of the
fraction is less than 0.1, the acid molecules lie at the model developed here for sulphuric acid and water mix-
cluster surface and that the cluster structure is strongly tures. A demonstration of the model is given in Section
dependent on the number of bisulphate ions present. IV and in Section V we summarise and discuss future
However,althoughmixturesofdissociatedspecieswere applications. TheAppendixdescribessometechnicalas-
studied by Toivola et al. [24], the dynamics of the trans- pects of the particle swarm optimisation (PSO) fitting
fer of protons between species was not considered in the procedure used to parametrise the model against ab ini-
model. It is important to note that a classical potential tio calculations at the DFT level of theory.
suitable for free energy computations for clusters of sul-
phuric acid and water ought to be designed to describe
thetransferofaprotonfromtheacidtoawatermolecule. II. THE EVB METHOD
Withoutsuchacapability,itisnotpossibletoallowpro-
ton transfer to take place naturally; instead, the degree A. EVB overview
ofdissociation would have tobefixed by handin agiven
simulation. In an EVB model the potential energy of the reactive
This brings us to the aim of this study, which is to de- systemisconstructedasaninterpolationbetweenenergy
velopandparametriseaclassicalpotentialthatdoesoffer surfaces defined for specific choices of bonding pattern.
this capability. Such a potential would allow us to test This is done in a fashion reminiscent of the quantum
the harmonic approximation in calculations of the free mechanics of a multilevel system, though it should be
energyinsulphuricacid/waterclustersatrelevantatmo- emphasised that the dynamics considered are entirely
spheric temperatures. Using a classical potential that classical. The system is represented mathematically as
can describe the dynamics beyond structural vibrations a superposition of basis states, each representing a pos-
is a crucial requirement. sible bonding pattern. For example, Figure 1 illustrates
There are a few classical schemes available which al- two basis states of a system correspondingto the possib-
low reactions to occur. These include the Gaussian Ap- ilityinthiscasethatthecentralhydrogenatomcanform
proximation Potential (GAP) [25] where a potential en- part of either a sulphuric acid molecule on the left or a
ergysurfaceisconstructedbyfittingaGaussianbasisset hydronium ion on the right.
to reference data following a Bayesian statistics proced- Once the basis states of the system are specified, a
ure; the ReaxFF approach [26] which uses the bond or- matrixHthatresemblesaHamiltonianisconstructedas
dermethodology; andtheempiricalvalencebond(EVB) follows:
model [27] based on a superposition of underlying clas-
sical states of the system.
H =(cid:104)i|H|j(cid:105), (1)
The EVB methodology is attractive for use in the sul- ij
phuric acid/water system since, as we shall see, it is where (cid:104)i|H|i(cid:105) is the potential energy of the system ac-
based on classical potentials for each species plus mix- cording to the bonding pattern defined by |i(cid:105) and an as-
ingterms,andisrelativelystraightforwardtoimplement sociated classical potential. (cid:104)i|H|j(cid:105), where i (cid:54)= j, rep-
into existing MC or MD schemes. The methodology was resents the coupling between basis states |i(cid:105) and |j(cid:105) re-
introduced by Aqvist and Warshel [28] and is reviewed sponsible for mixing. In practice, (cid:104)i|H|j(cid:105) is specified by
by Kamerlin and Warshel [27]. It was first developed an empirical function chosen to reproduce the behaviour
in order to model proton transfer between hydronium of a higher level theoretical model.
and water species. Schmitt and Voth [29] designed a Once the H matrix has been constructed, the ei-
so-called multi-state EVB (MS-EVB) model for the sim- genvectors and eigenvalues are found. The eigenvector
ulation of systems of water molecules containing excess |Ψ(cid:105)=(cid:80) c |i(cid:105) with the lowest eigenvalue represents the
i i
protons. This work was further developed into the MS- ground state of the system, with coefficients c specify-
i
EVB2 [30] and the MS-EVB3 [31] models. Our strategy ingtheappropriatesuperpositionofthebasisstates. The
is to use the MS-EVB3 model as a framework for con- energyisgivenbyE =cTHcandtheforcescanbecom-
structinganEVBmodelforthesulphuricacidandwater puted via the Hellmann-Feynman theorem
system. As a simplification, we limit the multiplicity of
statestotwo(i.e. theprotonmaybeattachedeithertoa
bisulphate or to the nearest appropriate water molecule) ∂H (cid:88) ∂(cid:104)i|H|j(cid:105)
F =−(cid:104)Ψ| |Ψ(cid:105)=− c c , (2)
making it a two level EVB approach. Multiple poten- n ∂xn i j ∂xn
i,j
tial proton transfers within the system are accommod-
ated by employing the self-consistent iterative MS-EVB where x and F indicate the position and the force,
n n
(SCI-MS-EVB)model,developedbyWangandVoth[32] respectively, for atom n.
3
Figure 1. Illustration of two possible basis states of bonding for the same atomic configuration. State |1(cid:105) is composed of a
sulphuric acid and a water molecule whereas state |2(cid:105) contains bisulphate and hydronium ions. The EVB approach provides
an interpolation between the energy surfaces appropriate to each case.
B. The MS-EVB3 model of water systems with an
Table I. Parameters specifying the off-diagonal matrix ele-
excess proton
ments of the Hamiltonian used in the MS-EVB3 model of
water-hydronium transfer [31]; α, β and γ from [31], all oth-
The MS-EVB3 model [31] was originally designed for ers from [29].
simulating systems of water molecules with one excess
Vij -23.1871874 kcal/mol β 4.5 Å−1
proton. The key process of interest is the mobility of the const
excessprotonacrosshydrogenbondsformedbetweenwa- γ 1.85 Å−2 RO0O 3.1 Å
ter and hydronium. Diagonal elements of the matrix H P 0.2327246 P(cid:48) 10.8831327
are specified by a suitable potential for water and hy- k 9.562153 Å−2 α 15 Å−1
dronium bonding patterns. Off-diagonal components of D 2.94 Å r0 1.8136426 Å
OO OO
H are calculated in the following way:
(cid:104)i|H|j(cid:105)i(cid:54)=j =(Vciojnst+Veixj)·A(ROO,q), (3) whereγ,P,k,DOO,β,RO0O,P(cid:48),αandrO0Oareempirical
parameters(seeTableI)andq= 1(r +r )−r where
where R is the oxygen-oxygen distance in the hydro- 2 O1 O2 H
OO r ,r andr arethepositionsofthetwooxygensand
gen bond that includes the transferring proton, Vij is O1 O2 H
const hydrogen, respectively, that participate in the hydrogen
a constant and Vij is a representation of Coulombic in-
ex bond. The model is somewhat empirical in construction
teractionsbetweentheH O+ Zundelcationconsistingof
5 2 but was designed to have a number of desirable features
the hydronium and water between which the proton is
[31].
considered to be shared, and the remaining waters. It is
given by
C. The SCI-MS-EVB procedure
Vij =(cid:88)7 NH(cid:88)2O−1(cid:88)3 qnHk2Oqmex, (4) TheEVBmethodhasthecapacitytomodelsystemsof
ex R n > 1 hydroniums (excess protons) with water, but this
m k nk mnk extension increases the size of H to order mn where m is
where m is a label for the seven atoms in the Zundel the size of the basis set for a single proton system. The
cation, k is a label for the remaining water molecules SCI-MS-EVB model was developed by Wang and Voth
(N is the total number of water molecules in the sys- [32] to improve the scaling of the procedure with respect
H2O
tem) and n labels the three constituent atoms of each to n. This is achieved by treating the system as n single
k
watermolecule. qH2O arethepartialchargesforthecon- excess protons in the presence of an environment with
stituent atoms inntkhe water molecule and qex are partial a given state of bond mixing. Employing a single pro-
m
charges describing the Zundel cation. R is the dis- ton EVB methodology, the state of mixing of the entire
mnk
tance between atoms labelled by m and n . system is then given by a set of n eigenvectors each of
k
The function A(R ,q) has the form length m. The eigenvectors are corrected iteratively as
OO
the effects of the other excess protons on the Hamilto-
nianmatrixforeachoftheZundelcationsaretakeninto
A(R ,q)=exp(−γq2)·(cid:8)1+P exp(cid:2)−k(R −D )2(cid:3)(cid:9) account. Essentially, the environment affects how the
OO OO OO
competing energy surfaces are specified. The problem of
×(cid:8)1(cid:8)1−tanh(cid:2)β(R −R0 )(cid:3)(cid:9)+ analysing one matrix of order mn is hence reduced to
2 OO OO (5)the consideration of n matrices of size m. The number
ofiterationsrequiredtoreachconvergenceinenergyand
P(cid:48)exp(cid:2)−α(ROO−rO0O)(cid:3)(cid:9), forces is typically small. We similarly use this procedure
for our system.
4
III. EVB MODEL FOR SULPHURIC ACID AND viour of sulphuric acid in clusters containing relatively
WATER few water molecules.
Inthefollowingsubsectionswedescribeinsomedetail
the construction of an EVB model for a system contain- B. Algorithms
ingsulphuricacid,bisulphate,hydroniumandwaterspe-
cies. Itisacomplexnarrative,butasummaryofthepro- A central part of the EVB methodology is the con-
cedure is provided at the end in Subsection IIIG. In the structionofappropriatebasisstatesofthesystem,which
discussionanamingconventionisusedwheretheground is nontrivial as reactions can cause atoms to associate
statereferstothebondingpatternthathasbeenidenti- with different molecules in the course of a simulation.
fied as geometrically the most ‘natural’, based upon the This section describes two algorithms which are applied
configurationoftheatoms. Excited statesarebonding at every time step in order to identify the ground state
configurationsthatareidenticaltotheground stateex- andexcitedstatesofthesystembygeometricarguments.
ceptfortherepositioningofonebondasaconsequenceof
a proton transfer. Therefore, it is possible to identify an
excited state based upon the ground state and the 1. Ground state selector algorithm
two species involved in a proton transfer. The species
are then referred to as a donor (the molecule to which
An algorithm has been designed which constructs the
the hydrogen atom is bonded in the ground state and
ground state from a list of atomic positions. For con-
which can either be a sulphuric acid molecule or a hy-
venience,fourmolecularlistsaredefinedandaredenoted
dronium ion), and an acceptor (the molecule to which
by SA (sulphuric acid), BS (bisulphate), HY (hydronium)
the hydrogen atom is bonded in the excited state and
andWA(water)andanassignmentreferstotheidentifica-
which can either be a bisulphate or a water). For clar-
tionofanatomasaconstituentofaparticularmolecule.
ity,thetermsneutralandionisedareusedtoidentifythe
The algorithm is performed in the following way:
state of the sulphur-bearing species (either as a neutral
sulphuric acid or an ionised bisulphate). 1. The four oxygens closest to each sulphur atom are
identified. They are assigned to a sulphuric acid
molecule(alongwiththesulphuratom)andplaced
A. Basis set size in the SA list.
2. The remaining oxygen atoms are assigned to hy-
The basis set size m for each transferable proton was
dronium ions and added to the HY list.
chosen to be two, in order to keep the model as simple
as possible, in line with the SCI-MS-EVB approach im-
3. Minimumnumbersofhydrogenatomsareassigned
plemented for water and described in Section IIC. This
to each molecular species. For each member of the
approach allows each donor and acceptor species in the
HY list the two hydrogen atoms closest to the asso-
systemtobeinvolvedinamaximumofoneprotontrans-
ciatedoxygenareidentifiedandassignedtothelist.
fer in a given configuration, but it is assumed that the
The same process is used to assign one hydrogen
possibility that such a species might be involved in two
to each member of the SA list.
proton transfers simultaneously can be safely neglected.
The small basis set leads to occasional ambiguities con- 4. The remaining hydrogen atoms are placed (in no
cerning which proton is considered to be shared but a particular order) in a temporary list named H. The
modification of the off-diagonal term specified in Eq. (5) closest oxygen to each of these hydrogens is iden-
has been made to limit the impact of these issues and is tified and the hydrogen is assigned to the molecu-
described in Section IIID1. lar species of which the oxygen is a constituent.
Our focus here has been to develop the simplest pos- The oxygen cannot be part of a molecule which
sible scheme of the process of proton transfer between hasalreadyacceptedoneoftheseremaininghydro-
sulphuricacidandwater,andforeaseofimplementation gens and cannot be a member of the SA list that
wehaveusedasimilartwo-statemodelforthewaternet- was assigned a bond to a hydrogen atom in step 3.
work. Largerbasissetshavepreviouslybeenemployedin If there is an attempt to assign a hydrogen where
studies of the water/hydronium system in order to cap- the bond length is greater than 1.2Å then this as-
ture effects such as the formation of highly correlated signment is rejected and the atom is moved to the
structures such as the H O + complex discussed in [33]. top of the H list. If the H list is rearranged at any
9 4
Ourapproachislessrefinedatpresent,butclearlyamore time,thenstep4isimmediatelyrestartedusingthe
elaboratemodelofstatemixingwithinthewaternetwork modifiedHlist. Thehydrogenswhichwentthrough
could be implemented at a later stage alongside the two- a rejected assignment are no longer subject to the
statemodelofacid/watermixinginordertocapturethe 1.2Å constraint, but if the bond length is over 2Å
finer details of proton transfer in the water sector. Our the hydrogen is again returned to the top of the H
principal focus has been to describe the reactive beha- list. If problems with assignment persist, then the
5
Figure 2. The ground state selector algorithm. Image zero shows the positions of atoms, and in subsequent steps bonds are
inserted according to the procedure described in Section IIIB1. Step three is shown in two parts; first, the attachment of two
hydrogen atoms to each oxygen and second, the attachment of one hydrogen atom to each sulphate.
algorithmisrunwiththebondlengthcheckturned 2
off. In practice, there is only need for these checks
in the circumstances of proton transfer events, in
3
which case the two state EVB mixing between the
states will ensure that the appropriate forces are 1
appliedregardlessoftheassignmentofthemolecu-
lar species. This step is completed when all hydro-
gensintheHlisthavebeenassigned. Theprocedure
followed in this step was found to work smoothly
in trials of the model.
5. AnysulphuricacidmoleculeintheSAlistforwhich
onlyonehydrogenatomhasbeenassignedismoved
to the BS list. Similarly, any hydronium molecule
Figure 3. A simultaneous consideration of three potential
in the HY list that does not have three assigned
proton transfer events. The groups labelled 1, 2 and 3 refer
hydrogen atoms is moved to the WA list.
to [hydronium/bisulphate], [sulphuric acid/water] and [hy-
Figure2illustratesthestepsofthealgorithmastheyare dronium/water] proton transfer, respectively. One ground
performedonasystem. Attheendoftheprocedure, the state (the bonding pattern shown) and three excited
SA,BS,HYandWAlistsrepresenttheground stateofthe states (where the central hydrogen in each group is bonded
system. differently) have therefore been identified.
2. Excited state identification algorithm
no further acceptor or donor molecules for which
As stated earlier, an excited state is defined as a
it is possible to construct an O−H separation dis-
changeinthebondingassignmentofonehydrogeninthe
tance.
ground state of a given configuration. An algorithm
identifies the excited states in a system using the fol-
lowing steps:
1. The shortest distance between a hydrogen belong-
The ground state selector and the excited state iden-
ingto adonor species, andan oxygenbelongingto
tification algorithms produce one ground state and n
an acceptor species is identified. The oxygen in a
excited stateswherenisthenumberofpotentialpro-
bisulphate ion which has an attached hydrogen is
ton transfer events in the configuration under consider-
notconsideredhere. Ifthisdistanceislessthan2Å
ation. When n = 0 the system is in a non-reactive con-
thenan excited statemaybeconstructedwhere
figuration. In the configuration shown in Figure 3, there
thehydrogenisreassignedtotheacceptorspecies.
arethreegroupswherethebondingofaprotonisunclear
2. Step1isrepeateduntileithertheshortestdistance (ringed), such that a ground state and three excited
identified is greater than 2Å in length or there are states have been identified in this system.
6
C. Diagonal terms only come about when the proton in question occupies a
specific region between the donor and acceptor species.
Thepotentialsneededtospecifythediagonalelements Outside this region, the interactions revert to those of a
of H are those developed by Loukonen et al. [34], and single pattern of bonding, and crucially, the energy sur-
the SPC/EF potential is used for the water molecules face is continuous at the boundary of the region. This is
[35]. The hydronium energy is represented by a set of encoded in the specification of off-diagonal terms which
harmonic angle potentials and Morse bond potentials as vanish for proton positions outside an ellipsoidal volume
follows: centred on the mid-point between the donating and ac-
cepting oxygen atoms.
V = 1(cid:88)3 k (θ −θ )2+(cid:88)3 D(cid:104)1−e−α(rj−r0)(cid:105)2, Specifically, the factor in Eq. (5) involving q becomes
hyd 2 i=1 θ i 0 j=1 (cid:40)(cid:104)exp(−(γq20−1))−1(cid:105) q2 ≤γ−1
(6) exp(−γq2)→ e−1 0 (7)
0 q2 >γ−1
where i is summed over the three hydrogen-oxygen- 0
hydrogen angles and j over the three oxygen-hydrogen
where q2 is related to q in the following way. The q
bonds. The k and θ parameters are taken from 0
θ 0 vector is expressed in a new orthogonal coordinate sys-
Loukonen et al. [34]. The Morse potential is similar
tem q(cid:48) in which the x-axis is parallel to r , where
in form to the model used in Wu et al. [31], and em- OO
r = r −r is the vector separation between oxy-
ploys the same value of D, but the α parameter has OO O1 O2
gens. The other two axes take arbitrary orientations. q2
been changed so that the second order Taylor expansion 0
is then defined as q2 =(q(cid:48)/τ )2+(cid:0)q(cid:48)/τ (cid:1)2+(q(cid:48)/τ )2
around r0 matches the strength of the harmonic spring 0 x x y yz z yz
where the τ and τ are unitless scaling factors repres-
used to represent the OH bond in Loukonen et al. [34], x yz
−1 enting semi-major axes of an ellipsoid with respect to
namely α = 2.327Å . The use of a Morse potential
the x axis, and the y and z axes, respectively. In effect,
rather than a harmonic spring was found to improve the
a hydrogen atom experiences mixed bonding within an
matchbetweenourmodelofhydroniumandthatusedin
ellipsoidal volume between the two appropriate oxygen
the original MS-EVB3 model.
atoms. We find that this eliminates occasional ambigu-
One of the diagonal terms is augmented by an energy
ities with regard to which oxygen pair is appropriate for
shift,∆,toaccountforthedifferenceinzerotemperature
the mixed bonding experienced by a particular hydro-
ground state energy between the sulphuric acid and wa-
gen, introducing discontinuities in the potential energy
ter bonding pattern and the bisulphate and hydronium
landscape. The two variables τ and τ are treated as
version. Dingetal. [23]calculatedthevalueofthispara- x yz
additional fitting parameters. The form of Eq. (7) en-
meter to be 144.0kcal/mol (602.5kJ/mol), but for our
sures two features: the expression is equal to unity at
purposes it was considered to be a free parameter and is
q2 = 0 and goes to zero on the surface of the ellipsoid
fitted to higher level theory data. 0
at q2 = γ−1, beyond which the off-diagonal term van-
0
ishes to enforce ordinary unmixed behaviour outside the
ellipsoid.
D. Off-diagonal terms
In principle, a jump in energy can still occur if ellips-
oidal mixing regions between two donor-acceptor pairs
The form chosen for the off-diagonal terms in the come into overlap as a consequence of configurational
Hamiltonian describing sulphuric acid-water proton change. Suchsituationscanoccasionallyariseduringthe
transfer is based upon the MS-EVB3 model (Eq. (5)), moleculardynamics,butwehaverestrictedtheregionsof
withtwomodifications. Theserelatetotheexpressionin- mixing, through the parameters τ and τ , specifically
x yz
volvingtheqparameter(SectionIIID1)andthemethod to avoid this. For most of the time there is no ambiguity
for calculating Veixj (Section IIID2). in the choice of basis set to use for a configuration and
fewoccasionswhenajumpinenergyoccursaccordingto
our scheme. These jumps are unphysical, and a break-
1. The q dependence age of NVE conditions, but are the price to pay for the
simplicity of the implementation.
Thereisaproblemwithusingalimitedbasisset,espe-
ciallywhentherearelong-rangeoff-diagonaltermsinthe
Hamiltonian matrix. If configurational evolution brings 2. Interpolated charges
about a change in identity of the species that might ac-
cept the proton from a given donor species, then the The Vij parameter in Eq. (4) is also modified since
ex
replacement of one energy surface by another will bring thechargesprovidedbytheMS-EVB3modelforthecal-
about an abrupt jump in energy that will need to be culation of this term are specific to the Zundel cation
damped out by the thermostat. However, we have de- and do not transfer to a sulphuric acid and water model.
signedthetwo-stateEVBschemetoavoidsuchjumpsas We recall that it represents external electrostatic inter-
much as possible. Mixing between bonding patterns will actions for the group of atoms involved in the mixing
7
of bonding patterns. Our approach is to interpolate the acceptor/donor pair m in the ground state and
charges of these atoms according to the c coefficients of excited state of m respectively. The expression
i
the mixture. The procedure is as follows: Vij intheoff-diagonaltermisupdatedinthesame
ex
fashion;thisincludesrecalculatingthechargesqac-
1. Initiallythechargesareassumedtobeq = 1(qgs+
qes), where qgs and qes are the partial c2harges cordingtotheprocedureinSectionIIID2. Thecn
values are then recalculated. This allows the coef-
for that atom according to the classical potential
ficients c for the acceptor/donor pair n to be
n
provided, in accordance with whether the bonding
correctedfortheintermolecularinteractionsarising
is identified as the ground state or the excited
from acceptor/donor pair m. This step is re-
state, respectively.
peated until each c vector has been corrected for
n
each acceptor/donor pair m (where m(cid:54)=n).
2. A self-consistent iteration is performed where q is
updated according to q = c2 qgs + c2 qes where
gs es 3. Step 2 is repeated until all c eigenvectors have
c and c are the eigenvector coefficients for the n
gs es converged. This is tested by defining C =
ground state and excited state, respectively, (cid:80) (cold−cnew)2 wherecold andcnew arethseum(nor-
obtained from the two-level Hamiltonian describ- n n n n n
malised)c eigenvectorscalculatedbeforeandafter
ing the mixing. n
step 2 is performed. The system is considered to
This procedure is performed alongside a separate self- be converged when Csum < 10−5, or when step 2
consistentiterativeprocessdescribedinthenextsection, has been cycled 10 times. In practice only a few
thepurposeofwhichistoallowconsiderationofmultiple iterations are required.
proton transfers.
F. Energy and force calculation
E. Multiple proton procedure
Once the self-consistency iterations have been per-
The modelling of systems undergoing multiple proton formed, the energy of the system can be computed in
transfers follows the SCI-MS-EVB procedure developed the same fashion as in the SCI-MS-EVB method. Eq.
by Wang and Voth [32]. Figure 3 shows a system where (13) in reference [32] gives a deconstruction of the total
thisextensionisrequiredasithasthreeexcited states energy of a system containing two excess protons in the
available. The starting point in this procedure is to form
neglect overlap between two different excited states (i.e.
(cid:104)i|H|j(cid:105) = 0 where i and j represent two different ex-
cited states). This allows for the reduction of the EVB E =E +E +E +E +E +E , (8)
Hamiltonian matrix of the system, which is of order nnc, total AA BB AB AR BR RR
d
to n matrices of order n where n is the number of
d c d where E refers to the total energy of a system ac-
total
donor molecules being considered and n is the number
c cording to the EVB method. A and B refer to sep-
basis states considered per donor molecule. An iterative
arate acceptor/donor pairs (see Figure 1 in reference
procedure is then used to correct these matrices for the
[32]). E and E refer to the independent energy
AA BB
(cid:104)i|H|j(cid:105) contributions. The method can be described as
contributions of the acceptor/donor pairs and includes
follows:
their off-diagonal contributions. E describes the en-
AB
ergy contribution due to interactions between A and
1. We determine the eigenvector for each 2×2 mat-
B. R refers to the rest of the system, and there are
rix H corresponding to the ground state and an
three further contributions due to interactions between
excited state. For clarity, n is used as a la-
A and R, B and R and the independent energy con-
bel such that c and H are defined for the nth
n n tribution of R, referred to as E , E and E re-
excited state. In addition, n is used to label AR BR RR
spectively. Eq. (8) can be generalised to E =
the donor and acceptor species that define the total
(cid:16) (cid:17)
(cid:80) (cid:80)
excited state. ERR + n Enn+EnR+ m(cid:54)=nEnm where n and m
refer to the identified donor/acceptor pairs. Once
2. EachH isthencorrectedfortheeffectof excited
n E has been constructed it simply remains to apply
state m where m (cid:54)= n. This is performed total
the Hellmann-Feynman theorem to determine the force
by updating the intermolecular energy contribu-
acting upon each atom within the system.
tions to the diagonal terms in H , i.e. Egs →
(cid:0)c2 Egs +c2 Ees (cid:1) where c2 nand c2 ainrteerthe
gs inter es inter gs es
squared coefficients of the c eigenvector repres-
m
enting the associated ground state and excited G. Model overview
state. Egs and Ees are the intermolecu-
inter inter
lar energy contributions resulting from interac- In summary, the model is an algorithm consisting of
tions between the acceptor/donor pair n and the the following series of steps:
8
1. Construct a ground state from a set of atomic
positions (Section IIIB1).
2. Identify the excited states (Section IIIB2).
3. Calculatetheon-diagonalandoff-diagonaltermsof
the 2×2 matrix H for each excited state (Sec-
tions IIIC and IIID1).
4. Optimise the c vectors for each excited state
i
by performing a self-consistent iterative proced-
ure which revises the off-diagonal terms and al-
lowsmultipleprotontransferstobeaccommodated
(Sections IIID2 and IIIE).
5. Calculate the energy of the system and the forces
acting upon each atom (Section IIIF).
Once the parameters of the model have been selected,
the procedure can be implemented in a classical mo-
lecular dynamics code, and we have done this using the
DL_POLY 4.03 package [36].
IV. RESULTS
We have parametrised the EVB model against refer-
ence configurational energies calculated from a DFT ap-
proachusingaparticleswarmoptimisation(PSO)fitting
scheme, the details of which are described in the Ap-
pendix. We focus in this section on investigating various
features of the parametrised model.
Figure 4 shows the potential energy of a sulphuric
acid/watersystemasafunctionofthepositionofatrans-
ferring hydrogen atom, comparing the performance of
the EVB model with a set of reference data. It is im-
portant to recognise that the EVB scheme does not op-
erate for hydrogen positions beyond a distance ±0.3Å Figure 4. Sets of positions of a hydrogen atom in spatial
in the r direction starting from 1(r + r ), the regions of mixed bonding (a), together with corresponding
OO 2 O1 O2 contour plots of the potential energy in units of eV (b), for
mid-pointbetweentheoxygensatomsinvolvedinthehy-
the I-n configuration (labelled according to reference [37]) of
drogen bond, and so energies at hydrogen positions la-
asulphuricacidandawatermolecule. Thecentreofthecon-
belled±0.4Åinther directionareforunmixedbond-
OO tourplotsliesatthemid-pointbetweenthetwoparticipating
ing and are not affected by the EVB parametrisation
oxygens and the labels 1, 2 and 3 in (b) refer to the arrays
scheme. The comparison demonstrates an overall cor- of points shown in (a). The cyan, red and green spheres in
respondence between the EVB model and the DFT res- (a) match the equivalent top left, top right and the bottom
ults. The scheme is limited in that it is only active for right of each plot in (b), respectively. The energy when the
hydrogen positions lying in an ellipsoidal spatial region hydrogenisatposition(−0.4Å,−0.4Å,−0.4Å)relativetothe
around the mid-point between the oxygens, that it only mid-pointissettoareferencepointofzeroforboththeDFT
represents two possible bonding states per proton and and EVB calculations.
that ultimately it is underpinned by a classical potential
fit to the reference data and an empirical form for the
off-diagonal term. Nevertheless, the EVB model seems 1 ns after a 20 ps equilibration period. The simula-
to capture the shape of the energy surface experienced tions employed a Langevin thermostat with target tem-
by the hydrogen atom. peratures of 300K or 200 K using a modified version of
the DL_POLY 4.03 program [36]. There was no con-
Next, we demonstrate how the propensity for proton straint on the centre of mass motion or the rotation of
transfer from an acid molecule is affected by the extent the clusters. The results show that at 300 K there is
of the surrounding water network. Figure 5 tracks the very little propensity for the proton to transfer to the
population of water molecules in various simulations of closest water and remain there more than momentarily
[H SO ]+[H O] with n=2−4, over a time period of in the n=2 and n=3 cases. The hydrogen bond rarely
2 4 2 n
9
work performed at zero temperature using a higher level
of theory [15, 37, 38], according to which the first dis-
sociation event of a sulphuric acid molecule occurs when
hydrated by between 3 and 6 water molecules.
When running the EVB model at 300 K we ob-
servechangesinconfigurationthatinvolveintermolecular
bondmakingandbreaking,characteristicoftheexpected
liquid-likenatureofthesystematrelativelyhightemper-
atures. We also observe the Grotthuss mechanism where
protons shuffle around the water network [39].
V. CONCLUSIONS
A two-state EVB model for the sulphuric acid, bisul-
phate, hydronium and water system has been presen-
ted, which is based upon the MS-EVB3 and SCI-MS-
EVB models developed for hydronium/water alone [29,
30]. The approach essentially provides an interpolation
betweentwoenergysurfacescorrespondingtothechoices
of bonding pattern available to a proton in a particu-
lar spatial region between a bisulphate and a water. It
involves a Hamiltonian represented in a basis of bond-
ingpatterns,withdiagonaltermsbasedonclassicalforce
fields,andoff-diagonaltermsconstructedempiricallyand
fitted to higher level calculations.
The limitation of the basis to two states implies that
the model does not match the state of the art in the de-
scriptionofprotonatedwatersystems[33]butitprovides
the simplest scheme for representing proton transfer
betweensulphuricacidandwaterspecies,thefocusofour
attention here and a new application of the EVB frame-
work. An extension to include more elaborate basis sets
would be possible.
Theparametrisationoftheoff-diagonaltermshasbeen
performed using the particle swarm optimisation (PSO)
technique, with reference to DFT data for specific cases
of sulphuric acid/water proton transfer events. The res-
Figure 5. Evolution in the population of water mo- ulting scheme has been integrated into a modified ver-
lecules in simulations starting with the values in- sion of the DL_POLY 4.03 MD code [36]. A quantit-
dicated above each frame. A decrease of one in ative comparison of the potential energy surface for the
the water population indicates an ionisation event transferring proton within a geometry-optimised static
[H SO ]+[H O] →[HSO ]−+[H O]++[H O] , and
2 4 2 n 4 3 2 n−1 configuration of one sulphuric acid and one water mo-
an increase by one refers to the reverse of this reaction.
lecule gives good agreement. The level of hydration re-
The simulations were performed at 200 K (top) and 300 K
quired for the first ionisation or dissociation of a single
(bottom).
sulphuric acid molecule at 300K has been investigated,
revealingthatasystemwithfourwatermoleculesisbest
describedbyacomplexmixturebetweenneutralandion-
fluctuatestosuchanextentthatthegroundstateconfig- ised bonding. At 200K the four water case remains pre-
uration consists of a bisulphate/hydronium rather than dominantly neutral, indicating that an increase in tem-
the sulphuric acid/water bonding pattern. However, for peraturepromotesprotontransfer,aswouldbeexpected.
then=4case,thesystemspendsthemajorityofitstime In highly hydrated structures, the proton becomes very
in the dissociated state. Proton transfers are more fre- mobile within the water network.
quentandquasistable. Whenthen=4simulationisrun The main motivation for developing this model is to
at 200 K, however, we see that it remains undissociated provide a tool for the fast simulation of cluster struc-
for a longer period of time suggesting that the dissoci- tures. The key performance details are as follows. The
ated state is less stable at lower temperature, which is run time for a 1 ns simulation with six water molecules
tobeexpected. Theseconclusionsareinagreementwith and one acid is 430seconds on one Intel Core i5-460M
Molecular Physics 10
processor and it is, therefore, easily possible to gener- ACKNOWLEDGEMENTS
ate trajectories several tens of nanoseconds in length for
quite substantial systems using this model. The para- SMK was supported in part by the U.S. Department
metrisationprocedureandtestswehavereportedgiveus of Energy, Office of Science, Office of Basic Energy Sci-
confidence that the physics of proton transfer has been ences, Division of Chemical Sciences, Geosciences, and
well captured by the scheme. Our computationally in- Biosciences; JLS and IJF were supported by the IM-
expensive method is suitable for performing calculations PACT scheme at University College London (UCL). We
ofthethermodynamicpropertiesofclustersofwaterand acknowledge the UCL Legion High Performance Com-
sulphuric acid molecules at a range of temperatures rel- puting Facility, and associated support services together
evanttotheatmosphere. Weshallreporttheseresultsin with the resources of the National Energy Research Sci-
future publications. entific Computing Center (NERSC), which is suppor-
ted by the U.S. Department of Energy under Contract
No. DE-AC02- 05CH11231. JLS thanks Dr. Gregory
Schenter, Dr. Theo Kurtén and Prof. Hanna Vehkamäki
for important guidance and discussions.
[1] R. Zhang, A. Khalizov, L. Wang, M. Hu and W. Xu, [21] I.Kusaka,Z.G.WangandJ.H.Seinfeld,J.Chem.Phys.
Chem. Rev. 112, 1957 (2012). 108, 6829 (1998).
[2] M. Kulmala, T. Petäjä, M. Ehn, J. Thornton, M. Sip- [22] S.M. Kathmann and B.N. Hale, J. Phys. Chem. B 105,
alä, D. Worsnop and V.M. Kerminen, Annu. Rev. Phys. 11719 (2001).
Chem. 65, 21 (2014). [23] C.G. Ding, T. Taskila, K. Laasonen and A. Laaksonen,
[3] J. Kirkby et al., Nature 476, 429 (2011). Chem. Phys. 287, 7 (2003).
[4] F. Riccobono et al., Science 344, 717 (2014). [24] M. Toivola, I. Napari and H. Vehkamäki, Boreal Envir-
[5] I.J. Ford, Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. onment Research 14, 654 (2009).
Sci. 218, 883 (2004). [25] A.P.Bartók,M.C.PayneandG.Csányi,Phys.Rev.Lett.
[6] V.I. Kalikmanov, Nucleation Theory (Springer, Heidel- 104, 1 (2010).
berg, 2013). [26] A.C.T.vanDuin,S.Dasgupta,F.LorantandW.A.God-
[7] M.J.McGrath,T.Olenius,I.K.Ortega,V.Loukonen,P. dard, J. Phys. Chem. A 105, 9396 (2001).
Paasonen, T. Kurtén, M. Kulmala and H. Vehkamäki, [27] S.C.L.KamerlinandA.Warshel,WileyInterdisciplinary
Atmos. Chem. Phys. 12, 2345 (2012). Reviews: ComputationalMolecularScience1,30(2011).
[8] S.L. Girshick and C.P. Chiu, J. Chem. Phys. 93, 1273 [28] J. Aqvist and A. Warshel, Chem. Rev. 93, 2523 (1993).
(1990). [29] U.W. Schmitt and G.A. Voth, J. Phys. Chem. B 102,
[9] A. Laaksonen, I.J. Ford and M. Kulmala, Phys. Rev. E 5547 (1998).
49, 5517 (1994). [30] T.J.F. Day, A.V. Soudackov, M. Čuma, U.W. Schmitt
[10] V.I. Kalikmanov and M.E.H. van Dongen, J. Chem. and G.A. Voth, J. Chem. Phys. 117, 5839 (2002).
Phys. 103, 4250 (1995). [31] Y.Wu,H.Chen,F.Wang,F.PaesaniandG.A.Voth,J.
[11] H.Vehkamäki,ClassicalNucleationTheoryinMulticom- Phys. Chem. B 112, 467 (2008).
ponent Systems (Springer, Berlin, 2006). [32] F. Wang and G.A. Voth, J. Chem. Phys. 122, 144105
[12] D. Brus, K. Neitola, A.P. Hyvärinen, T. Petäjä, J. Van- (2005).
hanen, M. Sipilä, P. Paasonen, M. Kulmala and H. Li- [33] C. Knight and G.A. Voth, Acc. Chem. Res. 45, 101
havainen, Atmos. Chem. Phys. 11, 5277 (2011). (2012).
[13] A. Kakizaki, H. Motegi, T. Yoshikawa, T. Takayanagi, [34] V. Loukonen, T. Kurtén, I.K. Ortega, H. Vehkamäki,
M.ShigaandM.Tachikawa,J.Mol.Struc-THEOCHEM A.A.H. Pádua, K. Sellegri and M. Kulmala, Atmos.
901, 1 (2009). Chem. Phys. 10, 4961 (2010).
[14] S.Sugawara,T.Yoshikawa,T.Takayanagi,M.Shigaand [35] J. Lopez-Lemus, G.A. Chapela and J. Alejandre, J.
M. Tachikawa, J. Phys. Chem. A 115, 11486 (2011). Chem. Phys. 128, 174703 (2008).
[15] C. Arrouvel, V. Viossat and C. Minot, J. Mol. Struc- [36] I.T. Todorov, W. Smith, K. Trachenko and M.T. Dove,
THEOCHEM 718, 71 (2005). J. Materials Chem. 16, 1911 (2006).
[16] J.L. Stinson, S.M. Kathmann and I.J. Ford, J. Chem. [37] S.Re,Y.OsamuraandK.Morokuma,J.Phys.Chem.A
Phys. 140, 024306 (2014). 103, 3535 (1999).
[17] D. Frenkel and B. Smit, Understanding Molecular Sim- [38] A.R.BandyandJ.C.Ianni,J.Phys.Chem.A102,6533
ulations: from Algorithms to Application, 2nd ed. (El- (1998).
sevier, Amsterdam, 2001). [39] A.Y.Lozovoi,T.J.Sheppard,D.L.Pashov,J.J.Kohanoff
[18] I.J. Ford and S.A. Harris, J. Chem. Phys. 120, 4428 and A.T. Paxton, J. Chem. Phys. 141, 044504 (2014).
(2004). [40] A.Banks,J.VincentandC.Anyakoha,Nat.Comput.6,
[19] H.Y.TangandI.J.Ford,Phys.Rev.E91,023308(2015). 467 (2007).
[20] C. Chipot and A. Pohorille, (Eds.), Free Energy Calcu- [41] A.Banks,J.VincentandC.Anyakoha,Nat.Comput.7,
lations (Springer, Heidelberg, 2007). 109 (2007).
[42] R. Poli, J. Artif. Evol. App. 2008, 685175 (2008).