One-Dimensional Projection of Collective Variables for Effective

Subscriber access provided by UNIV OF DURHAM

Dynamics

One-Dimensional Projection of Collective Variables for Effective Sampling of Complex Chemical Reaction Coordinates Yong Su Baek, and Cheol Ho Choi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01139 • Publication Date (Web): 09 Apr 2018 Downloaded from http://pubs.acs.org on April 9, 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 36 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

One-Dimensional Projection of Collective Variables for Effective Sampling of Complex Chemical Reaction Coordinates

Yong Su Baek and Cheol Ho Choi* Department of Chemistry and Green-Nano Materials Research Center, College of Natural Sciences, Kyungpook National University, Sangyeok, Bukgu, Daegu 702-701, South Korea *E-mail: [email protected] ABSTRACT A one-dimensional projection (ODP) technique is introduced to overcome the challenges in multi-dimensional sampling. With the help of a projected “advancement of reaction ()” control parameter, it was demonstrated that multi-dimensional samplings could be performed with a single parameter, thus dramatically reducing the computational overhead. In our test studies, the ODP technique successfully yielded the free energy surface of the double proton transfer reaction of the acetic acid dimer and the hydrolysis of the methyl diazonium ion. In short, the ODP can be a promising protocol for the samplings of complex chemical reaction dynamics.

1 ACS Paragon Plus Environment

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

I.

Introduction Molecular dynamics (MD) can provide informative free energy surfaces that explicitly

include important aspects of entropic effects. A simple sampling strategy assumes that the evolution of a trajectory will sufficiently sample the configurational space of interest if the simulation time is long enough. However, such an approach is computationally inefficient and undesirable, since transitions between local states are typically rare. Considerable progress has been made in the development of more sophisticated methods to explore a configurational space, which include transition path sampling1, targeted molecular dynamics,2 constrained dynamics,3 umbrella sampling, essential molecular dynamics,4 replica exchange MD,5 integrated tempering sampling, and metadynamics.6 In recent years, ab initio molecular dynamics (AIMD) has emerged as a promising tool for accurate free energy calculations of chemical reactions.7 AIMD in a quantum mechanical/molecular mechanics (QM/MM) framework has been successfully applied to the study of a variety of systems, including isolated molecules, condensed matter, and solid-state systems.8 However, typical configurational transitions or chemical reactions occur on time scales that are significantly longer than those accessible using standard AIMD methodologies. In the specific context of AIMD, the two most popular contemporary free energy methods are umbrella sampling and metadynamics. In umbrella sampling, a series of simulations are performed using a predefined internal degree of freedom. In metadynamics, the system is destabilized along a small set of predefined collective variables (internal degrees of freedom) by adding Gaussian potentials onto the PES in a history-dependent fashion. The successful application of these enhanced sampling methods is dependent on the appropriate definition of a reaction coordinate or a set of collective variables9 (CV), and therefore, requires at least some a priori understanding of the underlying PES. Since the CV typically is a subset of


Page 2 of 36

Page 3 of 36 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


variables, its free energy is obtained by the marginal probability distribution. Therefore it is dependent on the choice of CV. The CV approach has been extremely successful due to its ability to sample large regions of the configurational space.10 In addition, no assumption is made on the final state. However, it is usually difficult to define a set of CV that can take into account the entire relevant configurational space, especially in the case of large biological systems. In such cases, the number of CVs becoming too big, requiring a large computational effort. In addition, its interpretation of the configurational space would be also difficult. To overcome the abovementioned challenge, Brandaurdi et al.11 introduced two variables that can describe the position of a point in configurational space relative to a preassigned path. With these variables, they could combine the CV approach with metadynamics, allowing global searches in the reaction path space. However, the two new variables are constructed based on trajectories. Therefore, simulations are still needed to explore a multidimensional space specified by the CVs. Additional reparameterization is required during simulations to optimize the reaction path. To further improve the sampling efficiency, it is necessary to reduce the explicit multi-dimensional searches while keeping the configurational space sufficiently large. Arguably, the umbrella sampling (US) method12 would be very efficient if the control parameter can be well established. However, because of the high dimension of the CVs, application to conventional umbrella sampling methods remains difficult, especially in the context of QM/MM simulations. To overcome this issue,13 an efficient approach combining umbrella sampling (US) with integrated tempering sampling (ITS) has been introduced. The US method is applied to the main pathway, which is chemically more relevant, and the ITS method is employed to facilitate the sampling of other degrees of freedom. This combined method significantly improved the sampling efficiency in complex reaction systems, but it



also increased the simulation time per window. Another approach involves introducing flexible coordinates that are redefined on the fly in each MD step to reduce the sampling dimensions. Our group successfully introduced flexible asymmetric coordinates (FAC) for reactions involving multi-protonated acids, where the dissociating proton can be reabsorbed by the other oxygen of the same molecule.14 Although the flexible coordinate approach can act as a good control parameter for umbrella sampling, it is restricted to special cases. When the system is too complicated, high-dimensional umbrella sampling cannot be avoided, with significantly increased computational overhead. Herein, we introduce an efficient strategy to address these issues by redefining the reaction coordinates, which effectively reduces highdimensional sampling into a one-dimensional problem in conjunction with umbrella sampling. II.

Theoretical Background and Computational Details

a) Definition of Reaction Coordinates As discussed above, one of the most important elements for successful sampling is the correct definition of reaction coordinates. Even if the reactant and product are well known, defining a reaction coordinate between two minima is challenging. Examples of reactions coordinate are presented in Figure 1. For NaCl dissociation (1a) in solvents, a single distance coordinate, , is sufficient for US simulations. An asymmetric coordinate, , defined as the difference between bond distances, is potentially a good reaction coordinate for the intramolecular proton transfer of glycine (1b). However, the intermolecular double proton transfer of the acetic acid dimer requires two asymmetric coordinates (1c), which requires at least two-dimensional coordinates and significantly increases the computational overhead. As discussed in the introduction, a flexible asymmetric coordinate, , can be adopted for multi-protonated acid systems (1d), for which US sampling along a particular proton transfer coordinate is complicated by the additional proton transfer channel as compared to the case of 4 ACS Paragon Plus Environment

Page 4 of 36

Page 5 of 36 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


a single-proton system. In a carboxylic acid, the two oxygen atoms are in principle equivalent, and all hydrogen atoms near the carboxyl group in the QM region are possible candidates for the proton transfer. Thus, the sampled O–H atomic pairs in the asymmetric coordinate in the umbrella window are interchangeable. Considering all possibilities, the flexible asymmetric coordinate (FAC) is defined as =

√

− ∗ ,

(1)

where the subscript i denotes the atomic index of the deprotonating donors and j represents all possible proton candidates. The superscript min indicates the minimum distance among the sets. Let us assume that an acetic acid and a water molecules are in the QM region of the QM/MM model, as shown in 1d. The two oxygen atoms of acetic acid are the donors, while one carbonyl hydrogen of acetic acid and two hydrogen atoms of the QM water molecule are the possible transferring protons. This particular definition of the reaction coordinate (FAC) effectively removes problematic re-adsorption reactions, and thus, uniquely defines the proton transfer path. For methyl diazonium hydrolysis (1e) by solvent water, a good definition of reaction coordinates is not obvious. As the molecular system becomes sophisticated, the number of relevant internal coordinates naturally increases. Consequently, the computational overhead of umbrella sampling becomes prohibitively expensive because of the corresponding multi-dimensional sampling. b) Linear Synchronous Transit (LST) Path The asymmetric coordinate ( ) approach reduces the two-dimensional sampling of two distance coordinates to a one-dimensional problem. Likewise, the flexible asymmetric coordinate ( ) approach also significantly reduces multi-dimensional sampling problems to one dimension. These coordinates not only reduce computational overhead but also allow sampling specific regions of the free energy surface that are otherwise too difficult to access. 5 ACS Paragon Plus Environment


Therefore, reducing high dimensional sampling systems to low- or one-dimensional ones is the main key to the efficient sampling. We thus propose a transition state searching technique based on conventional ab initio calculations. The particular method of interest is the LST (Linear Synchronous Transit) approach, which is a classical multi-dimensional searching method for transition states.15 A model LST path is shown in Figure 2, which is drawn between the two limiting structures, R and P. The LST path is defined by the condition that all coordinates (Cartesian, internal, or distance matrix) vary linearly between the two points R and P, and thus, confine any intermediate points to the LST path. As illustrated in Figure 2, the LST typically does not pass through the real transition state. Assuming that the transition state is displaced laterally from the LST maximum, one can further laterally optimize the LST maximum under the constraint that the refined structure retains the same relative distance between R and P in the geometric space. The success of the LST method as well as the much improved QST (Quadratic Synchronous Transit)15b approach is based on the availability of an analytical gradient at each point of the geometric space. However, the corresponding analytical gradient on the free energy surface is not available during MD simulations, which is the main obstacle to the adoption of an LST approach in molecular dynamics simulations. c) One-Dimensional Projection of Collective Variables Although MD simulations cannot produce analytical gradients of free energy on the fly, exploring the free energy surface presents some advantages over transition state searching methods. Unlike the LST approach, which requires the exact geometric structures of a given system along the path, the free energy surface (FES) is determined by local probability histograms of the structures. Therefore, only a handful of reaction coordinates are needed as control parameters, as long as the control parameters can faithfully connect the reactant () and product ( ), as seen in Figure 1. Therefore, if one can consistently define a general 6 ACS Paragon Plus Environment

Page 6 of 36

Page 7 of 36 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


control parameter of “the advancement of reaction” (), one can reliably obtain the FES along the . Since the FES is dependent on the choice of reaction coordinates, one should design a systematic way of defining reaction coordinates. A geometric structure, or set of collective variables (!), between the reactant ( ) and product ( ) is presented in Figure 3. !, , and can be represented by either Cartesian, internal, or distance matrix coordinates. However, unlike the " and # of the LST method, where all the vibrational degrees of freedom of the molecular structures need to be specified, only a handful of coordinates are necessary to define (Figure 3). Therefore, the vector !, representing a structure, does not need to lie on , a perpendicular vector the − . Using geometric vectors, we define a tangent vector $ !% , and “the advancement of reaction ()” as = & − '/) − ) $

(2a)

) = )&! − ' ∙ $

(2b)

!% = &! − ' − $

(2c)

Since the maximum length of the reaction path is ) − ), lies between 0 and ) − ). represents the one-dimensional direction of the reaction. Consequently, The tangent vector $ the remaining vibrational degrees of freedom are contained in !% , which requires special attention when the system size increases. The represents the magnitude of ! projected onto the direction of the reaction and can therefore be interpreted as the advancement of the reaction. For convenience, is redefined such that its values lie between 0 and 1 according to the following expression: )/) − ) . = )&! − ' ∙ $

(2d)

Due to the particular definition, = 0 and 1 correspond to the reactant and product, 7 ACS Paragon Plus Environment


Page 8 of 36

respectively. During simulations, it can sample beyond these two limits. Therefore the actual value of can be other than the values between 0 and 1. Our procedure appears to have some resemblances with the method by Brandaurdi et al.11, since the two methods take the path defined by only the initial and final states, and utilize only a subset of coordinates. However, there is a significant difference in the reaction coordinates. Brandaurdi et al. expressed it parametrically in the form of R(t), in which the parametrization is such that R(0)= RReactant and R(1)=RProduct. The initial R(t) is guessed and is optimized by rather involved reparameterization scheme during simulations. On the other hand, our general control parameter is determined by a simple vector operation, which does not require reparameterization. In principle, the one-dimensional direction could be combined with various sampling techniques. In the current paper, due to its simplicity, the umbrella sampling was chosen for convenience. Correspondingly, the center of the umbrella sampling window (+ ) can vary along the . With these definitions, the biasing one-dimensional harmonic restraint potential in each US window can be simply written as

, = - . − + / ,

(3)

where the - is the biasing force constant for the i-th window. The is made dimensionless, since the geometric vectors of ! , , and can contain both bond distance and angle coordinates. Therefore, the dimension of - is kcal/mol. For the intermolecular double proton transfer of the acetic acid dimer shown in Figure 1c, !, , and consist of two asymmetric coordinates. In fact, it is desirable that the control parameter is as simple as possible to provide more freedom to the MD simulations. By defining the control parameter of US as expressed in eq. (2d), the two-dimensional problem of intermolecular double proton transfer in the acetic acid dimer reaction (Figure 1c) can be reduced to a one-dimensional sampling problem. This overall procedure is called the “one-dimensional projection (ODP)” 8 ACS Paragon Plus Environment

Page 9 of 36 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


procedure. As a result, the histogram along the control parameter is obtained naturally.11 d) Various Scenarios of Reaction Paths It should be emphasized that the biasing potential is only applied to , which represents . In other words, there is no the advancement of the reaction along the directions of $ restriction on the perpendicular direction along !% , implying that the minimum free energy points along !% at each are simply determined by the local histogram in each window. Therefore, whether the true free energy minimum along the reaction corresponds to the linear path (Figure 4a) or deviates from it (Figure 4b), the ODP procedure can faithfully locate it. Issues arise when there is more than one path, as shown in Figure 4c (paths A and B). Suppose there is a barrier C separating paths A and B. If the barrier is low, the ODP procedure can sample both paths. However, if the barrier is high, the ODP procedure will most likely only find the lower one. However, by investigating the perpendicular directions of the reaction path (!% ), multiple hidden paths can be also detected if the barriers between paths are low. However, the ODP procedure may not find any meaningful paths, implying that the initial reactant and product states are not correct, and one should then design new reactant and product states, for which a good chemical insight is necessary. With the degrees of freedom along the reaction path being reduced to one, the remaining degrees are contained in the perpendicular directions of the reaction path (!% ). Therefore, as the number of degrees of freedom of the reaction increases, the number of degrees of freedom along !% also significantly increases to the point where the deviation from the linear path becomes much greater than the length of the reaction ()!% ) ≫ ) − )), as shown in Figure 4d, implying that there are several different conformations along !% . This deviation can also arise if the free energy surface along !% is too flat. In this particular case, 9 ACS Paragon Plus Environment


sampling becomes inefficient because too much computational time is spent sampling unnecessary configurational space. An efficient approach to circumvent this issue is to introduce a constant biasing potential along !% , so that the MD simulation only samples the vicinity of the linear reaction path between and . Alternatively, one can divide the path into two smaller sections, 1 − and − 1, or one can also impose mild restraints on the relative positions of the molecules in the reaction. III.

Applications

Two examples were adopted to study the efficiency of the ODP technique: double protonation of the acetic acid dimer and hydrolysis of the methyl diazonium ion. The double protonation of the acetic acid dimer requires only two asymmetric coordinates, as shown in Figure 5 (a). A spherical system of QM molecules surrounded by 398 MM water molecules with a radius of 14 Å was prepared for QM/MM-MD simulations. The QM and MM subsystems were calculated at the HF/6-31G(d) level of theory and with TIP5P, respectively. To prevent water evaporation during long simulations, we applied a harmonic restraint potential with a force constant of 2.0 kcal/mol/Å2 for the boundary solvent molecules, as implemented in CHARMM.16 The center of mass of the solute was also constrained to remain in the center of the sphere with a restraint force of 1.0 kcal/mol/Å2. To equilibrate the systems QM/MM-MD NVT simulations with umbrella sampling potentials were performed on each window at 300 K for 30 ps. Then, QM/MM-MD production simulations were performed for 100 ps on the final structures obtained after equilibration. The free energy surfaces were obtained by the weighted histogram analysis method (WHAM)17. The bin size, tolerance, and temperature for the analysis were 0.02, 0.0001, and 300 K, respectively. a) Double Protonation of the Acetic Acid Dimer


Page 10 of 36

Page 11 of 36 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 two asymmetric coordinates ( , ) are defined in Figure 5(a). They represent the transfer of H1 from O1 to O2 and the transfer of H2 from O3 to O4, respectively. The ODP coordinate ( ! ) is defined from the two asymmetric coordinates as ! = 2 Correspondingly, we have defined the reaction and product vectors as = 20.451 and = 2−0.451

3 .

0.4513

−0.4513, respectively. The reaction path () along − was divided

into 28 US umbrella sampling windows with a biasing force constant of 56~70 kcal/mol. The one-dimensional FES along obtained by WHAM is presented in Figure 6. The FES has a symmetric free energy curve with two transition states, TS1 and TS2, and a shallow intermediate I. The forward reaction barrier is calculated to be 11.9 kcal/mol. From intermediate I, it can be considered that the double proton transfer reaction occurs via a stepwise path. One of proton transfers to the oxygen atom of the other molecule via TS1, forming a weakly bound Zwitterionic intermediate state, [CH3C(OH)2+⋅⋅⋅⋅CH3COO−], followed, by the second proton transfer to the newly deprotonated oxygen via TS2. To further explore the FES near the intermediate, a two-dimensional FES with reaction path and perpendicular direction )!% ) was generated and is presented in Figure 7a. The reaction path from Figure 6 is included as a black solid line. The intermediate appears at (, )!% )/ = .0.5, 0.5/, indicating a slight deviation from the − line and that the ODP procedure can faithfully follow the lowest free energy path. Two-dimensional umbrella sampling with the and coordinates was also performed for comparison, and the resulting FES is shown in Figure 7 (b). The two coordinates ( , ) were divided into 49 umbrella sampling windows from (−1.0 Å, −1.0 Å) to (1.0 Å, 1.0 Å) with a biasing force constant of 100 kcal/mol/Å2. Intermediate states appear at ( , ) = (−0.3, 0.3) and (0.3, −0.3) and exactly correspond to the intermediate shown in Figure 7a. It is noted that the absolute value of )!% ) does not distinguish the two identical reaction paths of these two possible 11 ACS Paragon Plus Environment


intermediates. In addition, the forward energy barrier in Figure 7b is 12.0 kcal/mol, which is nearly identical to the value obtained with the ODP method. Therefore, the two-dimensional umbrella sampling approach can be effectively replaced by the one-dimensional ODP procedure for this particular system. b) Hydrolysis of the Methyl Diazonium Ion The hydrolysis of the methyl diazonium ion (CH3N2+) was also studied with the ODP procedure. The overall SN2 reaction can be summarized as 2H2O + CH3−N2+ → H3O+ + CH3−OH + N2, where one of the water molecules acts as a nucleophile by transferring its hydrogen atom to another solvent water molecule. As a result, at least two QM solvent waters are needed, which requires several internal coordinates for the description of the reaction path. In the simulations, the QM molecules comprise the methyl diazonium ion and two water molecules, surrounded by 398 MM water molecules. The coordinates of the nucleophilic attack ( = 9: : ), bond breaking ( = ;: : ), as well as proton transfer (< ) should be included in the ODP, as shown in Figure 5b. They are essential for the SN2 reaction. However, additional internal coordinates are necessary to align the reaction path along the − direction. Since the introductions of additional coordinates are rather arbitrary, we have prepared two sets of ODP coordinates. They are designated as ! and != (see Figure 5b). The first ODP ! includes two distances, 9> : and 9>;: , in addition to the three main coordinates, while the second set != adds one angular coordinate, ?;: : 9: . The ! and != values of the reactant and product are listed in Figure 5b. Those values are obtained from the gas-phase structure optimized at the ab initio MP2/CCT level of theory. Although the simulations were not particularly sensitive to the exact structures of the two states, accurate 12 ACS Paragon Plus Environment

Page 12 of 36

Page 13 of 36 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


geometries were determined. The overall reaction paths with the two ODP ! and != were divided into 25 windows, with biasing force constant of 56 kcal/mol for both ! and != . To prevent sampling far from the reaction path, a perpendicular restraint potential of 8.4 kcal/mol along !% was also applied for the two simulations. To study the effect of the perpendicular restraint potential, a smaller restraint potential of 5.6 kcal/mol. was applied to ′ the ODP != for comparison, which shall be denoted as ODP != . The resulting one-

dimensional FESs are shown in Figure 8. Sato et al.18 performed full ab initio MD simulations of the hydrolysis of the methyl diazonium ion using the FMO-MD methodology. However, owing to the limited simulations, they only reported snapshots rather than statistical FESs. Despite differences in the curvature, the FESs obtained from the ODP simulations with ′ the ! , != , and != settings are consistent with each other, indicating that the ODP

procedure is rather insensitive to the choice of reaction coordinate and perpendicular restraint. A transition state appears at = 0.15 with a forward free energy barrier (R TS) of ~3 kcal/mol for all three simulations. To analyze the transition state geometries, the = 9: : and = ;: : histograms, the proton transfer coordinate < , and the O1⋅⋅⋅⋅⋅C1⋅⋅⋅⋅⋅N1 angle in sampling windows near = 0.15 were calculated and are presented in Figure 9. The following analysis is based on the simulation with ! . The breaking = ;: : bond distances are the most narrowly distributed, with a maximum value of 2.14 Å. On the other hand, the newly forming = 9: : distances have a broad distribution of 1.5-2.5 Å. After careful observations, two major bond distances of 1.7 and 2.2 Å are identified, which suggests two types of bond-forming structures. The corresponding structures from the trajectories are presented in Figure 10(a) and (b), respectively. The former and latter can be 13 ACS Paragon Plus Environment


designated as tight and loose structures. According to the statistics, the loose structures are dominant and are consistent with the tight structure reported by Sato et al.18 However, their reported loose structure with = 2.680 Å and = 2.683 Å was not found in the ODP simulation. It should be pointed out that their simulations were performed at 700 K, a high temperature that may affect the overall reaction dynamics. The proton transfer coordinate < has a broad distribution with a maximum value at 1.0 Å. The positive < value indicates that the proton is still bonded to the nucleophilic oxygen O1 in the transition state, as shown in the snapshots in Figure 10. Therefore, the proton transfer to the other solvent water occurs much later during the reaction. The O1⋅⋅⋅⋅⋅C1⋅⋅⋅⋅⋅N1 angles lie between 150° and 180° with a maximum at 174°, indicating a nearly linear configuration. After the transition state, a large free energy release of ~40 kcal/mol is observed until = 0.7, and a mild energy release is observed after = 0.7. According to our simulations, the flat free energy surface near = 0.7 is due to the proton transfer. To analyze the structures at = 0.7, the histograms of the proton transfer coordinate of < were obtained, as presented in Figure 11. Two structures with < = ± 0.3 exist in the proton transfer region. These bond distances are much shorter than those in the TS (Figure 9a), indicating that the proton transfer is underway. After the proton transfer, a free energy release of ~4 kcal/mol is observed. To study the structural change of the solvent along the entire reaction pathway, the radial distribution functions (RDFs, g(r)) of O1, C1, and N1 with solvent oxygen (Ow), were computed in the reactant (R), transition state (TS), proton transfer (PT), and product (P) windows. The results are presented in Figure 12. For the oxygen (O1), the first solvation peak around 2.7 Å for all windows indicates that the solvent water is consistently bound to O1. For the carbon (C1), the first solvation peaks of the R and TS structures appear at 3.5 Å. However, the position is shifted to 3.8 Å for PT and P, indicating the dissociation of hydronium and 14 ACS Paragon Plus Environment

Page 14 of 36

Page 15 of 36 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


methyl alcohol. For the nitrogen (N1), the first solvation peaks are distributed around 3.5 Å, except for R, which appears at the shorter distance of 3.0 Å, indicating that N1 is well solvated as a methyl diazonium ion. However, once the = ;: : bond is weakened at the TS, N1 becomes less solvated, implying that the molecular nitrogen character already appears at the TS. Overall, it is concluded that the formation of molecular nitrogen and hydronium occurs sequentially. To analyze the free energy surface perpendicular to the reaction path , twodimensional plots of )!% ) and from simulations with ! are presented in Figure 13. The one-dimensional reaction path from Figure 9 is included as a black solid line. The value of )!% ) represents the deviation from the theoretical reaction path of − . The black solid line lies between 0.5~2, indicating that the true reaction path is somewhat different from the − trajectory. This discrepancy arises from the initial internal coordinates for the reactant and product, as shown in Figure 5. Even though the values are not accurate, the MD simulations can nonetheless yield consistent FESs, as shown in Figure 8. Therefore, a small variation in the reference structures would not affect the final results. Large curvatures in the black solid line along )!% ) are observed, especially near points A and B in the figure, indicating that the true reaction path is non-parallel to the − direction. Nonetheless, it is emphasized that the overall samplings were efficiently performed. c) Convergence Tests of ODP The equation of 9 (see the Figure 14) in the paper by Mones and Csányi19 were adopted to study the convergence behavior of our ODP. The is the general control parameter of ODP and @A.B//@ and @A.BC //@ are the free energy derivatives at time t and at the end of the simulation, respectively. This equation estimates the overall deviation of the free 15 ACS Paragon Plus Environment


energy derivatives at time t compared to the final one at the end of the simulation. The results of acetic acid as well as the three coordinate sets of methyl diazonium simulations are presented in the Figure 14. Overall, it is seen from all the simulations that the deviations significantly decrease upto 20 ps and then monotonically reduced afterwards. At about 80ps, the deviations become below 1 kcal/mol, which is much better than the Dis and DDis coordinates by Mones and Csányi, even though the molecular systems studied here are much more complex than theirs. The ODP convergency is as good as the energy-based reaction coordinates by Mones and Csányi. It is also pointed out that although there are some performance differences, the three choices of coordinates for methyl diazonium exhibit nearly parallel converging behaviors, indicating that the final performance is rather independent on the choice of the coordinates. It is also seen that the convergence performance is quite significantly improved by adding the angular coordinate, ?;: :9: in != as compared to ! to the degree that it converges better than the less complex system of acetic acid simulation. Therefore, one can achieve a good convergence performance with properly designed geometric vector, even if the system becomes much more complex. As shown in Figure 12, the water structure does rearrange in the course of chemical reactions, which can lead to hysteresis problems or very long convergence times. However, it does not appear to deteriorate the convergency in our test systems. Nonetheless, such effects can appear in the much more complex systems, for which the energy-based coordinates can be useful.19-20 IV.

Conclusion

We introduced an ODP technique that can reduce multi-dimensional sampling to one dimension with a single parameter, which represents the reaction advancement along the direction of the reaction. We demonstrated that can be a good control parameter for 16 ACS Paragon Plus Environment

Page 16 of 36

Page 17 of 36 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


umbrella sampling in our test systems. The one-dimensional FES of the double protonation reaction of the acetic acid dimer, a typical two-dimensional sampling problem, could thus be obtained. Specifically, the ODP procedure identified a very shallow zwitterionic intermediate, demonstrating its ability to reproduce such subtle surfaces. In the methyl diazonium ion hydrolysis, which is a very difficult multi-dimensional sampling problem, the ODP simulations consistently reproduced the one-dimensional FES, regardless of the internal coordinates, indicating that the ODP can be applicable to such a complex chemical reaction with a reasonable set of internal coordinates. However, one should be careful in applying the ODP technique to other systems. Our method is limited to cases where one can assume the reactant and product states. However, if this assumption fails, one should redesign the reactant and product states. It is not guaranteed that one will always find the lowest path; therefore, trial and error rounds are to be expected when one studies complex free energy surfaces. Although the validity of the ODP procedure should be further evaluated, it is a promising approach for exploring the free energy surface of complex reaction paths, which typically occurs in condensed phase chemical reactions.

Acknowledgments. This work was supported by the Samsung Science and Technology Foundations (SSTF-BA1701-12)



FIGURE CAPTIONS Figure 1: Umbrella sampling coordinates for the (a) NaCl dissociation, (b) intramolecular proton transfer of glycine, (c) intermolecular double proton transfer of the acetic acid dimer, (d) intermolecular proton transfer of acetic acid to water, and (e) methyl diazonium hydrolysis. Figure 2: Schematic of a linear synchronous transit (LST) path. Figure 3: Illustration of the one-dimensional projection (ODP) technique. !, , and are and !% are the the vectors of the current point, reactants, and products, respectively. $ normalized tangent vector and perpendicular vector to the direction of the reaction, is a unit vector. The magnitude of the red line respectively. It should be noted that $ represents the “advancement of reaction ().” Figure 4: Various possible free energy surfaces between the reactant ( ) and product ( ). (a) The maximum is on the linear path between and . (b) The maximum is outside the linear path between and . (c) There are more than one reaction paths between and . (d) The true reaction path is much deviated from the linear path between and . Figure 5: Definitions of !. (a) Double proton transfer in the acetic acid dimer. ! is defined with two asymmetric coordinates ( , ). (b) Hydrolysis of the methyl diazonium ion with a two-water system. Two sets of ! and != were prepared, consisting of five and four internal coordinates, respectively. The geometries of the reactants and products in gas phase were obtained by MP2/CCT. Figure 6: One-dimensional free energy curve along the reaction coordinate () obtained with the ODP procedure. There are two symmetrical transition states (TS1 and TS2) and one intermediated state I. Figure 7: (a) Two-dimensional free energy curve generated with the reaction path, , and the perpendicular vector, )!% ). (b) Two-dimensional PMF generated with and defined in Figure 5 by the original umbrella sampling technique. The black solid line represents the reaction path. Figure 8: One-dimensional free energy curve of the methyl diazonium ion hydrolysis obtained with the ODP procedure. Key configurations are inserted along the reaction path. 18 ACS Paragon Plus Environment

Page 18 of 36

Page 19 of 36 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 dotted and solid lines represent ! and != , which have five and four internal ′ coordinates, respectively, as described in Figure 5(b). The dashed line != represents the

simulations with a smaller perpendicular restraint of 5.6 kcal/mol. Figure 9: (a) Histograms of the nucleophilic attack ( ), bond breaking ( ), and proton transfer (< ) coordinates of ! . (b) Histogram of the angle coordinate (?;: : 9: ) of ! . These values were obtained at the transition state window of 0.15 . Figure 10: Two different snapshots at the transition state window of 0.15 . (a) Tight structure with = 1.7 Å, = 2.1 Å and (b) loose structure with = 2.2 Å, = 2.1 Å. Figure 11: Histograms of the proton transfer (< ) of ! at the 0.75 window. Figure 12: Radial distribution functions (RDFs) of (a) gO1–Ow, (b) gC1–Ow, and (c) g

N1–Ow

of

the methyl diazonium ion hydrolysis for the reactant (R), transition state (TS), proton transfer (PT), and product (P), respectively, from the ! simulation. Figure 13: Two-dimensional free energy curve generated with the reaction path, , and perpendicular vector, )!% ) from the ! simulation. The black solid line indicates the reaction path. Figure 14: Integrated errors of free energy derivatives of all simulations as the function of simulation time.



Figure 1


Page 20 of 36

Page 21 of 36 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 2


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

Figure 3

22

ACS Paragon Plus Environment

Page 22 of 36

Page 23 of 36 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 4



Figure 5


Page 24 of 36

Page 25 of 36 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


Figure 6

25



Figure 7


Page 26 of 36

Page 27 of 36 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


Figure 8

27



Figure 9


Page 28 of 36

Page 29 of 36 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 10



Figure 11


Page 30 of 36

Page 31 of 36 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 12



Figure 13


Page 32 of 36

Page 33 of 36 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 14



References 1. Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D., Transition path sampling and the calculation of rate constants. Journal of Chemical Physics 1998, 108 (5), 1964-1977. 2. Schlitter, J.; Engels, M.; Krüger, P.; Jacoby, E.; Wollmer, A., Targeted Molecular Dynamics Simulation of Conformational Change-Application to the T ↔ R Transition in Insulin. Molecular Simulation 1993, 10 (2-6), 291-308. 3. Sprik, M.; Ciccotti, G., Free energy from constrained molecular dynamics. Journal of Chemical Physics 1998, 109 (18), 7737. 4. Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C., ESSENTIAL DYNAMICS OF PROTEINS. Proteins-Structure Function and Genetics 1993, 17 (4), 412-425. 5. Sugita, Y.; Okamoto, Y., Replica-exchange molecular dynamics method for protein folding. Chemical Physics Letters 1999, 314 (1-2), 141-151. 6. (a) Laio, A.; Parrinello, M., Escaping free-energy minima. Proceedings of the National Academy of Sciences 2002, 99 (20), 12562; (b) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M., Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics. The Journal of Physical Chemistry B 2006, 110 (8), 3533-3539; (c) Laio, A.; Gervasio, F. L., Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Reports on Progress in Physics 2008, 71 (12); (d) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M., Free-Energy Landscape for β Hairpin Folding from Combined Parallel Tempering and Metadynamics. Journal of the American Chemical Society 2006, 128 (41), 13435-13441. 7. Kresse, G.; Hafner, J., ABINITIO MOLECULAR-DYNAMICS FOR LIQUIDMETALS. Physical Review B 1993, 47 (1), 558-561. 8. Friesner, R. A.; Guallar, V., Ab initio quantum chemical and mixed quantum mechanics/molecular mechanics (QM/MM) methods for studying enzymatic catalysis. In Annual Review of Physical Chemistry, 2005; Vol. 56, pp 389-427. 9. Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M., Icosahedral Bond Orientational Order in Supercooled Liquids. Physical Review Letters 1981, 47 (18), 1297-1300. 10. Iannuzzi, M.; Laio, A.; Parrinello, M., Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Physical Review Letters 2003, 90 (23), 238302. 11. Branduardi, D.; Gervasio, F. L.; Parrinello, M., From A to B in free energy space. The Journal of Chemical Physics 2007, 126 (5), 054103. 12. (a) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; WiorkiewiczKuczera, J.; Yin, D.; Karplus, M., All-atom empirical potential for molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry B 1998, 102 (18), 3586-3616; (b) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J., Development and testing of the OPLS allatom force field on conformational energetics and properties of organic liquids. Journal of the American Chemical Society 1996, 118 (45), 11225-11236. 13. Yang, M.; Yang, L.; Gao, Y.; Hu, H., Combine umbrella sampling with integrated tempering method for efficient and accurate calculation of free energy changes of complex energy surface. The Journal of Chemical Physics 2014, 141 (4), 044108. 14. Uddin, N.; Choi, T. H.; Choi, C. H., Direct Absolute pKa Predictions and Proton Transfer Mechanisms of Small Molecules in Aqueous Solution by QM/MM-MD. The Journal of Physical Chemistry B 2013, 117 (20), 6269-6275. 34 ACS Paragon Plus Environment

Page 34 of 36

Page 35 of 36 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


15. (a) Peng, C. Y.; Schlegel, H. B., COMBINING SYNCHRONOUS TRANSIT AND QUASI-NEWTON METHODS TO FIND TRANSITION-STATES. Israel Journal of Chemistry 1993, 33 (4), 449-454; (b) Halgren, T. A.; Lipscomb, W. N., The synchronoustransit method for determining reaction pathways and locating molecular transition states. Chemical Physics Letters 1977, 49 (2), 225-232. 16. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. Journal of Computational Chemistry 1983, 4 (2), 187-217. 17. Grossfield, A. "WHAM: the Weighted Histogram Analysis Method", Version 2.0.9, http://membrane.urmc.rochester.edu/content/wham 2012. 18. Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; Ishikawa, T.; Nakano, T., How Does an SN2 Reaction Take Place in Solution? Full Ab Initio MD Simulations for the Hydrolysis of the Methyl Diazonium Ion. Journal of the American Chemical Society 2008, 130 (8), 2396-2397. 19. Mones, L.; Csányi, G., Topologically Invariant Reaction Coordinates for Simulating Multistate Chemical Reactions. The Journal of Physical Chemistry B 2012, 116 (51), 1487614885. 20. Mones, L.; Kulhánek, P.; Simon, I.; Laio, A.; Fuxreiter, M., The Energy Gap as a Universal Reaction Coordinate for the Simulation of Chemical Reactions. The Journal of Physical Chemistry B 2009, 113 (22), 7867-7873.



TOC

(kcal/mol)

B �

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

A ξ


Page 36 of 36

One-Dimensional Projection of Collective Variables for Effective

Recommend Documents