Extended Adaptive Biasing Force Algorithm. An On-the-Fly

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

The extended adaptive biasing force algorithm. An on-thefly implementation for accurate free-energy calculations Haohao Fu, Xueguang Shao, Christophe Chipot, and Wensheng Cai J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00447 • Publication Date (Web): 11 Jul 2016 Downloaded from http://pubs.acs.org on July 15, 2016

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 33

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

The extended adaptive biasing force algorithm. An on-the-fly implementation for accurate free-energy calculations Haohao Fu,† Xueguang Shao,†, ‡, # Christophe Chipot,*,$,§,¶ and Wensheng Cai*,†,‡

†

Research Center for Analytical Sciences, College of Chemistry, Tianjin Key Laboratory of

Biosensing and Molecular Recognition, Nankai University, Tianjin 300071, China ‡

Collaborative Innovation Center of Chemical Science and Engineering (Tianjin), Tianjin

300071, China #

State Key Laboratory of Medicinal Chemical Biology (Nankai University), Tianjin 300071,

China $

Laboratoire International Associé Centre National de la Recherche Scientifique et University of

Illinois at Urbana-Champaign, Unité Mixte de Recherche No. 7565, Université de Lorraine, B.P. 70239, 54506 Vandœuvre-lès-Nancy cedex, France §

Theoretical and Computational Biophysics Group, Beckman Institute, University of Illinois at

Urbana-Champaign, Urbana, Illinois 61801 ¶

Department of Physics, University of Illinois at Urbana−Champaign, 1110 West Green Street,

Urbana, Illinois 61801, United States

ACS Paragon Plus Environment 1


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

ABSTRACT

Proper use of the adaptive biasing force (ABF) algorithm in free-energy calculations needs certain prerequisites to be met, namely, that the Jacobian for the metric transformation and its first derivative be available and the coarse variables be independent and fully decoupled from any holonomic constraint or geometric restraint, thereby limiting singularly the field of application of the approach. The extended ABF (eABF) algorithm circumvents these intrinsic limitations by applying the time-dependent bias onto a fictitious particle coupled to the coarse variable of interest by means of a stiff spring. However, with the current implementation of eABF in the popular molecular dynamics engine NAMD, a trajectory-based post-treatment is necessary to derive the underlying free-energy change. Usually, such a post-hoc analysis leads to a decrease in the reliability of the free-energy estimates due to the inevitable loss of information, as well as to a drop in efficiency, which stems from substantial read-write accesses to file systems. We have developed a user-friendly, on-the-fly code for performing eABF simulations within NAMD. In the present contribution, this code is probed in eight illustrative examples. The performance of the algorithm is compared with traditional ABF, on the one hand, and the original eABF implementation combined with a post-hoc analysis, on the other hand. Our results indicate that the on-the-fly eABF algorithm (i) supplies the correct free-energy landscape in those critical cases where the coarse variables at play are coupled to either each other, or to geometric restraints or holonomic constraints, (ii) greatly improves the reliability of the freeenergy change, compared to the outcome of a post-hoc analysis, and (iii) represents a negligible


Page 2 of 33

Page 3 of 33

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


additional computational effort compared to regular ABF. Moreover, in the proposed implementation, guidelines for choosing two parameters of the eABF algorithm, namely the stiffness of the spring and the mass of the fictitious particles, are proposed. The present on-thefly eABF implementation can be viewed as the second generation of the ABF algorithm, expected to be widely utilized in the theoretical investigation of recognition and association phenomena relevant to physics, chemistry and biology.



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 4 of 33

INTRODUCTION Molecular dynamics (MD) has been broadly used for studying proteins,1 nucleic acids,2 membranes,3 organic polymers,4 and supramolecules5 at the atomic level. A host of processes, notably of biological relevance are, however, not amenable to classical MD owing to the timescales that they span. Considering the inherent complexity of the modeled objects at play, alternate avenues, therefore, need to be examined to sample the underlying, characteristically slow degrees of freedom. Two groups of methods, namely generalized-ensemble schemes6 and importance-sampling schemes,7 emerged as inescapable tools towards achieving this goal. In generalized-ensemble strategies, uniform sampling of configurational space is expected to be achieved by introducing an auxiliary probability weight factor. In importance-sampling numerical schemes, some coarse variables are believed to be more important than others and sampling along these variables is artificially encouraged, notably by means of biasing potentials. Umbrella sampling,8-9 adaptive biasing force (ABF),10-12 and metadynamics13-15 pertain to the latter class of approaches and are commonly used to explore complex free-energy landscapes. The classical ABF algorithm, developed by Darve et al,10-11 is based on thermodynamic integration,16 i.e., (1) where

is the free-energy difference,

transition coordinate, , in each bin , of the Jacobian,

is the average force acting along the chosen is the potential energy function,

is the Boltzmann constant and

is the determinant

is the temperature.

Classical ABF has been used successfully to reveal the mechanisms that underlie many


Page 5 of 33

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


processes of biological, physical and chemical relevance, like the extension of mechanical proteins,17 membrane permeation,18 as well as recognition-association phenomena,19 owing to its rigorous mathematical foundations and remarkable convergence properties.20 To perform a classical ABF simulation, the following three requirements should be met: (i) Availability of the derivative of the Jacobian. As shown in eq. (1), calculation of the biasing force

requires the analytical expression of the derivative of the Jacobian,

,

which involves the second derivative of the coarse variable and can sometimes take a mathematically complicated form.21 This requirement can be circumvented in those cases where the coarse variables can be described by a vector-field formalism,22 instead of the standard expression of eq. (1). In the Colvars module23 of the current version of NAMD (2.11),24 only the following types of coarse-variable components can be used in ABF calculations — Euclidean distance between two groups, projection of a distance vector onto an axis or a plane, angle between three groups, torsional angle between four groups, radius of gyration of a group of atoms, root mean square deviation (RMSD) from a set of reference coordinates and projection of the atomic coordinates on a vector. (ii) Independence of the coarse variables. In a multidimensional ABF calculation, each coarse variable should be mutually orthogonal. Similarly, when a coarse variable consists of a combination of several components, these components should also be mutually orthogonal. (iii) Orthogonality of the coarse variables and constraints or restraints. If there are geometric restraints or holonomic constraints in the assay, they should be fully decoupled from the coarse variables. Should holonomic constraints be involved in the model reaction coordinate, the



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 6 of 33

incorrect free-energy change stems from the contamination of the measured thermodynamic force by the constraint forces. When geometric restraints are coupled to the coarse variables, regular ABF is not able to discriminate the origin of the forces at play, resulting in the application of an erroneous time-dependent bias. Said differently, ABF is not biasing the correct free-energy hyperplane. The aforementioned three prerequisites markedly limit the practical application of the ABF method, especially in characterizing processes of complex systems. An alternate approach proposed by Lelièvre et al.7, 25 and Zheng et al.26-27 may address this issue.23, 28 This scheme is implemented by applying a biasing force to an added fictitious particle which couples the coarse variable by a stiff spring, and thus an extended potential is introduced to the system, (2) where

is a fictitious degree of freedom,

is the fictitious particle, and

> 0 is the force

constant of the introduced spring. This approach is called extended ABF, or eABF. On account of the orthogonality of the extended degree of freedom produced by the stiff spring, the biasing force becomes nearly always independent. In addition, compared with classical ABF, the measurement of the biasing force in eABF is greatly simplified because it can be directly obtained by Hooke's law. The prerequisites of classical ABF, therefore, are expected to be removed.23, 29 However, the free-energy change along ,

, is only an approximation of the

original PMF.26 To calculate the PMF rigorously, therefore, a method to perform deconvolution of the harmonic-spring contribution was proposed for the first time by Zheng and Yang, using the following equation,26


Page 7 of 33

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) where

is the number of samples along

collected from the

-restrained ensemble.

In practice, deconvolution of the harmonic-spring contribution is carried out by means of a post-hoc analysis of the output data generated by the Colvars free-energy calculation engine.30 To perform the post-treatment, recording of the trajectory of the value of the coarse variables is needed. This strategy is, therefore, prone to result in a loss of information because such data usually cannot be stored at every time step. It is also likely to affect the reliability of the freeenergy estimated. In addition, the post-treatment process itself has proven under certain circumstances to be time-consuming. To address these shortcomings, we have developed a code to perform on-the-fly eABF calculations, resting on the strategy proposed in orthogonal space tempering, wherein eABF is employed as the recursion kernel.26, 31 The computational cost and performance of the new algorithm were probed in eight examples, and the results were compared with those obtained by traditional ABF10-11 and eABF with post-hoc analysis,30 on the other hand. We primarily focused on the performance of on-the-fly eABF in those particular cases where the coarse variables are not independent, and on the effect of the different data-processing methods on the reliability of free-energy profiles. Moreover, the choice of the spring constant and fictitious mass appearing in the eABF algorithm was investigated, leading to the suggestion of practical guidelines for eABF calculations.



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

METHODS Implementation Details A code to carry out the deconvolution of the harmonic-spring contribution in eABF calculation, summarized in eq. (3), is proposed. This code is called by the popular molecular dynamics program NAMD24 at every time step to perform on-the-fly eABF calculations. Communication between this code and the free-energy calculation engine is implemented through the scripted interface of the Colvars module.23 The computationally intensive component of the code is written in C++. The user interface has been devised to reduce human intervention as much as possible. Apart from the inclusion of a number of specific parameters in the NAMD configuration file, no additional input is needed (see the Supporting Information). By design, onthe-fly eABF calculation is made completely transparent for the end-user.

Free-energy Calculations Six molecular systems were selected to probe our on-the-fly implementation, as shown in Figure 1, where the chosen model reaction coordinates are highlighted. Two other test cases are presented in Figure S1-S3 of the Supporting Information. A detailed description of the molecular assemblies and the simulation times is provided in Table S1 of the Supporting Information. It is worth noting that the coarse variables depicted in Figure 1A and 1D are coupled, and that appearing in Figure 1B is not orthogonal to the holonomic constraints used to freeze the geometry of the water molecules. Classical ABF10-12 was employed to determine the PMF associated to the first two cases in order to illustrate the inherent limitations of the algorithm.


Page 8 of 33

Page 9 of 33

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


Interestingly enough, in practice, violating the last two aforementioned requirements of regular ABF will not lead to runtime errors of the MD engine, but result in an incorrect free-energy landscape. For comparison, the eABF algorithm25-26 combined with a post-treatment strategy30 was utilized to obtain the free-energy profile characterizing the conformational equilibrium of Figure 1F. A multiple-walker strategy was used to enhance ergodic sampling in all the freeenergy calculations characterizing the isomerization of NANMA.32-35 Unless specified otherwise, the standard deviation between the fictitious particle and the connected atomic group was equal to the width of the bins in which

and

accrue, that is the bin width, and the oscillation period

of the fictitious particle was 200 fs. A detailed discussion of the choice of these two parameters is provided below.

Figure 1. Molecular assays and corresponding transition coordinates used in this study. (A) Conformational change of cyclohexane. The transition coordinate is defined as the average of three dihedral angles. (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

Page 10 of 33

Disassociation of a water dimer. The transition coordinate is defined as the intermolecular O···H distance. (C) Unfolding of deca-alanine. The transition coordinate is defined as the distance between two carbon atoms (Cα) at the termini of the peptide chain. (D) Interaction of an acetate ion with a sodium ion. The transition coordinate is defined as the sum of two distances. (E) Conformational change in n-butane. The transition coordinate is defined as the C-C-C-C torsional angle. (F) Isomerization of N-acetyl-N’-methyl-L–alanylamide (NANMA). The transition coordinate consists of two backbone dihedral angles. The cyan box is a stylized representation of the aqueous environment.

Simulation Details All the atomistic MD simulations presented herein were performed employing the parallel, scalable program NAMD 2.11.24 The TIP3P model,36 the CHARMM force field37-38 and the CHARMM general force field (CGenFF)39 were used to describe water, the peptides and the small organic molecules, respectively. For molecules examined in vacuum, a time step of 0.5 fs was utilized to integrate the equations of motion. The temperature was maintained at 300 K employing Langevin dynamics. Specifically, the SETTLE algorithm40 was used to constrain the covalent bonds to their equilibrium length only when simulating the assay of Figure 1B to show the effect of holonomic constraints on free-energy change estimated with regular ABF. For those simulations in a condensed phase, the r-RESPA multiple time-step algorithm41 was employed to integrate the equations of motion with a time step of 1 and 2 fs for short- and long-range interactions, respectively. Covalent bonds involving hydrogen atoms were constrained to their equilibrium length by means of the SHAKE/RATTLE42-43 and SETTLE algorithms40 for the (bio)organic and water molecules, respectively. The temperature and the pressure were maintained at 300 K and 1 atm, respectively, using Langevin dynamics and the Langevin piston method.44 Long-range electrostatic forces were taken into account by means of the particle mesh


Page 11 of 33

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


Ewald algorithm.45 A 12-Å cutoff was introduced to truncate van der Waals and short-range Coulombic interactions. Periodic boundary conditions (PBCs) were applied in the three directions of Cartesian space. Visualization and analyses of the MD trajectories were performed with the VMD program.46

RESULTS AND DISCUSSTION Applicability, Reliability and Efficiency of On-the-fly eABF Applicability of the eABF algorithm. The two cases depicted in Figure 1A and 1B were used to compare the applicability of eABF and classical ABF. The free-energy profile delineating the chair-chair conversion of cyclohexane was determined with these two methods. Figure 2A reveals an obvious difference between the two resulting PMFs. Periodicity considerations in cyclohexane would suggest that the two chair conformers at ±53° correspond to the same free energy.47-49 Yet, plain ABF yields a free-energy difference of about 13.8 kcal/mol. This counterintuitive result stems from the fact that the three components of the transition coordinate are not independent, thereby violating the second requirement of the algorithm. A similarly erroneous result can be found in the separation of two water molecules. As highlighted in Figure 2B, the result supplied by regular ABF incorrectly suggests that the binding of the two water molecules is thermodynamically unfavorable, at variance with both experimental and theoretical findings.50-51 This unexpected result can be explained by the coupling of the biasing force in ABF and the constraint force from the SETTLE algorithm. In stark contrast, Figure 2A indicates that the eABF algorithm predicts an energetic



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 12 of 33

equivalence of the two chair conformers (±53°). Figure 2A also shows that the free energy of the twist-boat conformer lies approximately 5.3 kcal/mol above that of the global minimum, i.e., the chair motif, and that a free-energy barrier of approximately 10.8 kcal/mol ought to be overcome for cyclohexane to isomerize between these two conformational states, in good agreement with the experimental values of about 5.5 and 10.4 kcal/mol, respectively.47-48 These results are consistent with previous theoretical findings of 5.6 and 9.0 kcal/mol, respectively.49

Figure 2. Free-energy landscapes for (A) the isomerization of cyclohexane and (B) the unbinding of two water molecules in the gas phase, derived from ABF and eABF simulations. The Error bars correspond to the


Page 13 of 33

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


standard error inferred from ten independent 2-ns runs. The error bars for the red curve are too small to be visualized.

Disassociation of the water dimer has been extensively studied, using different force fields and model reaction coordinates.52-54 In these theoretical investigations, the free-energy difference between the free and binding states is reported to range between 2.2 and 3.6 kcal/mol, depending on the water model used in the simulation. The free-energy difference inferred from the PMF of Figure 2B amounts to 3.1 kcal/mol, congruent with the reported values. Two additional free-energy calculations, characterizing, on the one hand, the variation of the coordination number of a potassium ion in an aqueous solution and, on the other hand, the isomerization of deoxycytidine, were carried out. In the former case, which corresponds to a discrete coarse variable, the Jacobian and its first derivative are not available, whereas, in the latter case, the two coarse variables are coupled. In other words, neither assay meets the prerequisites of the classical ABF algorithm — which, therefore, is unlikely to supply the correct answer. The PMFs obtained from on-the-fly eABF simulations, however, coincide with previous theoretical studies, thereby, suggesting that eABF can be safely utilized in the present cases (see Figure S1-S3 of the Supporting Information). Figure 2 and S1-S3 cogently illuminate that eABF is perfectly adapted to those examples wherein the three prerequisites for ABF simulations are not met, hence, greatly expanding the range of applications amenable to ABF-based free-energy calculations. Reliability of on-the-fly eABF calculations. The isomerization of N-acetyl-N’-methyl-L– alanylamide (NANMA) was investigated by the eABF algorithm using both the post-hoc



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

analysis and the on-the-fly scheme in order to explore the effect of data processing on the reliability the free-energy estimates. It has been demonstrated in a previous study that isomerization of such a short peptide can be properly studied by means of the regular ABF algorithm.29 Use will, therefore, be made of the free-energy landscape mapped from a long ABF calculation as a reference here. Figure

3A

shows

the

two-dimensional

free-energy

surface

characterizing

the

conformational equilibrium of NANMA obtained from a 120-ns ABF simulation. The freeenergy difference of the right-handed α-helical, αR, and the β-strand, β, conformations is calculated using the following equation, (4) where

is the probability distribution function associated with

the transition coordinate

.

and

, the PMF along

are the probabilities that system is in αR and β state,

respectively. Integration over the basins delineating the αR and the β conformations yields a free-energy difference of about 0.2 kcal/mol, in favor of the former state, congruent with previous computer simulations with the same potential energy function,55-57 suggestive that our ABF calculation is likely to be adequately converged. The free-energy surface derived from the eABF simulation combined with a post-hoc analysis is depicted in Figure 3B. That obtained using the on-the-fly scheme is shown in Figure 3C. The root-mean-square (RMSD) between the free energy of Figure 3B and 3C and that of Figure 3A amounts to 0.4 and 0.1 kcal/mol, respectively, confirming the


Page 15 of 33

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


greater reliability of the on-the-fly eABF calculation. The difference in reliability between the simulation followed by a post-hoc analysis and that wherein data processing proceeds on the fly stems from a loss of information in the former. Traditionally, in eABF calculations, the trajectory of

and

appearing in eq. (3) is recorded in a text file. Taking into account, on the one hand,

the capacity and read/write speed of hard drives and, on the other hand, the computational effort involved in the post treatment (see Table 1), one usually saves the relevant information every hundreds or thousands of MD steps, thus leading to significant data loss. In the example proposed herein, the time interval for writing the trajectory was 200 MD steps. The average number of samples per bin in the on-the-fly eABF calculation is 23,000, while that in the traditional calculation followed by a post-hoc analysis only amounted to 115.

Figure 3. Two-dimensional free-energy landscape underlying the isomerization of NANMA, obtained from a 120-ns (A) classical ABF calculation, utilized here as a reference, (B) eABF calculation followed by a posthoc analysis, and (C) on-the-fly eABF calculation. The basins delineating the αR and the β conformations are defined in Figure 3A.

Table 1. Comparison of on-the-fly eABF and eABF with post-processing of the accrued data. A 120 ns freeenergy calculation was performed for each method. Method Time interval for Computational a data writing cost (MD steps) On-the-fly eABF 8 hours eABF with post200 8 hours

Size of produced files d

20 MB 60 MB


c

Post-treatment b Time

RMSD (kcal/mol)

4 hours

0.1 0.4


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 16 of 33

hoc analysis eABF with post1 10 hours 12 GB 4 days hoc analysis d Classical ABF 8 hours 2 MB a Computer Configuration: Intel Xeon E5-2670v3*2/64GB memory/Geforce GTX 980Ti *2. b

0.1 0

We contend that the current post-treatment suffers from low efficiency also because of the slow AWK

scripting language, especially for two-dimensional free-energy calculations. c

RMSD of the PMF by eABF with respect to that by regular ABF.

d

Size of the restart files in text format.

The aforementioned results unambiguously indicate that on account of the loss of information, the post-treatment strategy in eABF calculations greatly affects the reliability of the computed free-energy estimates. Moreover, considering that the test case of Figure 3 is only a toy model, the observed drop in the reliability may become a far more acute issue for complex molecular assemblies. An on-the-fly strategy is necessary in the calculation of a more complex system such as a protein-ligand assembly. Efficiency of the on-the-fly eABF implementation. The efficiency of the on-the-fly eABF implementation can be regarded from two vantage points, namely (i) the additional computational effort induced by the code, and (ii) the convergence rate of the algorithm compared with classical ABF. As shown in Table 1, no perceptible difference is observed between the new eABF implementation and the traditional ABF scheme insofar as computational cost is concerned, which is because evaluation of eq. (3) is computationally cheaper than force evaluation in MD. Overall, the computational cost for an on-the-fly eABF simulation is not more than that of regular ABF. To compare the convergence rate of on-the-fly eABF and classical ABF, the time evolution of the RMSD of the gradients was monitored in the free-energy calculation of NANMA — i.e., the RMSD between

and


, that is the RMSD

Page 17 of 33

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


between the average force at time t and at the beginning of the trajectory. As shown in Figure 4, the eABF calculation converged within 10 ns, while classical ABF simulation thrice as much sampling to reach the same level of convergence. A similar trend was observed in three parallel runs, indicating that the convergence of on-the-fly eABF can be, at least in some favorable cases, faster than that of traditional ABF. One may rationalize this observation as follows. Since the fictitious particle is only affected by a hard spring, there is no hidden free-energy barrier in the space orthogonal to the direction of the transition coordinate. Hence, the fictitious particle diffuses along the latter direction much easier than the complex assembly, the diffusion of which is often thwarted by hidden barriers in traditional ABF. The real system is, therefore, often coerced by the hard spring to move along the transition coordinate, leading to a uniform sampling and a high convergence rate.

Figure 4. Time evolution of the root-mean-square deviation (RMSD) over the gradients of the free-energy surfaces depicted in Figure 3A and 3C. Inset: zoom-in of the main graph.



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 18 of 33

Choice of the Parameters of eABF Calculations As shown in eq. (2), stiff springs and fictitious particles are added to the system when performing an eABF free-energy calculation. Obviously, two additional parameters, i.e., the stiffness factor of the springs and the mass of the particles, should be defined as a preamble to any free-energy calculations. These two parameters are calculated by eq. (5) and (6),58 (5) (6) where

is the standard deviation between the fictitious particle and the collective variable to

which it is connected. The difference between the trajectory of the collective variable and that of the extended system is Gaussian distributed. In the Colvars module of NAMD,23

is the oscillation period of the fictitious particle.

is controlled by the keyword extendedFluctuation and

is

defined by keyword extendedTimeConstant, which are both set by the users. To investigate systematically the effects of these parameters on the calculated free-energy change, the five molecular assays depicted in Figure 1A-E were simulated repeatedly using different parameters. As an example, the PMFs characterizing the conformational change of cyclohexane determined with different combinations of parameters are shown in Figure 5 (the results for the four other assays are relegated in the Supporting Information, Figure S4-S7). eABF yields reasonable results only when extendedFluctuation ranges from 0.5° to 2°, as highlighted in Figure 5A. When extendedFluctuation is too large or too small, the free-energy profiles deviate from that of the reference simulation. Interestingly enough, we found in all test cases that the result of the free-energy calculations is always reasonable when the chosen


Page 19 of 33

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


extendedFluctuation coincides with the bin width, as shown in Figure S4-S7 of the Supporting Information. Keeping in mind that the fictitious particle and the molecular system are connected by a spring, the effect of extendedFluctuation on the calculated free-energy landscape can be explained as follows. A small extendedFluctuation implies a large spring constant, corresponding to a hard spring. A small oscillation of the fictitious particle would, therefore, cause a huge variation of the collected biasing force, leading to strong fluctuations of the average force, which can be seen from the black curve in Figure 5A. Conversely, a large extendedFluctuation corresponds to a soft spring. Under these premises, the movement of the system couples weakly to that of the fictitious particle, causing the real PMF (see the yellow curve in Figure 5A).


to depart significantly from


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 20 of 33

Figure 5. Free-energy profiles characterizing the isomerization of cyclohexane, using different values for the eABF parameters (A) extendedFluctuation and (B) extendedTimeConstant.

Compared to that of extendedFluctuation, the effect of extendedTimeConstant on the freeenergy landscape is relatively minute. As confirmed by Figure 5B and S4-S7, this parameter has only a marginal influence on the final free-energy profile when it is much larger than the MD time step. The default value in the Colvars module23 of 200 fs appears to be suitable for all the molecular assays examined here.


Page 21 of 33

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


Apart from extendedFluctuation and extendedTimeConstant, the parameters associated to the ABF algorithm, such as the bin width, the force constant of the potentials confining sampling within a window, the threshold number of samples per bin prior to application of the timedependent bias, are also inherited by the eABF method. A detailed discussion of the choice of appropriate values for these parameters can be found in a previous work.59

CONCLUSION The on-the-fly eABF strategy put forth in the present contribution is envisioned to expand in a significant manner the scope of applications amenable to ABF-based free-energy calculations, notably in those cases where the criteria on the model reaction coordinate that underlie regular ABF simulations are not met. We have shown in the test examples proposed herein that in those cases where the derivative of the Jacobian is not available and/or the coarse variables are not independent, or are coupled with a holonomic constraint or a geometric restraint, eABF will yield the correct PMF. Moreover, under certain circumstances and when both methods are suitable, eABF may sample the transition coordinate more efficiently than classical ABF. The code presented herein, designed for performing on-the-fly eABF calculations, also offers at negligible additional cost free-energy estimates of markedly increased reliability by obviating the need for a post-hoc analysis generally carried out on an incomplete dataset and, thus, eliminating possible losses of information can also prevent from the loss of information in the traditional eABF method with post-hoc analysis, at negligible additional cost. In many



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 22 of 33

respects, the on-the-fly eABF strategy detailed in this contribution can be viewed as the second generation of the ABF algorithm for free-energy calculations. To improve the reliability of the PMFs determined using the eABF scheme, we suggest that the additional parameters of the algorithm be chosen as follows, namely, that extendedFluctuation be consistent with the bin width and the extendedTimeConstant be the default value of the Colvars module. We observe that whereas the calculated free-energy profiles are sensitive to the value of the former, they appear to be by and large insensitive to that of the latter. Free-energy calculations aimed at unveiling with exquisite atomic detail the physical principles that underlie complex chemical and biological processes always constitute a significant challenge for theory. Developing novel algorithms or improving existing methods is of paramount importance. In the present work, we have implemented an on-the-fly eABF strategy to virtually eliminate the intrinsic shortcomings of ABF-based free-energy calculations. Notwithstanding the advances made possible by this development, a number of practical issues remain to be addressed. In particular, the definition of the transition coordinate — and, hence, the choice of the coarse variables at play, may be too rudimentary to characterize with the desired efficiency and accuracy processes of physical, chemical and biological relevance.60-62 A variety of methods were put forth to address the problem at hand, which can be roughly classified into two groups, namely (i) resort to a small number of coarse variables to model the reaction coordinate, while improving sampling in the orthogonal space, such as orthogonal space tempering,26,

31, 63-64

or (ii) take advantage of multidimensional free-energy calculations to

describe complex transformations, yet focusing sampling along energetically favorable pathways,


Page 23 of 33

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


as would be the case in self-learning adaptive umbrella sampling65-66 and the so-called on-thepath random walk.67 Another approach that pertains to the latter group is a variant of ABF, coined generalized ABF (gABF), which has proven to possess the ability to find approximate minimum-action pathways for complex transformations based on rather short simulations.68 Combination of eABF and gABF could constitute yet another route to address sampling of freeenergy landscapes defined by a limited number of coarse variables, albeit assessment of this alternative would require further exploration, which goes beyond the present contribution.

ASSOCIATED CONTENT Supporting Information The Supporting Information is available free of charge on the ACS Publications website. Details of the composition and the size of the molecular assemblies. Free-energy calculations charactering the variation of the coordination number of a potassium ion in aqueous solution and the pesudorotation of five-membered rings. The PMFs characterizing the process of Figure 1B1E, adopting different extendedFluctuation and extendedTimeConstant. On-the-fly eABF code with an example. Comparison of (e)ABF, umbrella sampling and metadynamics.



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 24 of 33

AUTHOR INFORMATION Corresponding Authors *

E-mail: [email protected] (W.C.).

*

E-mail: [email protected] (C.C.).

Notes The authors declare no competing financial interest.

ACKNOWLEDGMENTS The authors gratefully acknowledge Drs. Wei Yang, Tony Lelièvre and Jérôme Hénin for helpful discussions and insightful comments. The study was supported by National Natural Science Foundation of China (No. 21373117), MOE Innovation Team (IRT13022) of China. The GENCI and CINES, Montpellier, France, and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) of China, are gratefully acknowledged for provision of generous amounts of CPU time. The Cai Yuanpei program is also appreciatively acknowledged for its support of the international collaboration between the research groups of C.C. and W.C.

REFERENCES

1.

Karplus, M.; Kuriyan, J. Molecular dynamics and protein function. Proc. Natl. Acad. Sci.


Page 25 of 33

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


U.S.A. 2005, 102 (19), 6679-6685. 2.

Šponer, J.; Banáš, P.; Jurečka, P.; Zgarbová, M.; Kührová, P.; Havrila, M.; Krepl, M.; Stadlbauer, P.; Otyepka, M. Molecular dynamics simulations of nucleic acids. From tetranucleotides to the ribosome. J. Phys. Chem. Lett. 2014, 5 (10), 1771-1782.

3.

Milton Jr, H.; Stiles, C. D. Molecular dynamics in biological membranes. Springer Science & Business Media: 2012.

4.

Barrat, J. L.; Baschnagel, J.; Lyulin, A. Molecular dynamics simulations of glassy polymers. Soft Matter 2010, 6 (15), 3430-3446.

5.

Wipff, G. Computational approaches in supramolecular chemistry. Springer Science & Business Media: 2012; Vol. 426.

6.

Mitsutake, A.; Sugita, Y.; Okamoto, Y. Generalized-ensemble algorithms for molecular simulations of biopolymers. Biopolymers 2001, 60 (2), 96-123.

7.

Lelièvre, T.; Stoltz, G.; Rousset, M. Free energy computations: A mathematical perspective. World Scientific: 2010.

8.

Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23 (2), 187-199.

9.

Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13 (8), 1011-1021.

10. Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115 (20), 9169-9183.



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 26 of 33

11. Darve, E.; Wilson, M. A.; Pohorille, A. Calculating free energies using a scaled-force molecular dynamics algorithm. Mol. Simulat. 2002, 28 (1-2), 113-144. 12. Hénin, J.; Chipot, C. Overcoming free energy barriers using unconstrained molecular dynamics simulations. J. Chem. Phys. 2004, 121 (7), 2904-2914. 13. Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (20), 12562-12566. 14. Huber, T.; Torda, A. E.; van Gunsteren, W. F. Local elevation: a method for improving the searching properties of molecular dynamics simulation. J. Comput. Aid. Mol. Des. 1994, 8 (6), 695-708. 15. Grubmüller, H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys. Rev. E 1995, 52 (3), 2893. 16. Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3 (5), 300313. 17. Lee, E. H.; Hsin, J.; Mayans, O.; Schulten, K. Secondary and tertiary structure elasticity of titin Z1Z2 and a titin chain model. Biophys. J. 2007, 93 (5), 1719-1735. 18. Wei, C.; Pohorille, A. Permeation of membranes by ribose and its diastereomers. J. Am. Chem. Soc. 2009, 131 (29), 10237-10245. 19. Rodriguez, J.; Semino, R. o.; Laria, D. Building up nanotubes: Docking of “Janus” cyclodextrins in solution. J. Phys. Chem. B 2009, 113 (5), 1241-1244. 20. Lelièvre, T.; Minoukadeh, K. Long-time convergence of an Adaptive Biasing Force method: the bi-channel case. Arch. Ration. Mech. An. 2011, 202 (1), 1-34.


Page 27 of 33

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


21. Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128 (14), 144120. 22. Den Otter, W. Thermodynamic integration of the free energy along a reaction coordinate in Cartesian coordinates. J. Chem. Phys. 2000, 112 (17), 7283-7292. 23. Fiorin, G.; Klein, M. L.; Hénin, J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111 (22-23), 3345-3362. 24. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781-1802. 25. Lelièvre, T.; Rousset, M.; Stoltz, G. Computation of free energy profiles with parallel adaptive dynamics. J. Chem. Phys. 2007, 126 (13), 134111. 26. Zheng, L.; Yang, W. Practically efficient and robust free energy calculations: doubleintegration orthogonal space tempering. J. Chem. Theory Comput. 2012, 8 (3), 810-823. 27. Zheng, L.; Chen, M.; Yang, W. Random walk in orthogonal space to achieve efficient freeenergy simulation of complex systems. Proc. Natl. Acad. Sci. U.S.A. 2008, 105 (51), 2022720232. 28. Gumbart, J. C.; Roux, B.; Chipot, C. Standard binding free energies from computer simulations: What is the best strategy? J. Chem. Theory Comput. 2012, 9 (1), 794-802. 29. Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. J. Chem. Theory Comput. 2010, 6 (1), 35-47.



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 28 of 33

30. Gumbart, J. C.; Roux, B.; Chipot, C. Protein: ligand standard binding free energies: A tutorial for alchemical and geometrical transformations. http://www.ks.uiuc.edu/Training/ Tutorials/namd/PLB/tutorial-protein-ligand.pdf (accessed April 22, 2016). 31. Lv, C.; Aitchison, E. W.; Wu, D.; Zheng, L.; Cheng, X.; Yang, W. Comparative exploration of hydrogen sulfide and water transmembrane free energy surfaces via orthogonal space tempering free energy sampling. J. Comput. Chem. 2015, 37 (6), 567-574. 32. Minoukadeh, K.; Chipot, C.; Lelievre, T. Potential of mean force calculations: A multiplewalker adaptive biasing force approach. J. Chem. Theory Comput. 2010, 6 (4), 1008-1017. 33. Comer, J.; Phillips, J. C.; Schulten, K.; Chipot, C. Multiple-replica strategies for free-energy calculations in NAMD: Multiple-Walker adaptive biasing force and Walker selection rules. J. Chem. Theory Comput. 2014, 10 (12), 5276-5285. 34. Comer, J.; Gumbart, J. C.; Hénin, J.; Lelièvre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: Everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2014, 119 (3), 1129-1151. 35. Comer, J.; Roux, B.; Chipot, C. Achieving ergodic sampling using replica-exchange freeenergy calculations. Mol. Simulat. 2014, 40 (1-3), 218-228. 36. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926935. 37. MacKerell Jr, A. D.; Bashford, D.; Bellott, M.; Dunbrack Jr, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S. All-atom empirical potential for molecular


Page 29 of 33

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


modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102 (18), 3586-3616. 38. MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126 (3), 698-699. 39. Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I. CHARMM general force field: A force field for druglike molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31 (4), 671-690. 40. Miyamoto, S.; Kollman, P. A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13 (8), 952-962. 41. Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 1992, 97 (3), 1990-2001. 42. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23 (3), 327-341. 43. Andersen, H. C. RATTLE: A “Velocity” version of the SHAKE algorithm for molecular dynamics calculations. J. Comput. Phys. 1983, 52 (1), 24-34. 44. Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: the Langevin piston method. J. Chem. Phys. 1995, 103 (11), 4613-4621. 45. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N⋅ log (N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98 (12), 10089-10092. 46. Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph.



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 30 of 33

1996, 14 (1), 33-38. 47. Squillacote, M.; Sheridan, R.; Chapman, O.; Anet, F. Spectroscopic detection of the twistboat conformation of cyclohexane. Direct measurement of the free energy difference between the chair and the twist-boat. J. Am. Chem. Soc. 1975, 97 (11), 3244-3246. 48. Ross, B. D.; True, N. S. NMR spectroscopy of cyclohexane. Gas-phase conformational kinetics. J. Am. Chem. Soc. 1983, 105 (15), 4871-4875. 49. Bucher, D.; Pierce, L. C.; McCammon, J. A.; Markwick, P. R. On the use of accelerated molecular dynamics to enhance configurational sampling in ab initio simulations. J. Chem. Theory Comput. 2011, 7 (4), 890-897. 50. Curtiss, L.; Frurip, D.; Blander, M. Studies of molecular association in H2O and D2O vapors by measurement of thermal conductivity. J. Chem. Phys. 1979, 71 (6), 2703-2711. 51. Kim, K. S.; Mhin, B. J.; Choi, U. S.; Lee, K. Ab initio studies of the water dimer using large basis sets: The structure and thermodynamic energies. J. Chem. Phys. 1992, 97 (9), 66496662. 52. Liu, J.; Yang, L.; Doren, D. J. Free energy perturbation and dynamical nucleation study of water dimer and trimer through TIP5P water model. Chem. Phys. Lett. 2006, 417 (1), 63-71. 53. 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 (2), 773-777. 54. Tritzant-Martinez, Y.; Zeng, T.; Broom, A.; Meiering, E.; Le Roy, R. J.; Roy, P. N. On the analytical representation of free energy profiles with a Morse/long-range model: Application to the water dimer. J. Chem. Phys. 2013, 138 (23), 234103.


Page 31 of 33

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


55. Tobias, D. J.; Brooks III, C. L. Conformational equilibrium in the alanine dipeptide in the gas phase and aqueous solution: A comparison of theoretical results. J. Phys. Chem. 1992, 96 (9), 3864-3870. 56. Smith, P. E. The alanine dipeptide free energy surface in solution. J. Chem. Phys. 1999, 111 (12), 5568-5579. 57. Wang, Z. X.; Duan, Y. Solvation effects on alanine dipeptide: A MP2/cc-pVTZ//MP2/631G** study of (Φ, Ψ) energy maps and conformers in the gas phase, ether, and water. J. Comput. Chem. 2004, 25 (14), 1699-1716. 58. M. Bhandarkar, A. B., E. Bohm, R. Brunner, F. Buelens, C. Chipot, A. Dalke, S. Dixit, G. Fiorin, P. Freddolino, P. Grayson, J. Gullingsrud, A. Gursoy, D. Hardy, C. Harrison, J. Hénin, W. Humphrey, D. Hurwitz, A. Hynninen, N. Jain, N. Krawetz, S. Kumar, D. Kunzman, J. Lai, C. Lee, R. McGreevy, C. Mei, M. Nelson, J. Phillips, B. Radak, O. Sarood, A. Shinozaki, D. Tanner, D. Wells, G. Zheng, F. Zhu NAMD User's Guide (Version 2.11) http://www.ks.uiuc.edu/Research/namd/2.11/ug/ (accessed April 22, 2016). 59. Chipot, C.; Hénin, J. Exploring the free-energy landscape of a short peptide using an average force. J. Chem. Phys. 2005, 123 (24), 244906. 60. Bolhuis, P. G.; Dellago, C.; Chandler, D. Reaction coordinates of biomolecular isomerization. Proc. Natl. Acad. Sci. U.S.A. 2000, 97 (11), 5877-5882. 61. Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition path sampling: Throwing ropes over rough mountain passes, in the dark. Annu. Rev. Phys. Chem. 2002, 53 (1), 291-318.



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 32 of 33

62. Hazel, A.; Chipot, C.; Gumbart, J. C. Thermodynamics of deca-alanine folding in water. J. Chem. Theory Comput. 2014, 10 (7), 2836-2844. 63. Lu, C.; Li, X.; Wu, D.; Zheng, L.; Yang, W. Predictive sampling of rare conformational events in aqueous solution: designing a generalized orthogonal space tempering method. J. Chem. Theory Comput. 2016, 12 (1), 41-52. 64. Cao, L.; Lv, C.; Yang, W. Hidden conformation events in DNA base extrusions: a generalized-ensemble path optimization and equilibrium simulation study. J. Chem. Theory Comput. 2013, 9 (8), 3756-3768. 65. Wojtas-Niziurski, W.; Meng, Y.; Roux, B.; Bernèche, S. Self-learning adaptive umbrella sampling method for the determination of free energy landscapes in multiple dimensions. J. Chem. Theory Comput. 2013, 9 (4), 1885-1895. 66. Meng, Y.; Roux, B. Efficient determination of free energy landscapes in multiple dimensions from biased umbrella sampling simulations using linear regression. J. Chem. Theory Comput. 2015, 11 (8), 3523-3529. 67. Chen, M.; Yang, W. On-the-path random walk sampling for efficient optimization of minimum free-energy path. J. Comput. Chem. 2009, 30 (11), 1649-1653. 68. Chipot, C.; Lelièvre, T. Enhanced sampling of multidimensional free-energy landscapes using adaptive biasing forces. SIAM J. Appl. Math. 2011, 71 (5), 1673-1695.


Page 33 of 33

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


TOC graphic:


Extended Adaptive Biasing Force Algorithm. An On-the-Fly

Recommend Documents