Bicanonical Ab Initio Molecular Dynamics for Open Systems - Journal

Article pubs.acs.org/JCTC

Bicanonical ab Initio Molecular Dynamics for Open Systems Johannes Frenzel,*,† Bernd Meyer,‡ and Dominik Marx† †

Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, 44780 Bochum, Germany Interdisziplinäres Zentrum für Molekulare Materialien (ICMM) and Computer-Chemie-Centrum (CCC), Friedrich-Alexander-Universität Erlangen-Nürnberg, 91052 Erlangen, Germany

‡

ABSTRACT: Performing ab initio molecular dynamics simulations of open systems, where the chemical potential rather than the number of both nuclei and electrons is fixed, still is a challenge. Here, drawing on bicanonical sampling ideas introduced two decades ago by Swope and Andersen [J. Chem. Phys. 1995, 102, 2851−2863] to calculate chemical potentials of liquids and solids, an ab initio simulation technique is devised, which introduces a fictitious dynamics of two superimposed but otherwise independent periodic systems including full electronic structure, such that either the chemical potential or the average fractional particle number of a specific chemical species can be kept constant. As proof of concept, we demonstrate that solvation free energies can be computed from these bicanonical ab initio simulations upon directly superimposing pure bulk water and the respective aqueous solution being the two limiting systems. The method is useful in many circumstances, for instance for studying heterogeneous catalytic processes taking place on surfaces where the chemical potential of reactants rather than their number is controlled and opens a pathway toward ab initio simulations at constant electrochemical potential. simulations,1,2 such as Widom’s insertion method or alchemical transformations, whereas our focus in what follows is exclusively on MD simulations at a predefined constant chemical potential or average fractional number of a specific chemical species. In the framework of ab initio molecular dynamics (AIMD) simulations,8 where the full electronic structure is propagated concurrently with the nuclei mostly in the framework of approximate density functional theory (DFT), constant temperature and constant pressure (stress) simulations are available (see e.g., sections 5.2.2 and 5.2.3, respectively, in ref 8). However, significant complications arise in AIMD at constant chemical potential due to the need to also propagate the electrons, which renders both particle exchange with a reservoir and dealing with fractional particles challenging (see e.g., section 5.3.4 in ref 8). Yet, AIMD at constant chemical potential of the electronic subsystem has been realized by introducing free energy density functionals and “fractional electrons” (more precisely fractional occupation numbers) in the framework of dynamical propagation schemes.9,10 A conceptually elegant alternative method to treat fractional electrons in AIMD has been based on the idea of a grand canonical density functional approach using separate potential energy surfaces for the two distinct oxidation states of the same chemical system.11 This computationally appealing technique allows one to fix the electronic chemical potential, thus allowing the system to change the number of electrons without introducing fractional electrons explicitly, while not allowing

1. MOTIVATION Computer simulations should be carried out in well-defined statistical−mechanical ensembles.1,2 In the realm of pencil and paper calculations, the corresponding partition functions or thermodynamic potentials (free energies) can be easily interconverted using Laplace and Legendre transforms, respectively, that exchange the desired conjugate variables. This is in stark contrast to the numerical sampling of statistical ensembles using molecular simulations. The “natural ensemble” of Monte Carlo (MC) simulations is the canonical or (N, V, T) ensemble where particle number, volume, and temperature are kept fixed during sampling by default, whereas the microcanonical or (N, V, E) ensemble is straightforwardly generated by standard molecular dynamics (MD). However, in order to control temperature instead of the total energy in MD, for instance, ingenious methods such as extended Lagrangian techniques1,2 are required that suitably modify the plain Newtonian dynamics to establish the canonical ensemble computationally. The situation gets much more involved for MD upon replacing the constant particle number N, being an integer, by a constant chemical potential μ, which is the continuous conjugate variable according to the grand canonical or (μ, V, T) ensemble. Indeed, such grand canonical formulations of MD have been devised long ago3−7 by following the extended system dynamics approach to introduce a continuous particle number variation variable (and thereby something like fractional particles) that couples the simulated system to a chemical potential reservoir. It is mentioned in passing that a plethora of powerful methods has been introduced over the years to compute chemical potentials from computer © XXXX American Chemical Society

Received: March 13, 2017 Published: June 21, 2017 A

DOI: 10.1021/acs.jctc.7b00263 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

usual ab initio thermodynamic integration (AITI) approach in the latter case, as well as its zero-temperature limit which yields the well-established ab initio thermodynamics (AITD) method. The implementation as a novel computational method in the AIMD program package CPMD23 is addressed in section 4. In the results section, we start by validating our bicanonical AIMD approach with reference to its canonical ensemble limits in section 5.1 before applying it to compute solvation free energies for Na+ and K+ ions in aqueous solution in section 5.2. We close the presentation in section 6 where we also present some ideas for future research directions.

the number of nuclei to adjust according to a constant chemical potential of a full chemical species. This particular grand canonical scheme11 has been used very successfully to study redox reactions involving atoms and molecules in aqueous solution, thus adding the required number of electrons to reach the desired reduced state or vice versa (without simultaneously changing the number of nuclei and thus atoms in the system) as applied in refs 12 and 13 and in a wealth of subsequent work that is reviewed in ref 14. Alternatively, redox free energies in solution have also been obtained by devising an AIMD-based thermodynamic integration (AITI) approach that is explicitly based on propagating the two potential energy surfaces associated with the oxidized and reduced states of the system,14,15 again without changing the number of nuclei during the reaction. In order to access de/protonation reactions, being reactions where the number of protons changes (by ±nH+ in general) whereas the number of electrons remains constant in the system, another sophisticated thermodynamic integration technique has been developed within the AIMD framework, which is able to annihilate (or create) a proton in solution.16 It transmutes a protonated acid in aqueous solution, AH(aq), via a continuous switching parameter η = 0 → 1 into its conjugate base, A−(aq), such that the associated free energy change, and thus the pKa value of AH(aq), is obtained from AITI starting in the protonated initial η = 0 and ending in the deprotonated final η = 1 state of the solution. This powerful approach provides access to de/protonation free energies14,16 and thus to pKa/pKb values of aqueous solutes including the solvation free energy of the proton itself in water in the limit of its Eigen cation, i.e., the H−bonded H3O+(aq) complex being AH(aq). Finally, in such cases where both the number of protons and electrons must be changed, the full reaction is split into the separate redox and de/protonation elementary reactions in the spirit of a thermodynamic cycle applied to solvation (using Born−Haber cycle and Hess’s law concepts).17 The electrostatic contribution to the hydration free energy of ions in water has also been obtained from a particular thermodynamic integration technique that considers the ionic charge as the order parameter,18 which is analogous to what is customarily done in force field simulations to compute this free energy term19 (where typically a Lennard−Jones interaction center that is already embedded in the aqueous solution is increasingly charged until the full ionic charge is reached20). Last but not least, thermodynamic integration methods formulated in terms of “semiconductor defect language” in the realm of electrochemistry21 have been used more recently to show that dispersion (“van der Waals”) interactions on redox levels are small since the solvation properties of the ionic species are not much affected by these nonlocal interactions.22 Here, we introduce what we call bicanonical AIMD simulations where the number of chemical species, consisting of a set of predefined nuclei (being not necessarily protons) and a fixed number of electrons, can be changed during the simulation in the sense of an open system that is kept at constant chemical potential. Having explained our idea in words in section 2, we will derive in section 3 what we call the effective bicanonical potential and the resulting propagation scheme, being at the heart of our bicanonical AIMD simulation technique. This presentation includes the analysis of its behavior if either the chemical potential or the average fractional particle number corresponding to a predefined chemical species is kept constant, its direct connection to the

2. CONCEPT In the present work, having the aim to change both the number of nuclei and electrons at the same time as required when the chemical potential of a chemical species is kept constant, we go back to interesting ideas of Swope and Andersen who devised an original method for the calculation of chemical potentials of liquids and solids.24 To this end, they invented a new equilibrium statistical thermodynamic ensemble, called the “bicanonical ensemble”.24 It can be considered, on the one hand, to be essentially a grand canonical ensemble that is rigorously restricted to two states that differ in particle number at a fixed chemical potential (which they then sampled conventionally via insertion and deletion of particles with the aim to compute the chemical potential by explicitly propagating only one system by MC or MD). On the other hand, the bicanonical ensemble can be viewed in terms of two coupled canonical ensembles from which the chemical potential can be rigorously obtained via the derivative of the Helmholtz free energy with respect to particle number. Based on the rigorous bicanonical formalism, we demonstrate herein that also the standard ab initio thermodynamic integration approach can be recovered when continuously changing the coupling parameter of the two canonical ensembles, which represent the two limiting states in the sense of thermodynamic integration. We consider this concept to be of great interest for computer simulation since this peculiar ensemble is ideally suited to properly link two hitherto disconnected worlds at the level of practical computer simulations, namely AIMD simulations where either the chemical potential or the fractional particle number of a certain chemical species is kept constant. Transferring the bicanonical ensemble concept24 to the world of AIMD, the two systems that are coupled within the bicanonical ensemble differ in terms of the number of both nuclei and electrons contained in two separate simulation supercells (which have the same volume V and are kept at same temperature T) that are used for performing the underlying electronic structure calculations. At this point, an interesting analogy to what is called ab initio thermodynamics might come to mind (see e.g., the review in ref 25, which also contains the references to the original work on defect formation free energies in solids followed by the calculation of surface phase diagrams). In AITD, two systems that contain a different number of nuclei and electrons in their supercells are referenced to a constant chemical potential of those species by which the systems differ. One such system could be a solid that hosts an interstitial defect or a surface on which an admolecule is chemisorbed, whereas the other system, being the reference, does not contain these additional species. The presence of the interstitial or admolecule is governed by an associated chemical potential: being sufficiently high it is thermodynamically preferable for the additional species to be B


Article


thus avoids both dealing with fractional nuclei and/or electrons and explicit insertion or deletion of particles.

present in the solid or on the surface, whereas in the limit of low chemical potential these species prefer to be in the “reservoir” (which is not explicitly represented in AITD calculations). Thus, starting with a perfect solid or a clean surface, the additional species eventually enter the scene upon increasing the respective chemical potential, which depends on temperature and pressure. Replacing the free energies of the two systems by the total (electronic structure) energies in practical AITD calculations,26 thus neglecting entropy and temperature effects altogether, the relative stability of the two systems and thus defect or adsorption (free) energies can be obtained as a function of the thermodynamic state point using conventional static electronic structure calculations (see e.g., ref 25 and literature cited therein). Importantly, it is not necessary in AITD calculations to perform particle insertion or deletion manipulations since the two systems are treated independently from each other. In our bicanonical AIMD approach, we take the original bicanonical idea24 literal in the sense that we couple two independent systems individually described by density functional theory (but differing in the number of nuclei and electrons, both being integer numbers and not fractional) by defining what we call an “effective bicanonical potential”, being their weighted superposition which explicitly depends on the (arbitrary but constant) chemical potential. The nuclear gradient of this potential yields the “bicanonical force” (being a superposition of the forces stemming from the two independent systems that are weighted according to the specified chemical potential) that acts on the nuclei in both systems and, thus, is used to propagate the two underlying systems simultaneously in the sense of AIMD simulations. By construction, bicanonical AIMD offers the significant advantage of implicitly representing fractional particle numbers by a superposition of two canonical ensembles that exclusively contain the usual integer numbers of nuclei and electrons, instead of explicitly introducing fractional particles. Thus, in bicanonical AIMD, fractions of particles emerge at the level of the thermal averages as a result of sampling, rather than being present already at the level of each and every individual configuration that is sampled. It is in this sense that we employ the term “fractional particle number” in the context of bicanonical AIMD simulations. The advantage comes from the well-known fact that both using explicitly fractional nuclei and electrons in electronic structure calculations and any involvement of MC sampling procedures are very difficult to efficiently implement in practical AIMD schemes.8 This is in contrast to force-field-based MD and MC simulations where the use of a fraction of a particle is straightforwardly implemented by a simple scaling of force field terms with appropriate coupling parameters (see, e.g., ref 27). Our proof of principle application of bicanonical AIMD is the calculation of the Helmholtz solvation free energy of sodium and potassium ions in water at ambient conditions, Na+(aq) and K+(aq). It is shown that the bicanonical approach is able to continuously introduce a single such metal ion into pure bulk water upon increasing the chemical potential, thus slowly building up the typical radial distribution function of a solvated cation in water. We show that the solvation free energy can be computed via both the constant chemical potential and fractional particle number approaches, providing consistent numbers in accord with the literature. This is the proof of concept that bicanonical AIMD indeed works in practice and

3. THEORY Before we derive the basic relations of our bicanonical (AI)MD technique in section 3.3, we quickly review the basic definitions of the canonical and grand canonical ensembles in statistical mechanics in section 3.1 merely to establish our notation before we introduce the bicanonical ensemble concept of Swope and Andersen24 in section 3.2. Readers who are familiar with the bicanonical ensemble may thus jump immediately to section 3.3. All derivations will be carried out for multicomponent systems at constant volume V and constant temperature T and for particle exchange of the first chemical component, s = 1; generalization to more components follows from the presented formalism. 3.1. Multicomponent Thermodynamic Ensembles. We consider a multicomponent system that contains a mixture of r different particle components s = 1, ..., r (including the associated number of electrons in the context of an electronic structure description as required in the framework of AIMD). The number of particles and the chemical potential of each species s are given by Ns and μs , respectively. The grand canonical ensemble of this open multicomponent system is characterized by the parameters (μ1, ..., μr , V, T). The total number of particles is N = ∑rs=1 Ns , and qN and pN denote the coordinates and momenta of all particles in the system. Next, H(qN, pN) = K(pN) + U(qN) is the classical Hamiltonian for the N

p2

N particles with the kinetic energy K (p N ) = ∑i = 1 2mi and the i

potential energy U(qN); note that U(qN) can (and will be later on) obtained from electronic structure calculations given the configuration qN. The probability that a system chosen at random from the grand canonical ensemble has a total number of N particles and is at phase point (qN, pN) in the phase space volume element dqN dpN is given by P(qN, pN) dqN dpN with

P(q N , p N ) =

⎛ r ⎜⎜∏ ⎝ s=1

⎞

N N 1 λ Ns⎟e−βH(q , p ) h3NsNs ! s ⎟

⎠

Ξ(μ1 , ..., μr , V , T )

(1)

where λs = eβμs is the activity of the s-th component, β = 1/kBT is the reciprocal temperature, and Ξ is the grand canonical partition function Ξ(μ1 , ..., μr , V , T ) ∞

=

∞

∑ ∑ N1= 0 N2 = 0

r

∞

...

∑

Q (N1 , ..., Nr , V , T ) ∏ λsNs s=1

Nr = 0

(2)

The function Q in eq 2 is the multicomponent partition function of the canonical ensemble, characterized by the parameters (N1, ..., Nr , V, T), Q (N1 , ..., Nr , V , T ) ⎛ r 1 ⎞ = ⎜⎜∏ 3N ⎟⎟ ⎝ s = 1 h sNs! ⎠

∞

N

N

∫−∞ dp N ∫V dq Ne−βH(q ,p )

(3)

Integration over the momenta reduces the canonical partition function to C


Article

Journal of Chemical Theory and Computation ⎛ r 1 ⎞ Q (N1 , ..., Nr , V , T ) = ⎜⎜∏ 3N ⎟⎟ ⎝ s = 1 Λs sNs! ⎠

Q I = Q (N1I , N2 , ..., Nr , V , T )

N

∫V dq Ne−βU(q )

Q II = Q (N1II , N2 , ..., Nr , V , T ) (4)

The thermodynamic variables of this multicomponent bicanonical ensemble are (μ†1, N2, ..., Nr , V, T). The dagger has been introduced, following the original notation of Swope and Andersen,24 in order to clearly indicate that the bicanonical chemical potential μ†1 is related to but not equal to the proper thermodynamic chemical potential μ1 of the unconstrained grand canonical ensemble. In the bicanonical ensemble the number of particles of the unconstrained component s = 1 can only change between the two fixed integer values NI1 and NII1 . In the following we will always assume that II is the larger and I is the smaller system, i.e. NII1 − NI1 = ΔN1 > 0 without loss of generality. 3.2.1. Bicanonical Ensemble as Restricted Grand Canonical Ensemble. On the one hand, the bicanonical ensemble can be regarded as an approximation to a true grand canonical ensemble, in which only two systems I and II are important and all others can be neglected. Accordingly, following eq 5, we can introduce the multicomponent bicanonical thermodynamic potential

with the thermal wavelength Λs = h/(2πmskBT)1/2. From the partition functions, associated thermodynamic potentials can be derived. From the grand canonical partition function Ξ we obtain the multicomponent grand canonical potential Ω(μ1 , ..., μr , V , T ) = −kBT ln Ξ(μ1 , ..., μr , V , T ) = −p(μ1 , ..., μr , V , T )V

(5)

which is equal to minus the pressure p(μ1, ..., μr , V, T) times the volume V of the system. The average number of particles N̅ s of any species s in the grand canonical ensemble can be evaluated by taking the derivative of the grand canonical potential Ω with respect to the corresponding chemical potential μs ⎛ ∂Ω ⎞ ⎟⎟ Ns̅ (μ1 , ..., μr , V , T ) = −⎜⎜ ⎝ ∂μs ⎠ μ ,..., μ , V , T 1

r

(6)

Ω Γ(μ1† , N2 , ..., Nr , V , T )

Correspondingly, the Helmholtz free energy A(N1 , ..., Nr , V , T ) = − kBT ln Q (N1 , ..., Nr , V , T )

= − kBT ln Γ(μ1† , N2 , ..., Nr , V , T )

1

r

Therefore, the average number of particles of species s = 1 in the bicanonical ensemble can be calculated, following eq 6, as N1̅ †

Q II

(8)

Q

Q (N1 , N2 , ..., Nr , V , T )λ1N1 (9)

Next, we remove all but two canonical systems I and II with particle numbers NI1 and NII1 of species s = 1. The total particle numbers for the two systems are NI = NI1 + ∑rs=2Ns and NII = NII1 + ∑rs=2 Ns . After imposing this restriction, eq 9 becomes the partition function of the multicomponent bicanonical ensemble Γ(μ1† ,

N2 , ..., Nr , V , T ) =

I Q Iλ1N1

+

II Q IIλ1N1

(14)

If the average number of particles in the bicanonical ensemble is close to that of system II, the value of the average bicanonical fraction will become x†1 → 1, whereas in the case that the average number of particles is close to the one of system I, we obtain x†1 → 0. 3.2.2. Bicanonical Ensemble as Superposition of Two Canonical Ensembles. On the other hand, the bicanonical ensemble can be regarded as the superposition of two canonical ensembles with particle numbers NI and NII. Swope and Anderson24 used this point of view to calculate the proper thermodynamic chemical potential μ1 of component s = 1 by using the Helmholtz free energies AI = A(NI1, N2, ..., Nr , V, T) of system I and AII = A(NII1 , N2, ..., Nr , V, T) of system II and by approximating the derivative in eq 8 by the finite difference

∞

∑

(13)

†

e βμ1 ΔN1 N1̅ † − N1I QI † x1 ≔ = † Q II ΔN1 1 + I e βμ1 ΔN1

Ξ̃(μ1 , N2 , ..., Nr , V , T ) N1= 0

I II ⎛ ∂Ω ⎞ N1Iλ1N1 Q I + N1IIλ1N1 Q II Γ = = −⎜⎜ † ⎟⎟ I II λ1N1 Q I + λ1N1 Q II ⎝ ∂μ1 ⎠V , T

Taking system I as reference we define the average bicanonical fraction of particles

3.2. Multicomponent Bicanonical Ensemble. We now constrain the multicomponent grand canonical ensemble in such a way that only one component s = 1 can exchange particles with a reservoir with associated chemical potential μ1, whereas all other particle numbers are kept fixed at N2, ..., Nr . Thus, the resulting ensemble is only grand canonical with respect to the component s = 1, whereas it is canonical concerning all other components s = 2, ..., r. The grand canonical partition function of the multicomponent system, eq 2, then reduces to the following expression

=

(12)

N†1

(7)

is the associated thermodynamic potential of the multicomponent canonical ensemble. The chemical potential μs of any particle component s is obtained from the Helmholtz free energy A by taking the derivative with respect to the particle number Ns ⎛ ∂A ⎞ μs (N1 , ..., Nr , V , T ) = ⎜ ⎟ ⎝ ∂Ns ⎠ N ,..., N , V , T

(11)

μ1(N1II , N2 , ..., Nr , V , T ) =

k T Q II AII − AI ΔA = = − B ln I ΔN1 ΔN1 ΔN1 Q

(15)

Thus, ΔA = A − A is the change in the Helmholtz free energy for removing ΔN1 particles from system II. By introducing the so-called bicanonical ratio II

(10)

with D

I


Article

Journal of Chemical Theory and Computation R(μ1† , N2 , ..., Nr , V , T ) =

x1†

=

1 − x1†

Q II βμ1†ΔN1 e QI

and I

(16)

U I = U (q N ) = U (q1I, ..., q IN I)

which is the ratio of the probability that a system chosen at random from our bicanonical ensemble has NII particles to the probability that it has NI particles,24 and by using eq 14, the thermodynamic chemical potential μ1 can be related to the bicanonical chemical potential μ†1 via μ1 = μ1† −

kBT ln R(μ1† , N2 , ..., Nr , V , T ) ΔN1

II

U II = U (q N ) = U (q1II, ..., q IIN II)

are the potential energy functions of the two canonical systems. The factor γ1, eq 22, cancels the contributions arising from the ΔN1 fictitious ideal gas particles, which have been formally introduced via eq 20 merely as a mathematical trick to extend II

the configuration integral ∫ dq N in the bicanonical partition V function to account for the higher-dimensional configuration space of the larger system II. 3.3.2. Effective Bicanonical Potential and Forces. After having rewritten the bicanonical partition function eq 10 as a single configuration space integral, eq 21, we can introduce the II effective bicanonical potential UΓ(qN )

(17)

Thus, the true thermodynamic chemical potential is not simply μ†1, but it has to be calculated from the relative probabilities of the two particle numbers in the bicanonical ensemble.24 3.3. Derivation of the Effective Bicanonical Potential. In the following we will derive an effective bicanonical potential that will allow us to sample the multicomponent bicanonical ensemble via AIMD. The effective bicanonical potential will guarantee a smooth transition between the two potential energy surfaces of the two canonical systems I and II during the AIMD sampling. 3.3.1. Thermodynamic Extension of the Small System. To this end we extend the partition function QI of the small canonical system I of eqs 4 and 11 by ΔN1 = NII1 − NI1 additional coordinates. We introduce the canonical partition function Qid(ΔN1, V, T) of an ideal gas of ΔN1 particles of species s = 1 in volume V at temperature T, i.e. Q id(ΔN1 , V , T ) =

h

ΔN1

1

= (h / Λ1)3ΔN1

V

II

= V ΔN1

(19)

Using this idea, we then extend the bicanonical partition function eq 10 by introducing what could be called a “thermodynamic 1” Γ(μ1† , N2 , ..., Nr , V , T ) =

Q id Q

id

I

II

Q Iλ1N1 + Q IIλ1N1

(20)

For the nominator and denominator in this thermodynamic extension, we insert eqs 18 and 19, respectively, and the canonical partition functions QI and QII of system I and II will be expressed by their corresponding configuration integrals, eq I

f αΓ = −

II

4. By combining the integrals ∫ dq ΔN1·∫ dq N ... = ∫ dq N ..., V V V the bicanonical partition function eq 10 can be written as a single integral over configuration space of the larger system II

×

I

∫V dq N (γ1e−βU e

βμ1†N1I

II

+ e − βU e

βμ1†N1II

)

f αI = −

(21)

Λ13ΔN1 N1II! V ΔN1 N1I!

(23)

(

I

†

I

II

†

II

)

(24)

II ∂ UΓ(q N ) = (1 − x1)f αI + x1f αII ∂q α

(25)

∂ I NI U (q ) ∂q α

and

f αII = −

∂ II N II U (q ) ∂q α

Thus, the effective bicanonical forces f αΓ are obtained by a simple linear superposition of the forces f αI and f αII of the two canonical systems I and II. The associated weighting factor for the bicanonical averaging of the forces in eq 25 turns out to be II configuration-dependent, x1 = x1(qN ) and is given by

with the factor γ1 =

)

with

⎛ r 1 1 ⎞ ⎜⎜∏ 3N ⎟⎟ Γ(μ1† , N2 , ..., Nr , V , T ) = II Λ13N1 N1II! ⎝ s = 2 Λs sNs! ⎠ II

N II

Γ

In the bicanonical ensemble the effective bicanonical potential UΓ is the analogue to the potential energy surface U in the canonical ensemble; see eq 4. Thus, instead of actually removing and adding particles as in the original scheme of Swope and Andersen,24 the bicanonical ensemble can also be sampled by propagating simultaneously the two canonical systems I and II on the energy surface of the effective bicanonical potential UΓ, which depends on both, constant temperature T and constant bicanonical chemical potential μ†1. In this case, the coordinates of the particles in the small and the large system are no longer treated independently, but the effective bicanonical potential UΓ is always evaluated for the same positions of all particles that are common to both systems. The forces f αΓ acting on the individual atoms α in the coupled bicanonical system are obtained by differentiation of the effective bicanonical potential UΓ with respect to the nuclear coordinates (qα) at constant T and μ†1. After some algebra the following simple expression is obtained

ΔN1

Λ13ΔN1ΔN1!

II

∫V dq N e−βU (q

UΓ(q N ) ≔ kBT ln γ1e−βU e βμ1 N1 + e−βU e βμ1 N1

(18)

1

   

⎛ r 1 ⎞ ⎜⎜∏ 3N ⎟⎟ II Λ13N1 N1II! ⎝ s = 2 Λs sNs! ⎠ 1

with

1 ΔN1!

∫−∞ dpΔN e−βK(p ) ∫V dq ΔN

=

=

3ΔN1

∞

×

Γ(μ1† , N2 , ..., Nr , V , T )

x1 = (1 + γ1exp[β(ΔU − ΔN1μ1†)])−1

(22) E

(26)


Article


for example, if a chemical reaction is taking place, in which one species is consumed and has to be replenished from a reservoir. Thus, our bicanonical AIMD approach reduces to the wellestablished AITD method (see, e.g., ref 25 and literature cited therein) in the suitable limits and, therefore, can be considered to be a continuous extension of AITD that is accessible to molecular dynamics sampling. Practical AITD calculations, however, are typically carried out in the zero temperature limit,25,26 thus replacing in the above-given formulas the Helmholtz free energies at finite temperature, AI and AII, by the corresponding total (electronic) energies, UI and UII. The latter are easily obtained from static single-point electronic structure calculations that are performed for a finite set of preselected I II particle configurations {qN } and {qN } in order to find the preferred configuration at the thermodynamic state point of interest. This procedure neglects altogether thermal activation and entropic effects in such standard AITD calculations, whereas they are included in bicanonical AIMD simulations carried out at finite temperatures. 3.3.4. Sampling at Constant Chemical Potential. The first option for sampling the effective bicanonical potential UΓ is by running AIMD simulations at constant chemical potential, μ†1, where both systems, I and II, are propagated simultaneously. The instantaneous weight x1(t) of the two configurations along the trajectory, which can be regarded as a fractional particle number of species s = 1, is calculated from the energy gap II ΔU(qN ) at time t via eq 26, and the forces acting on the individual atoms at time t are given by the weighted average eq 25 with x1(t) as the force weight factor. Recall that what we call f ractional particle numbers in bicanonical AIMD, which refers to nuclei and electrons, emerges by superimposing and thus averaging two canonical systems that both contain integer numbers of nuclei and electrons. The average bicanonical fraction of particles x†1 can then be determined from the bicanonical AIMD trajectory by taking the ensemble (i.e., time) average of the fractional particle number x1

It is determined by the difference of the configurationdependent potential energies of the two canonical systems, II II I ΔU(qN ) = UII(qN ) − UI(qN ), and it represents the instantaneous fractions of particles in a sampling scheme where the bicanonical chemical potential, μ†1, is kept fixed. If system II is energetically much more favorable than system I, II i.e., ΔU < ΔN1μ†1 for a given configuration qN , the system weight will be mostly on II and the weight factor becomes x1 → 1, whereas in the case of ΔU > ΔN1μ†1, system I is energetically more favorable and x1 → 0. It is important, however, to note that the configurationdependent force-weighting factor x1 of eq 26 differs conceptionally from the average bicanonical fraction of particles, x†1, as thermodynamically defined in the bicanonical ensemble by eq 14, because the former relates the potential energies UI and UII of the two systems rather than the free energies AI and AII as in the latter case. This fundamental difference becomes obvious by combining eq 14 with eq 15 to yield x1† = (1 + exp[β(ΔA − ΔN1 μ1†)])−1

(27)

being the thermodynamic connection between the bicanonical fraction of particles of species s = 1, x†1, and the associated bicanonical chemical potential, μ†1, via the corresponding Helmholtz free energy difference, ΔA. 3.3.3. Limiting Cases and AITD. By using the configurationdependent weight factor x1, the effective bicanonical potential of eq 24 can be rewritten in the two following, equivalent ways, ⎧U I − N Iμ † − k T ln γ + k T ln(1 − x ) ⎪ II 1 1 B B 1 1 UΓ(q N ) = ⎨ ⎪U II − N IIμ † + k T ln x ⎩ 1 1 B 1 (28)

We now take a closer look at two limiting cases for the effective bicanonical potential of eq 28 ⎧U I − N Iμ † − k T ln γ for x = 0 ⎪ 1 1 B 1 1 UΓ(q ) = ⎨ ⎪U II − N IIμ † for x1 = 1 ⎩ 1 1 N II

x1† = ⟨x1(t )⟩μ1†

at constant By combining eqs 15−17 (or, equivalently, by inverting eq 27), finally the difference in the Helmholtz free energy ΔA between system I and II can be calculated from

Similarly, the bicanonical partition function becomes Γ(μ1† ,

(29)

μ†1.

⎧ I βμ1†N1I for x1 = 0 ⎪Q e N2 , ..., Nr , V , T ) = ⎨ † II ⎪Q IIe βμ1 N1 for x = 1 ⎩ 1

⎛ x† ⎞ ΔA = ΔN1 μ1† − kBT ln⎜⎜ 1 † ⎟⎟ ⎝ 1 − x1 ⎠

and the thermodynamic potential ΩΓ = −kBT ln Γ is

(30)

Thus, ΔA is directly obtained from the relative probabilities of the two particle numbers in the bicanonical ensemble.24 As already pointed out by Swope and Andersen,24 the right-hand side of eq 30 is independent of μ†1 for fixed values of (μ†1, N2, ..., Nr , V, T), although each of the two terms depend on it. Therefore, in principle, any choice of μ†1 will give the same result for ΔA; however, for an efficient sampling of the bicanonical ensemble μ†1 should be about ΔA/ΔN1 so that x†1 becomes close to 1/2. Since the appropriate choice of μ†1 is not known a priori, usually several bicanonical AIMD simulations for different values of μ†1 are performed and ΔA is determined11 from the graph of x†1(μ†1) at the value x†1 = 0.5. 3.3.5. Sampling at Constant Fractional Particle Number. As an alternative to running the bicanonical AIMD simulations at constant chemical potential μ 1† and evaluating the instantaneous weights x1(t) for the force calculation using eq

⎧ AI − N Iμ † for x = 0 ⎪ 1 1 1 † Ω Γ(μ1 , N2 , ..., Nr , V , T ) = ⎨ ⎪ AII − N IIμ † for x = 1 ⎩ 1 1 1

This is the known expression from standard AITD, which is used to deduce whether system I or II is thermodynamically more favorable at a given chemical potential μ†1. If x1 is close to zero, system I has the highest weight and the bicanonical ensemble is dominated by AI. In contrast, for x1 close to one, the bicanonical system is described by AII. In our extension of the thermodynamic formalism, the continuous weight function x1 now guarantees for a smooth transition between the two energy landscapes associated with systems I and II, if the relative thermodynamic stability of the two systems changes in the course of the bicanonical AIMD simulation. This can occur, F


Article

Journal of Chemical Theory and Computation 26, we can also keep the fractional weights x1 constant and determine the corresponding instantaneous bicanonical chemical potential μ†1(t) by inverting eq 26, which thereby becomes configuration-dependent

common to the two canonical systems. This implies in the practical implementation that at any bicanonical AIMD time step, all positions and momenta of the nuclei of system I can always be created from those of system II by removing the additional nuclei α = {NI1 + 1, ..., NII1 }, i.e., the M+ ion in the present case. Sampling the bicanonical potential at the AIMD level requires the propagation of the full electronic structure of both underlying canonical systems at each AIMD step. In order to guarantee stable sampling of the bicanonical ensemble at the AIMD level in terms of a superposition of two canonical systems, we must take special care for those situations in which the bicanonical force weight function x1 → 0 or becomes equal to zero. In this limit, the x1 weight is mostly or entirely on the smaller canonical system I containing only nuclei α = {1, ..., NI}. In this limit, system II does not contribute at all to f αΓ and the generated AIMD dynamics is effectively that of the small canonical system I. This implies for the AIMD propagation of system II that the additional nuclei α = {NI1 + 1, ..., NII1 }, which are only in the larger system, namely the M+ cation in the present case, are now fully decoupled from the rest of that system at the level of the weighted forces which act on them. This is the expected scenario, as α = {NI1 + 1, ..., NII1 } propagate according to free particles dynamics like they would do in a reservoir. However, technically these “free” particles still reside in the simulation box of system II, but given that they essentially feel no forces if x1 → 0, very close contacts or even overlapping positions of these additional particles with the other particles in that system, i.e., α = {1, ..., NI}, may be created. Clearly, such events must be avoided because they lead to instabilities in the AIMD sampling as the convergence of the electronic structure will become very time-consuming or will even fail completely. In order to prevent such artifacts, we apply a very short-ranged repulsive potential acting exclusively on the additional nuclei α = {NI1 + 1, ..., NII1 } that are only hosted by system II

⎛ x ⎞ ΔN1 μ1† (t ) = ΔU (t ) + kBT ln γ1 + kBT ln⎜ 1 ⎟ ⎝ 1 − x1 ⎠ (31)

By doing so, however, we no longer sample the proper effective bicanonical potential UΓ, but simply a specific fixed linear superposition (1 − x1)UI + x1UII of the potential energy surfaces of system I and II depending on a preselected x1 value between zero and unity. Nevertheless, we can calculate the instantaneous μ†1(t) along the AIMD trajectory from eq 31, and thus, we can compute its ensemble average ⟨μ†1(t)⟩x1 at the fixed value of the weight factor x1. By doing so, we recover from bicanonical AIMD the well-established thermodynamic integration scheme and thus AITI.16 The difference in the Helmholtz free energy ΔA between system I and II is, therefore, given by the integral ΔA =

∫0

1

⟨ΔU (t )⟩x1 dx1 1

=

∫0 (ΔN1⟨μ1†(t )⟩x

1

− kBT ln γ1) dx1

(32)

along the constraint coordinate x1 that transmutes the system from I to II upon increasing it from zero to unity. Here, ⟨ΔU(t)⟩x1 and thus the right-hand side of eq 32 is obtained 1

from eq 31; note that the full ∫ dx1 integral over the rightmost 0 term of eq 31 yields rigorously zero in eq 32 and, thus, is not present in practical AITI calculations being based on the latter equation.

4. COMPUTATIONAL DETAILS In order to validate the bicanonical AIMD approach, we have set up the two limiting canonical systems I and II consisting of 32H2O and M+32H2O where M+ is a single metal ion in order to represent the pure solvent and the solution, respectively. As usual, periodic boundary conditions were applied to a cubic supercell of linear dimension L = 9.8280 Å. Next, we have implemented the bicanonical ensemble as a novel method in the Bochum in-house branch of version 4.1 of the CPMD program suite.23 The implementation allows one to perform both types of AIMD sampling that are offered by the bicanonical ensemble via UΓ (see section 3). In the first mode, the bicanonical force weight x1 is set to a constant for the entire sampling run, which can be fractional and also numerically assume its two limits of x1 = 0 and unity, in order to compute bicanonical ensemble averages of the associated chemical potential ⟨μ†1⟩x1. In the second mode, we can set the bicanonical chemical potential μ†1 of species s = 1, which herein is a single metal ion M+, to a constant and calculate the instantaneous x1(t) value at each AIMD step following eq 26. In the two limiting cases of x1 = 0 and x1 = 1, the forces on nuclei α = {1, ..., NI} and α = {1, ..., NII} correspond to the 32H2O and M+32H2O systems, respectively, fully contribute to f αΓ and, as such, the bicanonical dynamics is effectively that of the corresponding canonical system I and II, respectively. Note that in our bicanonical AIMD sampling approach the positions of the nuclei are identical for the subset of atoms which is

Urep

⎧ NII NI ⎪ ∑ ∑C for |q α − q β| ≤ C3 1 ⎪ I ⎪ α=N +1 β=1 2 =⎨ ⎛ −C (|q − q |− C ) ⎞ ⎪ × ⎝⎜1 − e( 2 α β 3 )⎠⎟ ⎪ ⎪0 for |q α − q β| > C3 ⎩ (33) −1

where C1 = 1 mHa, C2 = 1 Å , and C3 = 0.1 Å. This shielding potential Urep adds additional force contributions to the forces f αII of nuclei α = {NI1 + 1, ..., NII1 }, which in our case is a single metal ion M+. As we will show later in our structural analysis in section 5.1, the excluded volume effect created hereby is negligible. The AIMD simulations were carried out within the Car− Parrinello28 propagation scheme.8 We used density functional theory with the PBE exchange correlation functional,29,30 a plane wave basis set up to a kinetic energy cutoff of 25 Ry, and ultrasoft pseudopotentials31 to describe the core electrons; the Na+ and K+ are represented using small cores, i.e., [1s2] and [1s2, 2s2, 2p6], respectively. A fictitious mass parameter of 900 a.u. for the electrons was employed to decouple their fictitious dynamics adiabatically8 from the real time dynamics of the nuclei. The mass of deuterium was used for hydrogen ions to allow for an AIMD time step of Δt = 4 a.u. (≈0.097 fs). The electronic structure within these AIMD simulations was G


Article


highlighted. Exploring again the limiting cases first, in Figure 2 we directly compare at the level of the Na+−O radial

calculated using the Γ-point of the Brillouin zone. Temperature control was achieved by Nosé−Hoover chains8,32 to thermostat both the nuclear and the orbital degrees of freedom. In order to most efficiently generate the canonical ensemble, in the case of nuclei, one independent Nosé−Hoover chain has been coupled to each of the Cartesian coordinates of all mobile ions (called “massive thermostatting”). As usual, the target temperature was set to 400 K in order to approximately generate the properties of liquid water and aqueous solutions at room temperature when using the PBE functional, which is known to overstructure these systems.33

5. APPLICATION: SOLVATION FREE ENERGIES 5.1. Structural Properties: Limiting Canonical Systems and Fractional Bicanonical Regime. As a first validation step, we demonstrate that our bicanonical AIMD implementation can indeed rigorously reproduce the two canonical limits, which yield the large and the small system for x1 → 1 and x1 → 0, respectively, where the runs are performed while keeping the (fractional) force weight x1 constant during propagation. To this end, we investigated the Na+(aq) solution as a function of x1, being a superposition of the 32H2O and Na+32H2O systems. With this setup, two AIMD trajectories of about 0.5 ns length were generated where the bicanonical ensemble was forced to sample these two pure canonical limits by setting x1 = 0 and x1 = 1, respectively. Analysis of the resulting bicanonical trajectories in direct comparison to the ones obtained from conventional canonical AIMD propagation of the respective 32H2O and Na+32H2O systems using the same computational setup showed that numerically identical results were produced. This demonstrates that our implementation of bicanonical AIMD in the CPMD package23 is able to correctly reproduce the exact canonical limits. In the next step, we went on to analyze the structural properties which have been generated by sampling aqueous Na+(aq) solutions that contain a fractional Na+ ion via the bicanonical ensemble at constant but fractional x1; also these trajectories have been propagated for about 0.5 ns each. A representative snapshot sampled from the bicanonical AIMD simulation of this system at x1 = 0.5 is depicted in Figure 1, wherein the water molecules of the first solvation shell are

Figure 2. Radial distribution functions of the Na+−O distances obtained from a canonical AIMD simulation of the Na+(aq) solution using the Na+32H2O system (red line with circles) and from bicanonical AIMD simulations using constant weights of x1 = 1 (black line) and x1 = 0 (blue line) corresponding to the Na+(aq) solution and to the pure solvent, respectively.

distribution functions (RDFs) the canonical Na+32H2O system to our bicanonical system with force weights of x1 = 1 and x1 = 0, thus setting full weight on the solvated Na+32H2O system and the pure 32H2O system, respectively. In the former case, besides very small deviations due to sampling statistics, these functions match essentially perfectly which demonstrates that our bicanonical AIMD method is able to correctly describe this canonical limit at the level of the solvation shell structure. The more challenging case is the opposite limit with a force weight of x1 = 0 where the forces according to the effective bicanonical potential come exclusively from the pure liquid 32H2O system. Yet, we can compute the RDF involving the Na+ ion in the large system II, which however does not contribute to the bicanonical force in this limit! Thus, the motion of the Na+ ion is expected to be completely decoupled from that of the water molecules, and thus, this solute species is expected to move like a noninteracting (i.e., free) particle in a pure solvent in the x1 = 0 limit. Indeed, this behavior is numerically recovered at the level of the corresponding RDF since g(rNa+−O) ≈ 1 for all but very short rNa+−O distances according to Figure 2. The pronounced fluctuations in the short-range limit are due to the sampling statistics, since this RDF is very hard to converge in the increasingly smaller radial volume elements as r → 0 even within the accuracy provided by the rather long bicanonical AIMD trajectories of about 0.5 ns; note that the distance regime on the order of 0.1 Å is never accessed in any molecular liquid that is subject to physical interactions. Recall from section 4, however, that there still is a very short-ranged and steeply repulsive potential, eq 33, that acts between the Na+ core and all O and H sites of all water molecules in order to prevent overlap of atoms that would no longer allow for converging the electronic structure calculations in the big system II in the x1 → 0 limit due to enormous Pauli repulsion interactions. This potential Urep(rNa+−X) acts on a length scale of only about 0.1 Å away from Na+ and, thus, imprints additional hard-sphere (Percus−Yevick) like positive and negative modulations on the bicanonical x1 = 0 RDF at very short distances. As a result, the corresponding RDF is unity only for ion−water distances beyond something like ≈0.2 Å as observed.

Figure 1. Snapshot from the bicanonical simulation of an aqueous Na+(aq) solution that hosts half a Na+ ion as obtained for constant x1 = 0.5. The oxygen atoms of the water molecules in the first solvation shell are highlighted by red color, while light blue is used for the more distant water molecules. H


Article

Journal of Chemical Theory and Computation We conclude this discussion by highlighting the finding that the Na+ ion, being the solute species in the present case, behaves essentially like a noninteracting free point particle in the x1 = 0 limit of bicanonical AIMD simulations as it should, and that also the other canonical limit, x1 = 1, is perfectly reproduced by the algorithm implemented in CPMD.23 Based on this validation, we move forward to demonstrate in Figure 3 the successive transformation of the hydration shell

Figure 4. Radial distribution functions of the O−O distances obtained from canonical AIMD simulations of the pure solvent using the 32H2O system and of the Na+32H2O system (red line with open and filled circles, respectively) and from bicanonical AIMD simulations using constant weights of x1 = 1 (black line) and x1 = 0 (blue line) corresponding to the Na+(aq) solution and to the pure solvent, respectively.

Figure 3. Radial distribution functions of the Na+−O distances obtained from bicanonical AIMD simulations using constant fractional weights of x1 = 0.25, 0.5, and 0.75 corresponding to a fractional Na+ solute in solution (see text) as well as the two limiting scenarios using constant weights of x1 = 1 (black line with circles) and x1 = 0 (blue line) corresponding to the Na+(aq) solution and to the pure solvent, respectively.

to achieving that goal which are rigorously connected via the underlying effective bicanonical potential as derived in sections 3.3.4 and 3.3.5. 5.2. Solvation Free Energies from Sampling the Bicanonical Ensemble. At this stage, we are ready to apply our bicanonical AIMD method to calculate solvation free energies, serving as our proof of concept that this technique is able to treat open systems subject to a variable number of both nuclei and electrons at the same time. To this end, we have chosen to compute the solvation free energies ΔA(M+) of wellstudied metal ions M+ in aqueous solution, namely Na+(aq) and K+(aq), given that very accurate experimental and computational reference data are available in the literature as comprehensively covered in ref 19. In the following we will explore both sampling schemes that are intrinsically offered by the bicanonical AIMD method, i.e., simulations (i) at constant bicanonical chemical potential μ†1 as described in section 3.3.4 with the solvation free energy given by eq 30 and (ii) at constant fractional force weight x1 according to eq 32 as presented in section 3.3.5. 5.2.1. Solvation Free Energy via Sampling the Bicanonical Chemical Potential. The first route to solvation free energies ΔA(M+) is the sampling of bicanonical ensemble averages of its chemical potential ⟨μ†1⟩x1 at constant fractional force weight x1. Using this approach, we have derived in section 3.3.5 that the bicanonical approach rigorously reduces to thermodynamic integration, AITI, along this continuous constraint variable x1; see eq 32. Thus, the solvation free energy is obtained from

pattern around Na+(aq) when “creating” a Na+ ion solvated by bulk water upon continuously changing x1 from 0 up to 1. Upon increasing the value of x1 starting from zero, i.e., from the neat liquid without any solute species, one can nicely see in the bicanonical AIMD simulations with fractional x1 > 0 how the Na+(aq) solvation shell with its first maximum at around 2.3 Å as well as the subsequent minimum at about 3.2 Å emerge in the solutions that contain a fractional Na+ ion until yielding the standard RDF of a Na+ ion in liquid water when reaching the final value of unity that corresponds to a full Na+ ion in bulk water; see Figure 3. Finally, we assess the oxygen−oxygen RDFs which are characteristic to the H-bonded structure of the solute. This function g(rO−O) is computed for both the neat solvent system and the Na+(aq) solution from bicanonical AIMD simulations in the limits x1 = 0 and 1, respectively, and compared to standard canonical AIMD simulations of these two limiting systems. The comparison in Figure 4 makes clear the point that the water structure of the pure solvent is indeed identical for both methods. Clearly, the water structure in the Na+ solution, as quantified by the oxygen−oxygen RDF, is only slightly different from that of the pure solvent given the total number of water molecules in these supercells. Yet, even this small difference is convincingly reproduced by the bicanonical AIMD simulations carried out at x1 = 1. In conclusion, the bicanonical AIMD method and its implementation in CPMD23 is shown to be able to “insert” solute species consisting of both nuclei and electrons into a pure bulk liquid or, in other words, it is able to ”alchemically generate” such a solvated solute starting from the neat liquid as the initial state. In what follows, we demonstrate that this method allows one to compute the free energy associated with this process, where bicanonical AIMD offers several pathways

ΔATI(M+) =

1

∫0 (⟨μ1†⟩x

1

− μ Mshift + ) dx1

(34)

+

which creates M in pure water upon integrating from x1 = 0 to unity; note that ΔN1 = 1 in the present case. Herein, the x1independent additive term μshift M+ collects all those additional I


Article

Journal of Chemical Theory and Computation contributions that are required to relate our bicanonical sampling result, which only represents the metal ion creation step within a finite sample of pure bulk solvent, to the full solvation process, i.e., (M+)vac + (H 2O)bulk → (M+H 2O)bulk

(35)

where a metal ion is transferred from vacuum into a macroscopic sample of pure water as explained in section 5.2.3. This term thus involves19,34 the total electronic energy of M+ in vacuum, finite size corrections for charged species treated within periodic supercells, as well as the surface potential of the vacuum/water interface as exposed in great depth in ref 19. Concerning the calculation of the surface potential, we follow section 6.2.1 of ref 19, and the finite size corrections for charged solutes subject to periodic boundary conditions are obtained within the so-called B-type correction presented in section 6.2.1 of ref 19 (note that we do not require the C1 correction since the vacuum energy of M+ is computed within finite cluster boundary conditions8 in our case). Following this validated scheme, the strong methodology dependence of the single metal ion solvation free energies computed from atomistic simulations are essentially eliminated,19 which enables us to compare our data to the authoritative literature on ion solvation energies.19 Finally, μshift M+ also contains the kBT ln γ1 term that is specific to sampling the bicanonical ensemble via our effective bicanonical potential. These aspects, however, are not key to the present investigation that focuses on what bicanonical AIMD simulations can offer. Thus, we decided to collect the detailed discussion of the μshift M+ term in section 5.2.3, where we present all individual contributions and also explain how we computed them in the present study. The bicanonical AIMD sampling to obtain the bicanonical ensemble averages ⟨μ†1⟩x1 at constant x1 was carried out for two monatomic species (s = 1), namely for M+ being Na+ and K+ in water. Using molecular dynamics, these ensemble averages were accumulated based on the instantaneous values of μ1†(t) obtained at the respective constant x1 by evaluating eq 31 at each configuration of the bicanonical AIMD trajectories. The results of these simulations are compiled and depicted in Table 1 and Figure 5, respectively, based on a statistics obtained from 82 ps of bicanonical AIMD simulation time for each of these

Figure 5. Bicanonical ensemble averages of the bicanonical chemical potentials ⟨μ†1⟩x1 − μshift M+ obtained from bicanonical AIMD simulations as a function of constant fractional force weight x1 for Na+ (crosses) and K+ (stars) cations in aqueous solution. Solid lines show the second order polynomial fits that are used to integrate analytically eq 34 to compute the full solvation free energies ΔATI(M+) as reported in Table 1.

constant x1 trajectories. The final integration from x1 = 0 to 1 to yield ΔATI(M+) was carried out analytically after fitting second order polynomials to the ⟨μ†1⟩x1 − μshift M+ data which are also displayed in Figure 5. The resulting full solvation free energies ΔATI(M+) for Na+(aq) and K+(aq) solutions are gathered in Table 2 and can be compared directly with the recommended experimental values (being the so-called standard intrinsic hydration free energies as taken from Table 5.6 of ref 19). Table 2. Solvation Free Energies of Single Metal Ions Na+ and K+ in Water from Sampling the Bicanonical Chemical Potentiala bicanonical ensemble AIMD ΔA (M ) ΔAmid(M+) ΔAshoot(M+) ΔsG⊖(M+) TI

0.10 0.20 0.25 0.40 0.50 0.60 0.75 0.80 0.90

M+ = K

−3.07 −3.50 −3.65 −4.15 −4.48 −4.67 −5.09 −5.14 −5.33

−2.52 −2.95 −3.12 −3.48 −3.65 −3.76 −4.06 −4.19 −4.29

M+ = K +

−4.33 −4.48 −4.44 −4.35

−3.54 −3.65 −3.66 −3.60

Results of three different approaches for calculating these energies are given, namely from full thermodynamic integration (TI) along constant but fractional force weight x1 as obtained within the full bicanonical AIMD simulation approach, from bicanonical AIMD sampling at a constant force weight of x1 = 0.5, and from sampling at constant bicanonical chemical potential via the shooting algorithm; see text. Experimental values of the recommended standard intrinsic hydration free energy ΔsG⊖(M+) are taken from Table 5.6 of ref 19. All energies are reported in electronvolts.

⟨μ†1⟩x1 − μshift M+ M+ = Na

eq 34 eq 36 eq 30 experiment

M+ = Na+

a

Table 1. Bicanonical Ensemble Averages of the Bicanonical Chemical Potentials ⟨μ†1⟩x1 Obtained from Bicanonical AIMD Simulations at Constant Fractional Force Weight x1 for Na+ and K+ Cations in Aqueous Solutiona x1

+

The total computer time for obtaining ΔATI(M+) is rather significant because one needs to sample several bicanonical trajectories at constant x1 values to generate sufficiently many quadrature points. However, a very reasonable estimation to our results of ΔATI(M+) is provided by the midpoint approximation ΔAmid (M+) ≈ ⟨μ1† ⟩x1= 0.5 − μ Mshift +

(36)

which is based on a single bicanonical simulation exactly halfway between the initial and final states at x1 = 0.5 (see vertical dotted line in Figure 5 and resulting data in Table 2); note that the rightmost term of eq 31 is strictly zero at x1 = 0.5. The relative deviation of ΔAmid(M+) from ΔATI(M+) is only about 3% for both ions in the present case.

a

The free energies of the complete solvation process (eV) are obtained after accounting for the ion-specific but constant μshift M+ term as discussed in section 5.2.3. J


Article


limiting states, i.e., system I or II, so that one can rather safely stop these simulations after a shooting time of only tshoot = 0.05 ps. Such an ultrafast time scale of commitment can be understood from the functional relationship between x1 and μ†1 according to eq 26, resulting essentially in a Heaviside-like projection onto the two limiting system states I and II at β values corresponding to room temperature; it is noted in passing that we extended a few such shooting runs to very long propagation times but did not find significant deviations from the short-time approximation. By counting how often the trajectories end up either in state I or state II after having reached tshoot, corresponding to x1(tshoot) = 0 or x1(tshoot) = 1, respectively, one obtains ⟨x1(tshoot)⟩μ†1 = ⟨x1(t)⟩shoot from the μ†1 ensemble average over all 1000 shooting trajectories. This data is plotted in Figure 6 including the standard error for each data point (again considering the μshift M+ term).

Before moving on to compare these results to the experimental reference data and the ones computed with bicanonical AIMD simulations at constant chemical potential in section 5.2.2, we discuss in what follows some shortcomings in our approach. The first one is the limited system size, which is only 32 water molecules because of the high demands of computer time due to using AIMD sampling, which is however in line with earlier such simulations. The finite size effect that comes from using such a small number of water molecules has been carefully quantified earlier.19,35 Using periodic boxes that host ≥512 water molecules, benchmark solvation free energies that have been converged up to 0.01 eV in that respect have been obtained using force field based MD (see Table 6.3 of ref 19). Compared to these computational reference data, the solvation free energy obtained using the same approach but with only 32 water molecules turns out to be overestimated by about 0.1 eV (see ibid). It is expected that a similarly small error occurs in the present case when using AIMD instead of force fields to represent the interactions. Upon applying the finite size and surface potential correction schemes via μshift M+ (as explained in section 5.2.3), our calculated solvation free energies correspond to what is called the “intrinsic” solvation free energies (see section 4.2.5 in ref 19). As demonstrated in refs 19 and 35, these energies can be meaningfully compared not only to those obtained with different computational methods but also to experimental solvation free energies (see Table 2) as a result of having corrected for systematic errors which are beyond the computational method as such. In addition, one has to correct by about 0.08 eV for the standard state (see eq 37 of ref 35) when referring to the corresponding standard solvation Gibbs free energies, ΔsG⊖(M+). Comparing now the values listed in Table 2, namely the recommended experimental solvation free energies and our computed values from bicanonical AIMD simulations, we find convincing agreement for both ions within the expected accuracy of about 0.1 eV. 5.2.2. Solvation Free Energy via Sampling the Average Fractional Number of Particles. In the second sampling mode offered by bicanonical AIMD simulations, the bicanonical chemical potential μ†1 is kept fixed and the average fractional number of particles x†1 = ⟨x1(t)⟩μ†1 , as thermodynamically defined within the bicanonical ensemble according to eq 29, is determined and is used to compute solvation free energies by utilizing eq 30 of section 3.3.4. However, an accurate sampling of ⟨x1(t)⟩μ†1 within a single AIMD run is tedious because many insertion and deletion events have to be generated before a reliable bicanonical ensemble average is obtained. As an alternative for sampling the force weight factor x1 at constant bicanonical chemical potential μ†1 we therefore devise the following shooting algorithm. In a first step, we generated at constant temperature T one long bicanonical trajectory at a fixed force weight of x1 = 0.5. After having produced about 0.5 ns of such a bicanonical trajectory for each ionic solution, Na+(aq) and K+(aq), 1000 equidistant snapshots of the trajectories were taken as initial conditions for new MD runs. Independent bicanonical AIMD simulations were started from the positions and momenta of each snapshot where now μ†1 is kept constant at some value as explained below; note that here and in the following the constant shift by −μshift M+ is always taken into account without saying. Analyzing these bicanonical trajectories individually, we find that the system quickly commits itself into one of the two

Figure 6. Ensemble average of the bicanonical fractional number of particles, ⟨x1(t)⟩shoot = x†1, computed via a shooting algorithm (see μ†1 text) as a function of the constant bicanonical chemical potential, μ†1, for Na+ (crosses) and K+ (stars) cations in aqueous solution; the error bars denote the standard deviation σ. The target line ⟨x1(t)⟩shoot μ†1 =0.5 that directly yields the respective solvation free energies ΔAshoot(M+), see text, is marked by a horizontal solid line.

= 0.5 Importantly, the value of μ†1 which yields ⟨x1(t)⟩shoot μ†1 exactly provides us directly with the solvation free energy ΔA(M+), see eq 30, which we denote by ΔAshoot(M+) to distinguish it from the previous estimates. Moreover, the bicanonical chemical potential μ†1 coincides with the usual chemical potential, i.e., μ†1 ≡ μ1, for x†1 = 0.5 according to eq 17, since R ≡ 1 in this case. Thus, bicanonical AIMD provides a rigorous procedure to estimate ΔA based on a simulation at constant chemical potential. In practice, a value of μ†1, which yields ⟨x1(t)⟩shoot = 0.5 to a μ†1 good approximation, is easily found from taking the ensemble average ⟨μ†1⟩x1=0.5 at fixed x1 = 0.5 as a first guess. Subsequently, one computes ⟨x1(t)⟩shoot for a set of constant μ†1 values close to μ†1 the initial guess, which allows one to fit the resulting curve of ⟨x1(t)⟩shoot as a function of μ†1 to determine the final μ†1 value μ†1 that produces 0.5 as depicted in Figure 6; note that we included in Figure 6 these data also in regions far away from the target value that yields ⟨x1(t)⟩shoot = 0.5 where they are not necessarily μ†1 needed for the present purpose. The μ†1 values, which were obtained by this procedure, are the solvation free energies ΔAshoot(M+) for the Na+ and K+ ions in water and are reported in Table 2. Their agreement with the results obtained from the K


Article

Journal of Chemical Theory and Computation corresponding thermodynamic integration approach ΔATI(M+) demonstrates the full functionality of our bicanonical AIMD sampling technique. 5.2.3. Determining the Chemical Potential Shift μshift M+ . In order to calculate the solvation free energy of the full solvation process, where a single metal ion is transferred from vacuum M+(vac) into bulk water M+(aq) across the water/vacuum interface as described via eq 35, we need to account for several additional contributions that come on top of the simulated process of creating a solute species in a small sample of bulk solvent, namely μ Mshift +

=

ES Evac (M+)

−

E Hfs2O(M+)

size of the solvated ion as well as the specific dielectric properties of the solvent, a small contribution of the charge remains in the supercell, which leads to spurious Coulomb interactions. Having applied periodic boundary conditions to the aqueous solutions that host a charged solute, we approximately correct for the resulting finite size errors as follows. Based on the techniques introduced and outlined in refs 19, 39, and 40, we calculated the corresponding finite size correction energy due to charged solutes from E Hfs2O(M+) =

+

+ χ ̌ − kBT ln γ1(M ) (37)

qI 2ξ 2εL

−

⎛ ⎛ + 2 4π R (M ) ⎞ (1 − 1/ε)⎜⎜ ⎜ I ⎟ 2L ⎝ 3 ⎝ L ⎠

qI 2ξ

16π ⎛ RI(M − ⎜ 45 ⎝ L

μshift M+

The resulting so-called chemical potential shift term is a constant for a specific solute−solvent system of a given size and, thus, independent from the specific sampling approach to simulate the process embodied in eq 35. The individual contributions in eq 37 are the total electronic structure energy + of the solute species M+(vac) in the vacuum state EES vac(M ), finite size corrections for representing charged species in the solvated state within periodic boundary conditions EHfs 2O(M+), the surface potential of the vacuum/water interface χ̌ that needs to be crossed upon immersing M+(vac) into the solvent, and the last term removes the contribution of the ficticious ideal gas particles, which was introduced in the reformulation of the bicanonical partition function in eq 21. For an all-encompassing and in-depth discussion as well as detailed analyses of the various contributions to be accounted for in order to correctly calculate solvation free energies (and other parameters), not only from computer simulations but also from experiment, we refer again to the monograph19 of Hünenberger and Reif as well as to their original work.34−37 Therein, chapters 6.3.3 and 6.3.4 also discuss those corrections needed when using ab initio simulation techniques,8 that apply on top of those required when using force fields, which covers the seminal work of Sprik et al. (see e.g., refs 16 and 17). In the remainder of this section, we are going to discuss, term by term, how we computed μshift M+ for the present case. 5.2.3.1. Vacuum State Energy of the Solute Species. The + purely electronic total energy EES vac(M ) of the single metal ion + M (vac) in the vacuum state, i.e., its total energy in the language of electronic structure theory, was calculated by optimizing the electronic structure of an isolated ion. Having employed the same computational method as for the AIMD simulation, we here use finite cluster boundary conditions,38 see e.g., section 3.2.3 in ref 8 for background and details. The + resulting values of the electronic energy EES vac(M ) amount to + −1299.334 and −764.265 eV for M = Na+ and K+, respectively, which depends on the particular pseudopotential and basis set that is used in our plane wave AIMD calculations. 5.2.3.2. Finite Size Correction in the Solvated State. The solvated ion (M+32H2O)bulk within system II introduces a positive unit charge, which interacts via Coulomb interactions with the uniform negative background charge of the same magnitude that is automatically included within the Ewald summation technique due to imposing periodic boundary conditions, see e.g., section 3.2 in ref 8. In stark contrast to the isolated ion in vacuum, M+(vac), the solvated ion M+(aq) is efficiently screened by the high κ dielectric medium provided by the water molecules in the simulation box. However, depending on the finite size of the simulation cell, charge, and

5⎞ )⎞ ⎟ ⎟⎟ ⎠⎠

+

(38)

with the net solute charge qI = +1, the effective cavity radius RI(M+) of the solvated ion in solution, the lattice constant of the cubic supercell L, the Madelung constant for a cubic lattice ξ = −2.837297 (as provided, for example, in ref 39) and the dielectric constant ε = 78.4 of bulk water.41 The value of RI(M+) is obtained here using the well-established approach to calculate consistent ion−solvent interactions (see section 6.2.2 of ref 19). This approach uses RI(M+) = (Rrdf(M+) + RG(M+))/ 2 with the effective ionic (Goldschmidt) radii RG(M+) with values of 0.98 and 1.33 Å for M+ = Na+ and K+, respectively (see set 4r in Table 5.4 of ref 19), and the first maximum in the calculated radial distribution function of the single metal ion wrt the oxygen sites of the water molecules, Rrdf(M+). The values of the latter we have calculated from the canonical trajectories of the aqueous M+(aq) solutions, Rrdf(Na+) = 2.37 Å and Rrdf(K+) = 2.79 Å, respectively (for Na+(aq) see Figure 2). Using these values in eq 38 provides us with the finite size corrections for these solvated ions of EHfs 2O(Na+) = −0.114 eV and EHfs 2O(K+) = −0.159 eV, respectively. 5.2.3.3. Surface Potential of the Water/Vacuum Interface. The so-called intrinsic surface potential (compare section 3.3.4 of ref 19) of the water/vacuum interface, χ̌, was calculated using the same computational setup and parameters as employed in the bicanonical AIMD simulations via averaging the electrostatic potential of the water/vacuum interface as follows. Following the computational approach of ref 42, this interface was modeled by a slab model which contained 108H2O molecules. Periodic boundary conditions were applied to the simulation box of a = b = 12.383 Å and c = 80 Å. A force field MD trajectory using the SPC/E water model was sampled within the (N, V, T) ensemble where the target temperature of 300 K was established by massive Nosé−Hoover chains.32 After an equilibration of 9 ns the system was propagated for 27 ns using a time step of Δt of 0.2 fs. Snapshots of this production run were taken every 0.2 ns to calculate the electrostatic potential. These single-point electronic structure calculations were performed with the identical setup of the electronic structure method, namely using the PBE density functional, as the one employed in bicanonical AIMD sampling (see section 4). The resulting average electrostatic potential, V̅ (z) = ⟨∬ dx dy V(x, y, z)⟩/A with the surface area of the interface A = ab and z being the direction perpendicular to the interface, is depicted in Figure 7. The surface potential χ̌ = −3.50 eV was calculated from the difference of the average electrostatic potential of V̅ aq(z) = −2.56 eV in the pure solvent versus that in L


Article


corresponding to a chemical species, to the reference system. The method has been implemented in the CPMD program package and demonstrated to work by computing the solvation free energies upon dissolving Na+ and K+ cations in water. By explicitly considering pure bulk water and the respective aqueous solutions as the two condensed phase systems that are superimposed according to the effective bicanonical potential, three approaches to run and evaluate bicanonical AIMD simulations are found to provide consistent estimates of the solvation free energies in accord with experimental data. Broad applications in different areas of molecular science are expected in the future, in particular after increasing the efficiency of the bicanonical AIMD simulations upon combining it with MD-based accelerated sampling methods that involve extended Lagrangian concepts.43−46 Clearly, more complicated solvation free energies could be accessed by using bicanonical AIMD, for instance those related to interfacial solvation of chemically complex liquid/solid interfaces. Our method can also be extended to study the deletion and insertion of whole molecules rather than single atoms. The intramolecular degrees of freedom of the inserted molecule have to be added to the bicanonical partition function, the short-ranged steeply repulsive potential must act on all nuclei of this molecule, and the fragmentation of the molecule at low configuration weight can be suppressed by adding force field terms (see ref 16) or by computing the intramolecular interaction on the fly by using the same electronic structure method as in bicanonical AIMD. In heterogeneous catalysis, for instance, so-called polar oxide surfaces display an intimate coupling of their catalytic properties and their coverage with adspecies.47−49 If they are reactants, such as H2 in stepwise reduction reactions, they get consumed along the catalytic cycle and thereby significantly degrade the reductive power of the catalyst.50,51 Keeping their chemical potential constant in the frame of bicanonical AIMD would counteract that effect and allow for realistic computer simulation of these complex catalytic processes at experimentally relevant thermodynamic conditions such as elevated temperatures and pressures. The bicanonical AIMD concept might also provide alternative perspectives on simulating elementary electrochemical processes at electrolyte/electrode interfaces, e.g., to control the external electrode potential52 or to extend the semiconductor defect approach to electrochemistry from the static ab initio thermodynamics viewpoint21 to explicitly including thermal solvent fluctuations and entropic effects.

Figure 7. Surface potential of the water/vacuum interface obtained from electronic structure calculations as described in the text (crosses). The values of the average electrostatic potential V̅ (z) in the vacuum region V̅ vac(z) and in liquid water V̅ aq(z) are depicted by horizontal dashed red lines. The limiting z-values defining the interval used for averaging V̅ vac(z) and V̅ aq(z) perpendicular to the interface are indicated by such vertical lines.

the vacuum region V̅ vac(z) = +0.94 eV. It is noted in passing that is has already been shown that the specific water model used in the force field-based MD sampling of interfacial configurations has only a very small influence on the final result of the computed surface potential.42 5.2.3.4. Contribution from the Bicanonical Partition Function. The last term contributing to μshift M+ according to eq 37 is the term −kBT ln γ1(M+). It accounts for the free energy contribution due to introducing fictitious ideal gas particle species in system I and also includes the corresponding ratio of combinatorial (classical particle permutation) factors as required for rewriting the bicanonical partition function in eq 21 and introducing our effective bicanonical potential UΓ according to eq 24, which is instrumental in order to sample the bicanonical ensemble within AIMD simulations. In the present case of a single metal ion inserted into water, the expression for γ1(M+) as given by eq 22 reduces to the partition function of a noninteracting particle in the ideal gas state, γ1(M+) = Λ1(M+)3/V, with the thermal wavelength of the M+ ion with mass mM+ being Λ1(M+) = h/(2πmM+kBT)1/2. Thus, the term −kBT ln γ1(M+) quantifies the corresponding free energy contribution. With M+ = Na+ and K+ and at T = 300 K the values of −kBT ln γ1(M+) amount to 0.299 and 0.319 eV, respectively.

■

6. SUMMARY AND FUTURE PERSPECTIVES We introduced the bicanonical ab initio molecular dynamics (AIMD) method that allows one to carry out simulations at constant chemical potential or constant average fractional particle number related to a specified chemical species, thus including both nuclei and electrons. The method is based on propagating two coupled condensed phase systems described using electronic structure theory, each containing an integer number of nuclei and electrons. Being a general framework, bicanonical AIMD reduces in the zero temperature limit to what is usually called ab initio thermodynamics (AITD), and therefore, adds entropy and thermal fluctuations to AITD. Moreover, bicanonical AIMD also embodies ab initio thermodynamic integration (AITI) if the fractional particle number is kept constant and used as the continuous switching parameter η = 0 → 1 that transmutes the system from the initial to the final state by adding a set of nuclei and electrons,

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 (0) 234 322 6750. Fax: +49 (0)231 4045. ORCID

Johannes Frenzel: 0000-0003-0458-5424 Bernd Meyer: 0000-0002-3481-8009 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

This work has been supported by the German Research Foundation (DFG) via the Clusters of Excellence EXC 1069 “RESOLV” in Bochum and EXC 315 “EAM” in ErlangenNürnberg. Computational resources were provided in Bochum M


Article


(24) Swope, W. C.; Andersen, H. C. A computer simulation method for the calculation of chemical potentials of liquids and solids using the bicanonical ensemble. J. Chem. Phys. 1995, 102, 2851−2863. (25) Meyer, B. In Computational Nanoscience; Grotendorst, J., Blügel, S., Marx, D., Eds.; Forschungszentrum Jülich: Jülich, 2006; pp 411− 418. (26) Meyer, B. First-principles study of the polar O-terminated ZnO surface in thermodynamic equilibrium with oxygen and hydrogen. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 045416. (27) Shi, W.; Maginn, E. J. Continuous Fractional Component Monte Carlo: An Adaptive Biasing Method for Open System Atomistic Simulations. J. Chem. Theory Comput. 2007, 3, 1451−1463. (28) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (30) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple (Errata). Phys. Rev. Lett. 1997, 78, 1396. (31) Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 7892−7895. (32) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé-Hover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635−2643. (33) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. J. Chem. Phys. 2004, 120, 300− 311. (34) Kastenholz, M. A.; Hünenberger, P. H. Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation. J. Chem. Phys. 2006, 124, 224501. (35) Reif, M. M.; Hünenberger, P. H. Computation of methodologyindependent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions. J. Chem. Phys. 2011, 134, 144103. (36) Kastenholz, M. A.; Hünenberger, P. H. Computation of methodology-independent ionic solvation free energies from molecular simulations. I. The electrostatic potential in molecular liquids. J. Chem. Phys. 2006, 124, 124106. (37) Reif, M. M.; Hünenberger, P. H. Computation of methodologyindependent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones interaction parameter sets for the alkali and halide ions in water. J. Chem. Phys. 2011, 134, 144104. (38) Hockney, R. W. The potential calculation and some applications. Methods Comput. Phys. 1970, 9, 135−211. (39) Hummer, G.; Pratt, L. R.; García, A. E. Ion sizes and finite-size corrections for ionic-solvation free energies. J. Chem. Phys. 1997, 107, 9275−9277. (40) Hünenberger, P. H.; McCammon, J. A. Ewald artifacts in computer simulations of ionic solvation and ion-ion interaction: A continuum electrostatics study. J. Chem. Phys. 1999, 110, 1856−1872. (41) Haynes, W. M., Ed. CRC Handbook of Chemistry and Physics, 96th ed.; CRC Press, 2015. (42) Pham, T. A.; Zhang, C.; Schwegler, E.; Galli, G. Probing the electronic structure of liquid water with many-body perturbation theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 060202. (43) Abrams, J. B.; Rosso, L.; Tuckerman, M. E. Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics. J. Chem. Phys. 2006, 125, 074115. (44) Wu, P.; Hu, X.; Yang, W. λ-Metadynamics Approach To Compute Absolute Solvation Free Energy. J. Phys. Chem. Lett. 2011, 2, 2099−2103. (45) Böckmann, M.; Doltsinis, N. L.; Marx, D. Adaptive Switching of Interaction Potentials in the Time Domain: An Extended Lagrangian Approach Tailored to Transmute Force Field to QM/MM Simulations and Back. J. Chem. Theory Comput. 2015, 11, 2429−2439.

by HPC@ZEMOS, HPC-RESOLV, Bovilab@RUB, and RVNRW.

■

REFERENCES

(1) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (2) Tuckerman, M. E. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: Oxford, 2010. (3) Ç ağin, T.; Pettitt, B. M. Grand Molecular Dynamics: A Method for Open Systems. Mol. Simul. 1991, 6, 5−26. (4) Ç ağin, T.; Pettitt, B. M. Molecular dynamics with a variable number of molecules. Mol. Phys. 1991, 72, 169−175. (5) Lo, C.; Palmer, B. Alternative Hamiltonian for molecular dynamics simulations in the grand canonical ensemble. J. Chem. Phys. 1995, 102, 925−931. (6) Lynch, G. C.; Pettitt, B. M. Grand canonical ensemble molecular dynamics simulations: Reformulation of extended system dynamics approaches. J. Chem. Phys. 1997, 107, 8594−8610. (7) Shroll, R. M.; Smith, D. E. Molecular dynamics simulations in the grand canonical ensemble: Formulation of a bias potential for umbrella sampling. J. Chem. Phys. 1999, 110, 8295−8302. (8) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: Cambridge, 2009. (9) Alavi, A.; Kohanoff, J.; Parrinello, M.; Frenkel, D. Ab Initio Molecular Dynamics with Excited Electrons. Phys. Rev. Lett. 1994, 73, 2599−2602. (10) Marzari, N.; Vanderbilt, D.; Payne, M. C. Ensemble DensityFunctional Theory for Ab Initio Molecular Dynamics of Metals and Finite-Temperature Insulators. Phys. Rev. Lett. 1997, 79, 1337−1340. (11) Tavernelli, I.; Vuilleumier, R.; Sprik, M. Ab Initio Molecular Dynamics for Molecules with Variable Numbers of Electrons. Phys. Rev. Lett. 2002, 88, 213002. (12) Blumberger, J.; Bernasconi, L.; Tavernelli, I.; Vuilleumier, R.; Sprik, M. Electronic Structure and Solvation of Copper and Silver Ions: A Theoretical Picture of a Model Aqueous Redox Reaction. J. Am. Chem. Soc. 2004, 126, 3928−3938. (13) Tateyama, Y.; Blumberger, J.; Sprik, M.; Tavernelli, I. Densityfunctional molecular-dynamics study of the redox reactions of two anionic, aqueous transition-metal complexes. J. Chem. Phys. 2005, 122, 234505. (14) Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox Potentials and Acidity Constants from Density Functional Theory Based Molecular Dynamics. Acc. Chem. Res. 2014, 47, 3522− 3529. (15) Blumberger, J.; Sprik, M. Free Energy of Oxidation of Metal Aqua Ions by an Enforced Change of Coordination. J. Phys. Chem. B 2004, 108, 6529−6535. (16) Sulpizi, M.; Sprik, M. Acidity constants from vertical energy gaps: density functional theory based molecular dynamics implementation. Phys. Chem. Chem. Phys. 2008, 10, 5238−5249. (17) Cheng, J.; Sulpizi, M.; Sprik, M. Redox potentials and pK[sub a] for benzoquinone from density functional theory based molecular dynamics. J. Chem. Phys. 2009, 131, 154504. (18) Leung, K.; Rempe, S. B.; von Lilienfeld, O. A. Ab initio molecular dynamics calculations of ion hydration free energies. J. Chem. Phys. 2009, 130, 204507. (19) Hünenberger, P.; Reif, M. Single-Ion Solvation; Royal Society of Chemistry: Cambridge, 2011. (20) Hummer, G.; Pratt, L. R.; García, A. E. Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100, 1206−1215. (21) Todorova, M.; Neugebauer, J. Extending the Concept of Defect Chemistry from Semiconductor Physics to Electrochemistry. Phys. Rev. Appl. 2014, 1, 014001. (22) Ambrosio, F.; Miceli, G.; Pasquarello, A. Redox levels in aqueous solution: Effect of van der Waals interactions and hybrid functionals. J. Chem. Phys. 2015, 143, 244508. (23) Hutter, J. et al. CPMD V4.1. http://www.cpmd.org/ (accessed Feb. 1, 2017). N


Article

Journal of Chemical Theory and Computation (46) Perego, C.; Giberti, F.; Parrinello, M. Chemical potential calculations in dense liquids using metadynamics. Eur. Phys. J.: Spec. Top. 2016, 225, 1621−1628. (47) Martínez-Suárez, L.; Frenzel, J.; Marx, D.; Meyer, B. Tuning the Reactivity of a Cu/ZnO Nanocatalyst via Gas Phase Pressure. Phys. Rev. Lett. 2013, 110, 086108. (48) Martínez-Suárez, L.; Frenzel, J.; Marx, D. Cu8/ZnO(0001̅) nanocatalysts in response to environmental conditions: Surface morphology, electronic structure, redox state and CO2 activation. Phys. Chem. Chem. Phys. 2014, 16, 26119−26136. (49) Martínez-Suárez, L.; Siemer, N.; Frenzel, J.; Marx, D. Reaction network of methanol synthesis over Cu/ZnO Nanocatalysts. ACS Catal. 2015, 5, 4201−4218. (50) Frenzel, J.; Kiss, J.; Nair, N. N.; Meyer, B.; Marx, D. Methanol synthesis on ZnO from molecular dynamics. Phys. Status Solidi B 2013, 250, 1174−1190. (51) Frenzel, J.; Marx, D. Methanol synthesis on ZnO(0001̅). IV. Reaction mechanisms and electronic structure. J. Chem. Phys. 2014, 141, 124710. (52) Bonnet, N.; Morishita, T.; Sugino, O.; Otani, M. First-Principles Molecular Dynamics at a Constant Electrode Potential. Phys. Rev. Lett. 2012, 109, 266101.

O


Bicanonical Ab Initio Molecular Dynamics for Open Systems - Journal

Recommend Documents