A Cost-Effective Method for Free-Energy Minimization in Complex

Subscriber access provided by AUSTRALIAN NATIONAL UNIV

Condensed Matter, Interfaces, and Materials

A Cost-Effective Method for Free-Energy Minimization in Complex Systems with Elaborated Ab Initio Potentials Carlos Bistafa, Yukichi Kitamura, Marilia T. C. Martins-Costa, Masataka Nagaoka, and Manuel F. Ruiz-Lopez J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00271 • Publication Date (Web): 09 May 2018 Downloaded from http://pubs.acs.org on May 10, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A Cost-Effective Method for Free-Energy Minimization in Complex Systems with Elaborated Ab Initio Potentials Carlos Bistafa,a Yukichi Kitamura, a Marilia T. C. Martins-Costa,b Masataka Nagaoka a,c,d* and Manuel F. Ruiz-Lópezb,e* a

Department of Complex Systems Science, Graduate School of Informatics, Nagoya University,

Chikusa Ku, Furo Cho, Nagoya, Aichi 4648601, Japan b

Laboratoire de Physique et Chimie Théoriques, UMR CNRS 7019, Faculté des Sciences et

Technologies, Université de Lorraine, CNRS, BP 70239, 54506 Vandoeuvre-lès-Nancy Cedex, France c

ESICB, Kyoto University, Kyodai Katsura, Nishikyo-ku, Kyoto 6158520, Japan

d

Core Research for Evolutional Science and Technology, Japan Science and Technology

Agency, Honmachi, Kawaguchi 3320012, Japan e

Future Value Creation Research Center, Graduate School of Informatics, Nagoya University,

Chikusa Ku, Furo Cho, Nagoya, Aichi 4648601, Japan

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 35

ABSTRACT

We describe a method to locate stationary points in the free-energy hypersurface of complex molecular systems using high-level correlated ab initio potentials. In this work, we assume a combined QM/MM description of the system although generalization to full ab initio potentials or other theoretical schemes is straightforward. The free-energy gradient (FEG) is obtained as the mean force acting on relevant nuclei using a dual level strategy. First, a statistical simulation is carried out using an appropriate, low-level quantum mechanical force-field. Free-energy perturbation (FEP) theory is then used to obtain the free-energy derivatives for the target, highlevel quantum mechanical force-field. We show that this composite FEG-FEP approach is able to reproduce the results of a standard free-energy minimization procedure with high accuracy, while simultaneously allowing for a drastic reduction of both computational and wall-clock time. The method has been applied to study the structure of the water molecule in liquid water at the QCISD/aug-cc-pVTZ level of theory, using the sampling from QM/MM molecular dynamics simulations at the B3LYP/6-311+G(d,p) level. The obtained values for the geometrical parameters and for the dipole moment of the water molecule are within the experimental error, and they also display an excellent agreement when compared to other theoretical estimations. The developed methodology represents therefore an important step towards the accurate determination of the mechanism, kinetics and thermodynamic properties of processes in solution, in enzymes and in other disordered chemical systems using state-of-the-art ab initio potentials.

.


2

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction The theoretical study of chemical reactions in solution and in complex molecular environments often requires the knowledge of the solute’s geometry at the stationary points of the free energy landscape of the system (reactants, products, intermediates, transition structures). Since the pioneer applications of implicit continuum models to describe reaction paths in solution,1-2 thousands of new studies have been carried out using those models, providing valuable insights into solvation effects on the kinetics, thermodynamics and/or selectivity of chemical processes. In continuum models (for further details see reviews reported before3-5), it is straightforward to compute the derivatives of the free energy of solvation with respect to the solute’s parameters6-7 and therefore to locate the relevant stationary points. However, due to crude approximations, implicit models are not always sufficiently accurate and large efforts have been made in the last years to develop more elaborated models. In this line, Molecular Dynamics (MD) using QM/MM (Quantum Mechanics / Molecular Mechanics) force-fields has become one of the most popular and successful techniques to study complex molecular systems.8 The price to be paid is a much higher computational cost and hence most applications assume semi-empirical or density functional theory QM methods and MD sampling limited to a few tens or hundreds of picoseconds (though recent algorithms9 have allowed to reach multi-nanosecond time scales). In MD simulations, the calculation of free energy reaction profiles is usually accomplished by the use of ad hoc techniques, which require the prior definition of a distinguished coordinate describing the reaction path, such as metadynamics, umbrella sampling, replica exchange, etc (for a review, see10-12). The FreeEnergy Gradient (FEG) method,13-16 in contrast, focuses on the calculation of the free energy derivatives and the direct location of the stationary points, much like in implicit solvation


3


Page 4 of 35

models. For other approaches related to the FEG method, the reader can refer to the work by Kuczera17 to locate stationary points of peptides using classical force-fields, and to the review by Hu and Yang18 on ab initio QM/MM simulations. According to FEG theory, the force acting on the free energy surface FFE(q), where q represents the solute’s atomic coordinates, is obtained by time averaging the forces on these coordinates along an MD simulation with a solute’s restricted structure:13-16 = −

= −〈

〉

(1)

where the brackets mean a configurational average. In the case of a QM/MM force-field, the potential energy is given by the following expression (for simplicity we do not indicate the explicit dependence with respect to solute’s coordinates):

= + / + / +

(2)

Here, is the QM subsystem wave function, while , and / correspond to the Hamiltonians for the QM subsystem, the MM subsystem and the interaction between them, respectively. The QM/MM interaction includes electrostatic and non-electrostatic terms, usually estimated with a Lennard-Jones potential-type function, which are indicated by the superscript. Note that in general obtaining the derivatives of the potential energy require the calculation of the molecular orbitals coefficients (and configuration interaction coefficients depending on method), but nowadays this is feasible using analytical procedures for many QM approaches. By minimizing though an iterative procedure, stationary points in the free energy surface (minima or saddle points) can be located. Starting from a guess solute’s geometry (which can be the gas phase geometry or the geometry optimized using other solvation models), a MD simulation is carried out, the average forces on the solute are calculated to obtain and a


4

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


new geometry is defined. The process is iterated until the change of the internal structure of the solute and the free energy gradients are smaller than predefined thresholds. The procedure may be computationally costly because several QM/MM MD simulations have to be carried out, and hence using high-level ab initio methods to describe the QM subsystem becomes prohibitive in general. Some approximate methods have been proposed in the literature that assume mean-field theory,19-20 according to which the average force in (1) is estimated as minus the derivative of the energy in an average solvent potential, which is obtained from a classical MD simulation. In the present work, we develop a new method in which thanks to the use of free energy perturbation (FEP) theory,12, 21 free energy forces are computed at very high ab initio levels after an appropriate modification of FEG equations. The procedure relies on the composite approach previously proposed,22-23 which has been already used to obtain accurate free energies in QM/MM simulations. The water molecule in liquid water was chosen to illustrate the methodology and to check its efficiency and accuracy, since many theoretical works have been devoted to describe this system and some experimental data are available as well. Pioneer QM/MM simulations24-26 showed that both the dipole moment and the structural parameters of the water molecule are quite sensitive to the solvation effect and undergo significant fluctuations in the liquid. In section 2, the basic equations of the FEG-FEP theory are developed. Computational details are provided in Section 3 and in Section 4 the results are discussed. Finally, some conclusions and perspectives are drawn in Section 5.


5


Page 6 of 35

2. Theoretical treatment Dual-level free energy calculations Dual-level methods are widely used in electronic structure calculations of isolated systems (usually represented by the double-slash notation A//B) where a single-point energy computation is made at a high-level A on the geometry optimized at a lower level B.27 A related approach was developed for correcting the total potential energy in QM/MM simulations of solutions by Gao and coworkers.28-29 In this case, the total energy of the system includes the gas phase solute’s energy (solvent independent) calculated at a high theoretical level, and the solute-solvent interaction energy computed in the simulation at a lower level. Variants of this QM/MM duallevel method were also developed by Tuñón and coworkers.30-31 More elaborated approaches have to take into account high-level corrections on the QM/MM subsystem interactions too, and to this aim FEP theory can be employed. FEP theory is one of the most useful techniques to calculate the free energy difference between a reference system described by Hamiltonian Ho and a target system described by Hamiltonian H1 = Ho + ∆H, where ∆H represents a perturbation.12, 21 FEP theory can be used, for instance, to calculate free energies of solvation, or the variation of the free energy as a function of a particular solute coordinate. In the context of dual-level QM/MM simulations, FEP theory has been used to estimate the free energy change associated to a change in theoretical level, i. e. to calculate the free energy of the system at a high QM theoretical level (High Level = HL) from MD simulations carried out at a lower QM theoretical level (Low Level = LL). Different approaches have been proposed in the literature and the reader is referred to references22-23, 32-40 and other works cited therein for further details on each technique.


6

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Let be the Gibbs free energy of a solute-solvent system described by the Hamiltonian LL at given thermodynamic conditions and for a given geometry of the QM solute . Let ! be the corresponding free energy when the same system is described by a Hamiltonian at theoretical level HL. We assume here that LL and HL correspond to QM/MM Hamiltonians with the same MM force-field and different QM level, and that the QM and MM subsystems correspond to the solute and solvent molecules, respectively. This is a particular case but extension to other schemes is straightforward. The change from LL to HL is considered to be a perturbation of the system so that FEP allows estimating ! from according to:22-23 "

! = − ln〈$ % & ;( 〉 #

(3)

where the configurational average indicated by the brackets is taken over the LL simulation sampling and Δ; ( is the difference between the HL and LL potential energies (given by equation 2) for a given solvent configuration (: Δ; ( = ! ; ( − ; (

(4)

Equation (3) can easily be deduced from the definition of the free energy in terms of partition functions (see Supporting Information). It states that the free energy difference can be obtained by sampling only equilibrium configurations of the reference state. In practical FEP applications, this means that the unperturbed Hamiltonian (the LL method) must provide a “suitable sampling” for the calculation of the HL free energy correction. In particular, the sampling must be such that the distribution of the HL vs LL potential energy difference Δ; ( in equation (3) is correctly mapped. Some discussion on this issue and techniques to minimize statistical errors can be found in previous reports.23, 32


7


Page 8 of 35

The FEG-FEP method From the FEP expression (3), we can obtain the derivatives of ! with respect to solute’s coordinates q as: *+,

=

*,,

"

− -# ln〈$ % & ;( 〉 .

(5)

Expanding the last term, and after multiplying and dividing by e% & ;( , one easily obtains: *+,

=

*,,

+〈

+, ;(

0; ( −

,, ;(

〉11

(6)

with the weighting factor: 2345;(

0; ( = 〈2345;( 〉

(7)

66

Equation (6) can be further simplified using the FEG equation (1) to: *+,

=〈

+, ;(

0; (〉11

(8)

This is the basic equation in the FEG-FEP method to obtain the HL free energy gradient, which requires the calculation of the HL potential energy ! ; ( and its derivatives,

+, ;(

, on a

selected set of configurations taken from the LL sampling. The following points deserve consideration: (1) If a non-polarizable MM force-field is used, the energy for a given configuration is independent of the QM atoms coordinates, and its derivative with respect to q is zero. In addition, if we assume the same MM force-field for the HL and LL, the MM solvent energy difference cancels out in Δ; ( when calculating w. (2) The HL vs LL energy difference for the non-electrostatic solute-solvent term (typically obtained with a Lennard-Jones potential) does also cancel out in Δ; ( if the same parameters are used for the two QM/MM levels.


8

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(3) It is worth reminding that a constant value can arbitrarily be added to Δ; ( without changing the actual value of the weighting factor w because the same quantity will appear in the numerator and the denominator and will cancel out. This is equivalent to change the reference of the potential energy of the system at HL and (or) LL. Since, in practice, it is suitable to work with small exponents to obtain w, the potential energy difference Δ; ( can be safely changed by the solvation energy difference Δ 7 ; ( defined as: Δ 7 = Δ; (−8!9 −8 9

(9)

where 8!9 and 8 9 are the HL and LL gas phase energies of the solute at optimized geometries. The quantity 8!9 −8 9 is therefore a constant, independent of the configurational sampling. The term Δ 7 contains the solute deformation (geometry change) and polarization energy, and the solute-solvent electrostatic and non-electrostatic interaction energies. If the HL and LL Hamiltonians involve different MM force-fields, an additional energy term for the pure MM liquid should appear in equation (9). (4) Finally, it is important to emphasize that in FEG-FEP calculations the solute-solvent interactions must be comparable in the HL and LL methods for any fixed geometry of the solute, because such terms determine the main features of the sampling. This condition, however, is not necessary for the potential energy of the solute.

3. Computational details Two different case studies have been conducted. In the first one, the aim was to discuss in some detail the numerical aspects of the FEG-FEP methodology, the convergence of the algorithm, errors related to limited sampling, etc. For this reason, the theoretical levels, HL and LL, were not unduly costly, specifically we assumed the HF/6-31G(d) and B3LYP41/6-


9


Page 10 of 35

311+G(d,p), respectively. A standard FEG minimization at the HL was carried out in parallel for comparison purposes and to check the accuracy of the FEG-FEP results. In case study 2, our aim was to reach a much higher accuracy for the structure of the water molecule. Hence, the B3LYP/6-311+G(d,p) was now assumed to be the LL, while QCISD42/aug-cc-pVTZ was chosen for the HL. The QCISD method was chosen rather than the CCSD43 method even though the later is generally considered to be slightly superior to the former.44 However, in the case of water molecule, the gas phase results are very close (see Supporting Information) but QCISD is slightly less costly, and considering the large number of energy and gradient calculations that had to be carried out, the computational time saving is not negligible. For simplicity, in the following we will refer simply to HF, B3LYP and QCISD calculations for calculations at HF/631G(d), B3LYP/6-311+G(d,p) and QCISD/aug-cc-pVTZ levels, and to B3LYP::HF and QCISD::B3LYP, for the two FEG-FEP minimizations described above. The QM/MM MD simulations were carried out using the previously developed interface45 connecting Gaussian 0946 and AMBER 09.47 They were performed in the NVT ensemble at 300 K; the temperature control was made through a Berendsen thermostat.48 We used a cubic simulation box, with a side of 26.25 Å, composed by 1 central QM water molecule surrounded by 592 TIP3P49 water molecules. The Lennard-Jones parameters for the QM water molecule were the same as those for the classical TIP3P waters in both case studies, as this choice provided satisfactory results in previous simulations of “water in water” at comparable computational levels.24 Errors in the effective intermolecular potential related to this approximation might be decreased by re-optimizing the parameters for each QM/MM level,50-52 possibly through an efficient iterative procedure based on mean-field QM/MM techniques.51 However, such a refinement was not contemplated in this work, since the main goal was to


10

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


illustrate the good performances of the FEG-FEP methodology in comparison to the standard FEG calculations, and the same parameters were used in both cases. The QM water molecule was placed in the ZY plane, the Z-axis being the symmetry axis, and its position was constrained during the simulation using the BELLY method implemented in AMBER 09.47 Periodic boundary conditions in the 3 directions were assumed and we used a strict cutoff radius of 11 Å for QM-MM interactions (MM-MM interactions were treated by the Ewald method with the same cutoff, as explained in reference45). The time step was 1 fs. At each iteration of the optimization procedure, we performed a 40 ps simulation, and 10% of the configurations were saved to calculate the gradient averages (4000 snapshots, or one configuration every 10fs). These values were chosen after preliminary convergence tests (see Supporting Information). They provide averages with an acceptable statistical error and preserve the symmetry conditions; for instance the average absolute differences in OH bond lengths in the optimization cycle was less than 2 10-4 Å, i. e. an order of magnitude smaller than the acceptable OH bond length error 10-3 Å. The optimization procedure follows the simple steepest-descent method used in previous

FEG applications.16 Although more elaborated optimization algorithms using the Hessian might be suitable for complex systems (such as the conjugate gradient and BFGS methods),53, in the present case only a limited number of degrees of freedom are involved (the OH distances and the HOH angle) and this technique performs satisfactorily. Note that the necessary equations to calculate the free energy Hessian in the FEG method are available54 but some developments will be required in the case of the FEG-FEP version. The geometry for the QM water molecule at iteration i+1 is obtained from the position at iteration i according to the expression:16 :;" = : + ∆:

(10)


11


Page 12 of 35

where the displacement vector has the form: ∆: = =: ∙ ? " ∙ : @

(11)

Here, =: is an adaptive constant with dimension (time)2, ?

is the diagonal matrix ?AB =

CD EAB , where CD is the atomic mass of nucleus a associated to the k coordinate, and : is the free energy force vector at iteration i. After some preliminary tests, we decided to use c=0.25 fs2 in all iterations, which as shown hereafter led to a smooth convergence of the geometry. Convergence was assumed to be achieved when the RMS displacement was below 0.00025 au (0.00013 Å), the RMS force below 2.5 10-3 au (3 kcal⋅mol-1⋅Å-1), and the predicted free energy change in the steepest descent method: H

δG: = −=: @ ∙ ? " ∙ : : @

(12)

was smaller than 5 10-7 au (0.0031 kcal⋅mol-1). In principle, the free energy change between two successive FEG iterations can accurately be obtained using free-energy perturbation theory16 but the additional computational cost is significant and we did not envisage its calculation in this work. The parameters above are similar to those used in previous FEG applications. One might note that the residual free-energy forces are larger than usual potential energy forces in gas phase calculations but it has been shown16 that this threshold is difficult to overcome, even when sampling is improved through the use of mean-field theory.20 The starting structure of the QM water molecule in the case study 1 was the optimized geometry in gas phase at the HF level (dOH=0.947 Å, αHOH=105.5°), which displays considerably shorter OH bond lengths that the optimized geometry at the B3LYP level in the same conditions (dOH=0.962 Å, αHOH=105.1°), though the HOH angle is quite similar. The gas phase experimental geometry is55 dOH=0.957 Å and αHOH=104.5°. In case study 2, since the computational cost is much higher, we sought to minimize the number of iterations; for this


12

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


purpose, the experimental geometry in the liquid,56 dOH=0.970 Å and αHOH=106°, was chosen for the initial structure of the water molecule in the calculations.

4. Results and discussion The efficiency of the FEG-FEP method is demonstrated below through the study of the water molecule in liquid water. However, before discussing in detail the obtained results in the two case studies, we briefly analyze the adequacy of the LL methods chosen to get the configurational sampling, as well as the statistical errors on the calculated free energy gradients at HL. At the end of this section, some comments will be also presented on the computational cost of the FEG-FEP method in comparison to calculations using standard FEG theory.

Comparison of solvation energy and solute-solvent radial distributions As in any FEP application, one must first of all corroborate that the sampling obtained with the reference Hamiltonian (LL here) is suitable to calculate average properties using the target Hamiltonian (HL here). To get a qualitative picture of this condition, we first make a comparison of the solvation energies and the solute-solvent radial distributions functions obtained with both methods. In the next subsection, we will make a quantitative evaluation of the errors in the FEGFEP results. The top panel of Figure 1 displays the distribution of Δ I9B in the two case studies for the corresponding initial point of each optimization procedure. As shown, both curves exhibit small 7 dispersion. In case study 1, the curve displays a maximum around ΔJDK =0.65 kcal⋅mol-1 and 7 92.9% of the 4000 calculated points lie within ΔJDK ±1 kcal⋅mol-1 (σ = 0.50 kcal⋅mol-1). In case 7 study 2, the Δ 7 distribution is even sharper; the curve has a maximum at ΔJDK =1.54


13


Page 14 of 35

7 ±1 kcal⋅mol-1 (σ = 0.22 kcal⋅mol-1 and 99.8% of the 4000 calculated points lie within ΔJDK

kcal⋅mol-1). The corresponding weighting factors w are plotted in the bottom panel of Figure 1, and as shown they also display a low dispersion (σ = 0.86 and σ = 0.46 for case studies 1 and 2 respectively; w≤8.7 and w≤6.7 for case studies 1 and 2 respectively). The solute-solvent radial distribution functions from LL and HL QM/MM simulations for the starting water geometry in case study 1 are compared in Figure 2. The distributions are very similar and in particular both show comparable solute-solvent hydrogen-bond distances. Such comparison is not possible in case study 2 because, for computational cost reasons, we could not carry out a QM/MM MD simulation at the QCISD level. But considering the even smaller 7 dispersion of the ΔJDK term with respect to case study 1, we can expect the HL and LL radial

distributions functions to be very similar in this case too. One can note, additionally, that the obtained RDF curves compare well with experimental data,57 which exhibit peaks at 2.88 Å (O···O) and 1.85 Å (O···H), close to the results reported in Figure 2. Overall, the energy difference distributions and the radial distribution functions suggest that sampling at LL is adequate to calculate properties at HL, which is a prerequisite to get accurate results with the FEG-FEP method. Quantitative errors are evaluated hereafter.


14

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Normalized distributions of Δ 7 (top panel) and w (bottom panel) in the two case studies for the saved 4000 configurations in the first iteration of the QM/MM MD simulation at LL of each optimization procedure. Δ 7 is the difference between the HL and LL QM/MM solvation energies of the water molecule in liquid water defined by equation (9). The weighting factor w is a scalar function of the solute and solvent coordinates defined by equation (7). Red, open circles: case study 1 (B3LYP::HF). Grey, full circles: case study 2 (QCISD::B3LYP).


15


Page 16 of 35

Figure 2. Radial distribution functions of a water molecule (QM) in liquid water (MM) from QM/MM simulations at HF (red, dashed) and B3LYP (blue, plain) levels in case study 1 (at the initial point of the optimization procedure).

Computational errors in FEG-FEP gradients In this subsection, we will focus on case study 1, for which a comparison between the FEGFEP and standard FEG results is possible. We will focus here on errors associated with the calculation of the free energy gradient correction at HL. Errors associated with the calculation of the zero-order term at LL are inherent to the FEG methodology and have been discussed before,14 so they will not be discussed here. The free energy forces and atom displacements at the initial step in the FEG-FEP and standard FEG calculations, which involves exactly the same QM water geometry, are collected in Table 1. The gradients below refer to derivatives with respect to Cartesian coordinates of the QM water molecule atoms. Average free energy forces at LL and HL, − 〈

LMN ,,

O

〉11 and − 〈

LMN +,

O

〉P1

respectively, are significantly different, as indicated in particular by the large root-mean-square error (RMSE), 10.8 kcal⋅mol-1⋅Å-1. However, when the LL gradient is corrected following our


16

Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


perturbation method (equation 8), the values compare very well with the reference ones so that the RMSE decreases drastically and amounts now only 0.51 kcal⋅mol-1⋅Å-1. Similar comments can be made for the atomic displacements associated with these forces.

Table 1. Analysis of the errors in the calculated free energy forces (in kcal⋅mol-1⋅Å-1) and atomic displacements (in Å) using the FEG-FEP method in comparison to the standard FEG method at HL (B3LYP/6-311+G(d,p)). In both cases, the values correspond to the first step of the optimization (same water geometry). For comparison, the values at LL (HF/6-31G(d)) are also provided. (MAV = maximum absolute value, RMSE= root mean square error, MAE= mean absolute error)

Forces

Displacements

MAV RMSE MAE MAV

RMSE

MAE

HL (standard FEG) 45.5

-

-

0.0034 -

-

HL (FEG-FEP)

46.5

0.5

0.4

0.0034 3.7 10-5 2.6 10-5

LL

26.5

10.8

7.8

0.0017 9.0 10-4 6.0 10-4

The excellent agreement between the FEG and FEG-FEP free energy gradients might be surprising considering their large difference with respect to the starting LL values. One must keep in mind, however, that a significant part of the free-energy gradient correction is associated to the internal potential energy of the solute molecule, which depends on the solute’s geometry but not on the solvent configuration around it.


17


Page 18 of 35

Figure 3 shows the convergence of some parameters along the two independent minimization procedures, the approximate FEG-FEP and the standard FEG ones. As clearly shown in the different panels of this figure, the water OH bond length, the forces, and the predicted freeenergy change in the FEG-FEP iterations follow very closely the values obtained in the (totally independent) standard FEG calculation at HL, demonstrating the adequacy of the perturbation approach. The angle convergence is not shown in this figure because it changes very little from gas phase to liquid phase, and also from LL to HL methods, see below. As a result, the final geometries for the water molecule are practically identical in the FEG and FEG-FEP calculations; they are discussed in deeper detail in the next subsection, where they will be compared to other theoretical data and to available experimental measurements.

Figure 3. Convergence of average OH bond length, RMS Force and predicted change in freeenergy (δG) for the water molecule in liquid water along the approximate FEG-FEP (red, circles, plain line) and standard FEG (blue, triangles, dashed line) optimization methods in case study 1 (HL=B3LYP/6-311+G(d,p), LL=HF/6-31G(d)).


18

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The structure of the water molecule in liquid water The results obtained for the water molecule in liquid water are summarized in Table 2, where they are compared to other QM/MM or ab initio MD simulations and to experimental measurements. Strictly speaking the experimentally measured geometrical parameters correspond to time averages, which may differ from values at the free energy surface minimum. An emblematic case, for instance, is hydrogen peroxide for which the probability distribution of the dihedral angle significantly deviates from a symmetric Gaussian form.9 In the case of the OH bond length and HOH bond angle of water, the deviation is expected to be small (see OH and HH distributions in reference58 for instance) but one should keep in mind that part of the differences reported below might be ascribed to the fact that vibrational (or librational) motions are inherently anharmonic and so asymmetric. All theoretical studies predict the equilibrium OH distance in liquid water to increase with respect to the gas phase. This trend is confirmed by experimental measurements although the liquid state neutron scattering data for the water structure in liquid water have large error bars. The experimentally observed OH bond length is 0.970 ± 0.005 Å, so that the OH bond length increase with respect to gas phase (0.957 Å) would lie between 0.008 and 0.018 Å. Previous theoretical estimations depend on method and model but the most elaborated calculations lead to an OH bond length increase between 0.018 Å and 0.020 Å, which is close to the experimental upper bound limit. The present calculations lead to an OH length increase of 0.020 Å in both case studies, which is equal to the prediction by previous FEG calculations at the MP2/631+G(d) level45 and very close to the ab initio MD simulation result.59


19


Page 20 of 35

The bond angle does not appear to be very sensitive to the solvation effect. Experimentally, an increase by a little more than 1° is foreseen, but again, the error is large. Our calculations predict a slight increase of the angle and in particular the QCISD/aug-cc-pVTZ model leads to a 1.1° angle increase, which is consistent with the experimental trend. Finally, the induced dipole moment of the water molecule in liquid water has been the subject of many discussions in the past and still there is no consensus (see discussion in60). It is an important quantity since it is related to the dielectric properties of water and to its solvation power. In gas phase, the dipole moment is 1.85 D.61 Though there is not a direct measurement of the dipole moment of water in liquid water, the commonly accepted value is 2.6 D, which is taken from the experimental value for the water monomer in ice.62 The experimental estimation for the induced dipole moment is therefore 0.75 D. As shown in Table 2, our calculations are in excellent agreement with this quantity and compare well also with other theoretical estimations gathered in the table or reported in the literature.58, 63-67

Table 2. Calculated properties of the water molecule in liquid water using the standard FEG method and the approximate FEG-FEP method developed in this work. OH distances (Å), bond angles (degrees) and dipole moments (D) are compared to previous calculations and to experimental data. Method

Gas

Liquid

∆

B3LYP/6-311+G(d,p)

dOH 0.962 0.982

FEG a

α

105.1 105.2

+0.1

µ

2.16

+0.75

B3LYP/6-311+G(d,p)

2.91

dOH 0.962 0.982

Ref

+0.020 This work

+0.020 This work


20

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FEG-FEP b

α

105.1 105.2

+0.1

µ

2.16

+0.79

2.95

QCISD/aug-cc-pVTZ

dOH 0.959 0.979

FEG-FEP c

α

104.4 105.5

+1.1

µ

1.86

+0.83

2.69

+0.020 This work

MP2/6-31G+(d,p)

dOH 0.963 0.983

FEG a

α

MP2/aug-cc-pV5Z

dOH 0.958 0.976

FEG d

α

104.3 106.6

+2.3

µ

1.86

+0.86

105.4 105.5

2.72

MP2/aug-cc-pVDZ

dOH 0.966 0.981

FEG d

α

103.7 103.7

Car-Parrinello simulation e dOH 0.972 0.991

DFT/MM simulationsg

+0.018

+0.015

20

67

0.0

+0.019

104.4 105.5

+0.1

µ

1.87

+1.08

2.95

dOH 0.958 0.976 104.8 105.3

+0.018

59

68

+0.5

dOH 0.974 1.001

+0.027

60, 69

0.973 0.978

+0.005

26

104.4 103.3

-1.1

60, 69

105.5 105.0

-0.5

26

2.10

2.92

+0.82

60, 69

2.09

2.77

+0.68

24

2.16

2.54

+0.38

26

α

µ

Experiment h

+0.1

α

α

45

+0.77

µ

CASSCF/ANO-L-VQZP f

+0.020

dOH 0.957 0.970 ± 0.005 +0.013

55-56, 61-62


21


α

104.5 ∼106±1.7

+1.5

µ

1.85

+0.75

2.6

a

Standard FEG optimization

b

Using HF/6-31G(d) as low level

c

Using B3LYP/6-311+G(d,p) as low level

d

Using mean-field theory

e

Average values in the simulation. The BLYP density functional70-71 was used.

Page 22 of 35

f

Values obtained as a perturbation of the gas phase geometry (electric field and field gradient from molecular dynamics simulations) g

Average values in the simulation. In all cases, the Becke-Perdew density functional70, 72 was used together with a triple-ζ quality basis set with polarization functions. See the original references for other details. h

Experimental values for the gas-phase geometry55 and dipole moment,61 and for the liquid water geometry56 and dipole moment of the water molecule. The commonly accepted dipole moment in the liquid is estimated from the dipole moment in ice.33 The bond angle in the liquid has been estimated here from experimental data for OH and HH distances.56

FEG-FEP vs standard FEG computational cost The CPU time reduction in FEG-FEP calculations depends on the LL and HL methods and on the number of configurations used to estimate the HL gradient corrections (10% here). In the present study, the CPU time was about 14% of the standard FEG cost (case study 2). It is worth noting, however, that the CPU cost can be further reduced by minimizing the number of correlated configurations through a judicious statistical analysis.73 For instance, in previous QM/MM Monte Carlo studies74 about only 0.1% of total configurations were used to recover the statistical averages. Hence, we can conclude that by using FEP theory, the computational cost of the FEG method can be reduced by 1-3 orders of magnitude.


22

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Though reducing CPU time is an important advantage of the FEG-FEP methodology, another (possibly more important) one is to allow for a drastic decrease of wall-clock time (WCT) since the expensive calculations (HL gradients on stored snapshots) can be efficiently parallelized. Thus, in our case study 2, the standard FEG calculation at QCISD level would have required 1.2 years of WCT, while the FEG-FEP calculation took roughly 16 days using 40 cores. This WCT is the sum of the time needed to complete the MD simulations at LL in one single core (about 15 days in total) and the HL computations in 40 cores (about 1 day in total). More generally, we can say that the total WCT in FEG-FEP calculations is comparable to the time required for obtaining the sampling at the LL level. This is expected to be a decisive advantage in practical applications of the method.

5. Conclusions The FEG-FEP approach developed here is a composite method allowing to minimize the freeenergy of a molecule in a complex environment at high ab initio levels. It is particularly well adapted to calculations using combined QM/MM force-fields but its extension to other theoretical schemes is not difficult and will be explored in the forthcoming work. The method allows a drastic decrease of the CPU and wall-clock times because only a reduced number of high-level computations have to be done, and they can be carried out in parallel, in contrast to standard free-energy gradient calculations. Errors have been shown to be small provided the lowlevel Hamiltonian is judiciously chosen. The unperturbed Hamiltonian at low-level must be selected so that a “suitable sampling” for the calculation of the high-level free energy correction


23


Page 24 of 35

is obtained, and we have shown that one condition for that is to get an appropriate mapping of the perturbation energies (Δ 7 here), whose distribution should display the smallest possible dispersion, as discussed in other works.23, 32 We have applied the method to study the structure of the water molecule in liquid water since abundant data, both theoretical and experimental, are available in the literature to which our results can be compared. Computations at the QCISD/aug-cc-pVTZ level reported in our work represent probably the most accurate theoretical result available today for this system. We predict an OH bond increase of 0.020 Å, a slight HOH angle increase of roughly 1°, and an induce dipole moment of about 0.8 D. All these results are consistent with the best theoretical estimations previously reported, and with experiment, considering that in the latter there is a large incertitude for the data in the liquid state. The FEG-FEP method opens the door to the accurate study of chemical reactions in solution, in enzymes, and in other complex molecular environments, and this is a promising direction for future work.

ASSOCIATED CONTENT Supporting Information. Development of the FEP equation (3), results for the structure of the water molecule in gas phase with different methods, and details on FEG-FEP convergence tests. This information is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Authors *E-mail: (MN) [email protected], (MFRL) [email protected]


24

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ACKNOWLEDGMENT The authors thank the CNRS and the JSPS for supporting the French-Japanese collaborative project CoSyDy. MFLR is indebted to the Future Value Creation Research Center and to Nagoya University for a Designated Professor position. MFRL and MTCMC acknowledge the French CINES for providing computational resources (Project LCT2550). This work was also supported partially by the Core Research for Evolutional Science and Technology (CREST) of the Japan Science Technology Agency (JST); by a Grant-in-Aid for Science Research from the Ministry of Education, Culture, Sport, Science and Technology (MEXT) in Japan; and also by the MEXT programs “Elements Strategy Initiative for Catalysts and Batteries (ESICB)” and “Priority Issue on Post-K computer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use). The calculations were partially performed using several computing systems at the Information Technology Center at Nagoya University. ABBREVIATIONS FEG, free-energy gradient; FEP, free-energy perturbation; MD, molecular dynamics; QM/MM, quantum mechanics and molecular mechanics; CPU, central processing unit; WCT, wall-clock time; MAV, maximum absolute value; RMSE, root mean square error; MAE, mean absolute error. REFERENCES 1.

Ruiz-López, M. F.; Assfeld, X.; García, J. I.; Mayoral, J. A.; Salvatella, L., Solvent

effects on the mechanism and selectivities of asymmetric Diels-Alder reactions. J. Am. Chem. Soc. 1993, 115, 8780-8787.


25


2.

Page 26 of 35

Pappalardo, R. R.; Sánchez-Marcos, E.; Ruiz-López, M. F.; Rinaldi, D.; Rivail, J. L.,

Solvent effects on molecular geometries and isomerization processes: a study of push-pull ethylenes in solution. J. Am. Chem. Soc. 1993, 115, 3722-3730. 3.

Tomasi, J.; Mennucci, B.; Cammi, R., Quantum mechanical continuum solvation models

Chem. Rev. 2005, 105, 2999-3093. 4.

Cramer, C. J.; Truhlar, D. G., Implicit Solvation Models: Equilibria, Structure, Spectra,

and Dynamics. Chem. Rev. 1999, 99, 2161-2200. 5.

Ruiz-López, M. F., The multipole moment expansion solvent continuum model: a brief

review. In Solvation effects on molecules and biomolecules: Computational methods and applications, Canuto, S., Ed. Springer: Dordrecht, 2008; Vol. 6, pp 23-38. 6.

Tuñón, I.; Bertrán, J.; Ruiz-López, M. F.; Rinaldi, D., Computation of hydration free

energies using a parametrized continuum model. Study of equilibrium geometries and reactive processes in water solution. J. Comput. Chem. 1996, 17, 148-155. 7.

Dillet, V.; Rinaldi, D.; Bertrán, J.; Rivail, J. L., Analytical energy derivatives for a

realistic continuum model of solvation: Application to the analysis of solvent effects on reaction paths. J. Chem. Phys. 1996, 104, 9437-9444. 8.

Rivail, J.-L.; Ruiz-Lopez, M.; Assfeld, X., Quantum Modeling of Complex Molecular

Systems. Springer: Cham Heidelberg New York Dordrecht London, 2015; Vol. 21. 9.

Martins‐Costa, M. T.; Ruiz‐López, M. F., Reaching multi‐nanosecond timescales in

combined QM/MM molecular dynamics simulations through parallel horsetail sampling. J. Comput. Chem. 2017, 38, 659-668. 10.

Abrams, C.; Bussi, G., Enhanced Sampling in Molecular Dynamics Using

Metadynamics, Replica-Exchange, and Temperature-Acceleration. Entropy 2014, 16, 163-199.


26

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11.

Zwier, M. C.; Chong, L. T., Reaching biological timescales with all-atom molecular

dynamics simulations. Curr Opin Pharmacol 2010, 10, 745-752. 12.

Chipot, C.; Pohorille, A., Free energy calculations. Theory and applications in Chemistry

and Biology. Springer: New York, 2007; Vol. 86. 13.

Okuyama‐Yoshida, N.; Nagaoka, M.; Yamabe, T., Transition‐state optimization on

free energy surface: Toward solution chemical reaction ergodography. Int. J. Quantum Chem. 1998, 70, 95-103. 14.

Okuyama-Yoshida, N.; Kataoka, K.; Nagaoka, M.; Yamabe, T., Structure optimization

via free energy gradient method: Application to glycine zwitterion in aqueous solution. J. Chem. Phys. 2000, 113, 3519-3524. 15.

Hirao, H.; Nagae, Y.; Nagaoka, M., Transition-state optimization by the free energy

gradient method: Application to aqueous-phase Menshutkin reaction between ammonia and methyl chloride. Chem. Phys. Lett. 2001, 348, 350-356. 16.

Nagaoka, M., Structure optimization of solute molecules via free energy gradient method.

Bull. Korean Chem. Soc. 2003, 24, 805-808. 17.

Kuczera, K., One‐and multidimensional conformational free energy simulations. J.

Comput. Chem. 1996, 17, 1726-1749. 18.

Hu, H.; Yang, W. T., Free energies of chemical reactions in solution and in enzymes with

ab initio quantum mechanics/molecular mechanics methods. In Annu. Rev. Phys. Chem., 2008; Vol. 59, pp 573-601. 19.

Galván, I. F.; Martín, M.; Aguilar, M., A new method to locate saddle points for reactions

in solution by using the free‐energy gradient method and the mean field approximation. J. Comput. Chem. 2004, 25, 1227-1233.


27


20.

Page 28 of 35

Georg, H. C.; Canuto, S., Electronic properties of water in liquid environment. A

sequential QM/MM study using the free energy gradient method. J. Phys. Chem. B 2012, 116, 11247-11254. 21.

Zwanzig, R. W., High‐Temperature Equation of State by a Perturbation Method. I.

Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420-1426. 22.

Retegan, M.; Martins-Costa, M.; Ruiz-López, M. F., Free energy calculations using dual-

level Born-Oppenheimer molecular dynamics. J. Chem. Phys. 2010, 133, Art. N. 064103. 23.

Martins-Costa, M. T.; Ruiz-López, M. F., Highly accurate computation of free energies

in complex systems through horsetail QM/MM molecular dynamics combined with free-energy perturbation theory. Theoret. Chem. Acc. 2017, 136, 50. 24.

Tuñón, I.; Martins-Costa, M. T. C.; Millot, C.; Ruiz-López, M. F., A Hybrid Density

Functional-Classical Molecular Dynamics Simulation of a Water Molecule in Liquid Water. J. Mol. Mod. 1995, 1, 196-201. 25.

Tuñón, I.; Martins-Costa, M. T. C.; Millot, C.; Ruiz-López, M. F.; Rivail, J. L., A

coupled density functional-molecular mechanics monte carlo simulation method : the water molecule in liquid water. J. Comp. Chem. 1996, 17, 19-29. 26.

Alhambra, C.; Byun, K.; Gao, J., The Geometry of Water in Liquid Water from Hybrid

Ab Initio-Monte Carlo and Density Functional-Molecular Dynamics Simulations. In Combined Quantum Mechanical and Molecular Mechanical Methods, American Chemical Society: 1998; Vol. 712, pp 35-49. 27.

Corchado, J. C.; Truhlar, D. G., Dual-Level Methods for Electronic Structure

Calculations of Potential Energy Functions That Use Quantum Mechanics as the Lower Level. In


28

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Combined Quantum Mechanical and Molecular Mechanical Methods, American Chemical Society: 1998; Vol. 712, pp 106-127. 28.

Gao, J. L., Origin of the solvent effects on the barrier to amide isomerization from the

combined qm/mm monte-carlo simulations. Proc. Indian Acad. Sci. Chem. Sci. 1994, 106, 507519. 29.

Byun, Y.; Mo, Y. R.; Gao, J. L., New insight on the origin of the unusual acidity of

Meldrum's acid from ab initio and combined QM/MM simulation study. J. Am. Chem. Soc. 2001, 123, 3974-3979. 30.

Ruiz-Pernía, J. J.; Silla, E.; Tuñón, I.; Martí, S.; Moliner, V., Hybrid QM/MM Potentials

of Mean Force with Interpolated Corrections. J. Phys. Chem. B 2004, 108, 8427-8433. 31.

Martí, S.; Moliner, V.; Tuñón, I., Improving the QM/MM description of chemical

processes: a dual level strategy to explore the potential energy surface in very large systems. J. Chem. Theor. Comput. 2005, 1, 1008-1016. 32.

Wood, R. H.; Yezdimer, E. M.; Sakane, S.; Barriocanal, J. A.; Doren, D. J., Free energies

of solvation with quantum mechanical interaction energies from classical mechanical simulations. J. Chem. Phys. 1999, 110, 1329-1337. 33.

Štrajbl, M.; Hong, G.; Warshel, A., Ab initio QM/MM simulation with proper

sampling:“first principle” calculations of the free energy of the autodissociation of water in aqueous solution. J. Phys. Chem. B 2002, 106, 13333-13343. 34.

Ming, Y.; Lai, G.; Tong, C.; Wood, R. H.; Doren, D. J., Free energy perturbation study of

water dimer dissociation kinetics. J. Chem. Phys. 2004, 121, 773-777. 35.

Rosta, E.; Klähn, M.; Warshel, A., Towards accurate ab initio QM/MM calculations of

free-energy profiles of enzymatic reactions. J. Phys. Chem. B 2006, 110, 2934-2941.


29


36.

Page 30 of 35

Kearns, F. L.; Hudson, P. S.; Woodcock, H. L.; Boresch, S., Computing converged free

energy differences between levels of theory via nonequilibrium work methods: Challenges and opportunities. J. Comput. Chem. 2017, 38, 1376-1388. 37.

König, G.; Hudson, P. S.; Boresch, S.; Woodcock, H. L., Multiscale free energy

simulations: An efficient method for connecting classical MD simulations to QM or QM/MM free energies using Non-Boltzmann Bennett reweighting schemes. J. Chem. Theor. Comput. 2014, 10, 1406-1419. 38.

König, G.; Mei, Y.; Pickard IV, F. C.; Simmonett, A. C.; Miller, B. T.; Herbert, J. M.;

Woodcock, H. L.; Brooks, B. R.; Shao, Y., Computation of hydration free energies using the multiple environment single system quantum mechanical/molecular mechanical method. J. Chem. Theor. Comput. 2015, 12, 332-344. 39.

Lu, X.; Fang, D.; Ito, S.; Okamoto, Y.; Ovchinnikov, V.; Cui, Q., QM/MM free energy

simulations: Recent progress and challenges. Mol. Simul. 2016, 42, 1056-1078. 40.

Polyak, I.; Benighaus, T.; Boulanger, E.; Thiel, W., Quantum mechanics/molecular

mechanics dual Hamiltonian free energy perturbation. J. Chem. Phys. 2013, 139, Art. No. 064105. 41.

Becke, A. D., Density-Functional thermochemistry. iii. The role of exact exchange. J.

Chem. Phys. 1993, 98, 5648-5652. 42.

Pople, J. A.; Headgordon, M.; Raghavachari, K., Quadratic configuration-interaction - a

general technique for determining electron correlation energies. J. Chem. Phys. 1987, 87, 59685975. 43.

Bartlett, R. J., Coupled-cluster approach to molecular structure and spectra: a step toward

predictive quantum chemistry. J. Phys. Chem 1989, 93, 1697-1708.


30

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


44.

Scuseria, G. E.; Schaefer III, H. F., Is coupled cluster singles and doubles (CCSD) more

computationally intensive than quadratic configuration interaction (QCISD)? J. Chem. Phys. 1989, 90, 3700-3703. 45.

Okamoto, T.; Yamada, K.; Koyano, Y.; Asada, T.; Koga, N.; Nagaoka, M., A minimal

implementation of the AMBER–GAUSSIAN interface for ab initio QM/MM‐MD simulation. J. Comput. Chem. 2011, 32, 932-942. 46.

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman,

J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Gaussian 09, Gaussian, Inc.: Wallingford, CT, USA, 2009. 47.

Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.;

Luo, R.; Merz, K. M., Jr.; ; Pearlman, D. A.; Crowley, M., et al. AMBER 9, University of California: San Francisco, 2006. 48.

Berendsen, H. J.; Postma, J. v.; van Gunsteren, W. F.; DiNola, A.; Haak, J., Molecular

dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684-3690. 49.

Jorgensen, W. L.; Chandrashekar, J.; Madura, J. D.; Impey, W. R.; Klein, M. L.,

Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. 50.

Freindorf, M.; Gao, J., Optimization of the Lennard-Jones parameters for a combined ab

initio quantum mechanical and molecular mechanical potential using the 3-21G basis set. J. Comp. Chem. 1996, 17, 386-395. 51.

Martin, M. E.; Aguilar, M. A.; Chalmet, S.; Ruiz-Lopez, M. F., An iterative procedure to

determine Lennard-Jones parameters for their use in quantum mechanics/molecular mechanics liquid state simulations. Chem. Phys. 2002, 284, 607–614.


31


52.

Page 32 of 35

Stanton, R. V.; Hartsough, D. S.; Merz, K. M., An examination of a density

functional/molecular mechanical coupled potential. J. Comp. Chem. 1995, 16, 113-128. 53.

Kitamura, Y.; Takenaka, N.; Koyano, Y.; Nagaoka, M., Free Energy Gradient Method

and its Recent Related Developments: Free Energy Optimization and Vibrational Frequency Analysis in Solution. In Quantum Modeling of Complex Molecular Systems, Rivail, J.-L.; RuizLopez, M.; Assfeld, X., Eds. Springer: Cham Heidelberg New York Dordrecht London, 2015; Vol. 21, pp 219-252. 54.

Takenaka, N.; Kitamura, Y.; Koyano, Y.; Asada, T.; Nagaoka, M., Reaction path

optimization and vibrational frequency analysis via ab initio QM/MM free energy gradient (FEG) method: application to isomerization process of glycine in aqueous solution. Theoret. Chem. Acc. 2011, 130, 215-226. 55.

Benedict, W.; Gailar, N.; Plyler, E. K., Rotation‐Vibration Spectra of Deuterated Water

Vapor. J. Chem. Phys. 1956, 24, 1139-1165. 56.

Ichikawa, K.; Kameda, Y.; Yamaguchi, T.; Wakita, H.; Misawa, M., Neutron-diffraction

investigation of the intramolecular structure of a water molecule in the liquid phase at high temperatures. Mol. Phys. 1991, 73, 79-86. 57.

Soper, A. K.; Phillips, M. G., A new determination of the structure of water at 25°C.

Chem. Phys. 1986, 107, 47-60. 58.

Shiga, M.; Tachikawa, M., Ab initio quantum mechanical/molecular mechanical

molecular dynamics using multiple-time-scale approach and perturbation theory. Mol. Simul. 2007, 33, 171-184. 59.

Silvestrelli, P. L.; Parrinello, M., Structural, electronic, and bonding properties of liquid

water from first principles. J. Chem. Phys. 1999, 111, 3572-3580.


32

Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


60.

Chalmet, S.; Ruiz-López, M. F., The reaction field of a water molecule in liquid water:

comparison of different quantum/classical models. J. Chem. Phys. 2001, 115, 5220-5227. 61.

Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. J., Dipole moment of water from

Stark measurements of H2O, HDO, and D2O. J. Chem. Phys. 1973, 59, 2254-2259. 62.

Coulson, C. A.; Eisenberg, D., Interactions of H2O Molecules in Ice. I. The Dipole

Moment of an H2O Molecule in Ice. Proc. Roy. Soc. London A 1966, 291, 445-453. 63.

Poulsen, T. D.; Ogilby, P. R.; Mikkelsen, K. V., Linear response properties for solvated

molecules described by a combined multiconfigurational self-consistent-field/molecular mechanics model. J. Chem. Phys. 2002, 116, 3730-3738. 64.

Coutinho, K.; Guedes, R.; Cabral, B. C.; Canuto, S., Electronic polarization of liquid

water: converged Monte Carlo-quantum mechanics results for the multipole moments. Chem. Phys. Lett. 2003, 369, 345-353. 65.

Millot, C.; Cabral, B. J. C., Electronic properties of liquid water by sequential molecular

dynamics/density functional theory. Chem. Phys. Lett. 2008, 460, 466-469. 66.

Tu, Y.; Laaksonen, A., The electronic properties of water molecules in water clusters and

liquid water. Chem. Phys. Lett. 2000, 329, 283-288. 67.

Galván, I. F.; Sánchez, M.; Martín, M.; Olivares del Valle, F.; Aguilar, M., Geometry

optimization of molecules in solution: joint use of the mean field approximation and the freeenergy gradient method. J. Chem. Phys. 2003, 118, 255-263. 68.

Nymand, T. M.; Åstrand, P.-O., Calculation of the geometry of the water molecule in

liquid water. J. Phys. Chem. A 1997, 101, 10039-10044.


33


69.

Page 34 of 35

Chalmet, S.; Rinaldi, D.; Ruiz-López, M. F., A QM/MM/continuum model for

computations in solution: Comparison with QM/MM molecular dynamics simulations. Int. J. Quantum Chem. 2001, 84, 559–564 70.

Becke, A. D., Density-functional exchange-energy approximation with correct

asymptotic behavior. Phys. Rev. A 1988, 38, 3098-3100. 71.

Lee, C.; Yang, W.; Parr, R. G., Development of the Colle-Salvetti correlation-energy

formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785-789. 72.

Perdew, J. P., Density-functional approximation for the correlation energy of the

inhomogeneous electron gas. Phys. Rev. B 1986, 33, 8822-8824. 73.

Coutinho, K.; Rivelino, R.; Georg, H. C.; Canuto, S., The sequential QM/MM method

and its applications to solvent effects in electronic and structural properties of solutes. In Solvation effects on molecules and biomolecules, Springer: Dordrecht, 2008; pp 159-189. 74.

Canuto, S.; Coutinho, K., From hydrogen bond to bulk: Solvation analysis of the n-pi*

transition of formaldehyde in water. Int. J. Quantum Chem. 2000, 77, 192-198.


34

Page 35 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents use only A Cost-Effective Method for Free-Energy Minimization in Complex Systems with Elaborated Ab Initio Potentials Carlos Bistafa, Yukichi Kitamura, Marilia T. C. Martins-Costa, Masataka Nagaoka and Manuel F. Ruiz-López


35

A Cost-Effective Method for Free-Energy Minimization in Complex

Recommend Documents