although the statistical - mechanical 
 foundations of free - energy 
 calculations were laid a long time ago , their practical applications 
 became possible only with the advent of modern computers . from the 
 inception of computer - based free - energy calculations it has been clear to theorists that direct boltzmann sampling of 
 rugged energy landscapes is inefficient . the subsequent development 
 of the field is a history of efforts to remedy this problem . in free - energy calculations , the quantity of interest is almost 
 always the free - energy difference between physical 
 states of the system rather than the absolute free energy of a given 
 state . from this standpoint , calculations can be categorized on the 
 basis of variables used to transform the system between states of 
 interest . then , two main classes can be distinguished , namely alchemical 
 and geometrical transformations . they 
 rely , respectively , on changes of a parameter in the hamiltonian or 
 a function of atomic coordinates . the first class encompasses structural 
 modifications of chemical species that rest upon the remarkable malleability 
 of the potential energy function in molecular - mechanics - based simulations , reminiscent of the fabled ability of alchemists to transmute base 
 metals into noble ones . alchemical transformations are often associated 
 with the free - energy perturbation method on account 
 of the progressive and perturbative nature of the change incurred 
 by the system of interest , although , strictly speaking , alchemical 
 free - energy calculations can be carried out by way of alternate approaches , 
 such as thermodynamic integration . the 
 very first application of alchemical transformations to a nontrivial 
 chemical problem was published nearly 30 years ago by william jorgensen , 
 to whom the present contribution is dedicated . in noteworthy agreement with experiment , this pioneering simulation 
 reproduced the relative hydration free energy of methanol with respect 
 to ethane . the second class of transformations embraces virtually 
 any geometric 
 modification in a molecule or a collection of molecules by means of 
 selected collective variables tailored to address the problem at hand , 
 which could vary from changes in the internal degrees of freedom in 
 a molecule to intricate recognition and association phenomena . such collective variables form the transition 
 coordinate , a low - dimensional representation of a multidimensional 
 mathematical object . the distinction between the two types of 
 transformations is theoretically 
 important . for geometric transformations , the transition coordinate 
 is , in effect , a generalized coordinate , the evolution of which is 
 usually described by hamilton s equations of motion . in contrast , 
 for a parameter in the hamiltonian , no equations of motion naturally 
 exist , although it is possible to extend the formalism of dynamics 
 to include such a parameter . as a consequence , 
 a number of methods for calculating free energy by way of geometric 
 transformations can not be applied to alchemical transformations without 
 such extension . a considerable number of 
 ingenious techniques have been developed 
 to improve the efficiency of mapping free - energy landscapes associated 
 with geometrical transformations along a transition coordinate . a common feature of these methods is their reliance on importance - sampling 
 techniques . the central idea of these techniques is to depart from 
 sampling from the boltzmann distribution defined by the original hamiltonian 
 and , instead , sample from another distribution that favors regions 
 of phase space that would be visited only infrequently but are important 
 to achieving reliable free - energy estimates . because this procedure 
 is clearly biased , it is essential to know how to correct , or unbias , 
 it to recover the true underlying distribution . importance sampling 
 is commonly used not only in statistical mechanics of condensed phases 
 but also in other fields of science , usually as a variance reduction 
 technique most frequently combined with the monte carlo method . probably the most popular , and also the oldest importance - sampling 
 technique used in free - energy calculations is umbrella sampling . it relies on introducing a bias in simulations 
 that favors states corresponding to large values of the free energy 
 along the transition coordinate . local elevation , conformational flooding , metadynamics , and the wang landau algorithm are examples 
 of more recent importance - sampling algorithms , united by the common 
 denominator that a memory - dependent potential disfavors regions of 
 conformational space that have already been frequently visited . in a sense , the adaptive biasing force method rewove the fabric of free - energy calculations of geometrical 
 transformations , as it is characterized by both conceptual and practical 
 simplicity and requires , at least in principle , little intervention 
 from the end user . in spite of apparent similarities with the local - elevation 
 and conformational - flooding strategies in its aim to sample efficiently 
 all values of the transition coordinate , its theoretical underpinnings 
 are quite different , as we will argue further in this paper . furthermore , 
 in contrast with seemingly similar strategies , the adaptive biasing 
 force algorithm requires no prior knowledge of the free - energy landscape 
 at hand . at its core , the adaptive biasing force method is an 
 adaptive importance - sampling 
 strategy in which the quantity being adjusted is the average force 
 acting along the transition coordinate . it helps the system under 
 study escape from kinetic traps in which it would otherwise have remained 
 for a very long time . the method constitutes a highly efficient route 
 to estimating free energies , which , since its inception , has been 
 used to tackle a number of challenging problems of chemical and biological 
 interest , such as mechanical proteins , transport phenomena , or protein  ligand and protein  protein 
 recognition and association . more generally , 
 it is a versatile , adaptive , importance - sampling strategy that can 
 be utilized in many fields , whenever sampling of a probability measure 
 is thwarted by metastability of the sampling dynamics . in this 
 self - contained contribution , the multiple facets of the 
 adaptive biasing force algorithm are discussed in an exhaustive manner , 
 tackling a number of issues that have not been addressed so far , or 
 only rarely so . in the following section , we present the theoretical 
 foundations of the method , discussed in the context of other free - energy 
 approaches . next , we briefly address a number of practical issues 
 related to a proper choice of the transition coordinate . then , an 
 analysis of the convergence properties of the method and approaches 
 to calculating and controlling statistical errors associated with 
 the calculated free - energy values are presented . the discussion of 
 convergence and errors continues with the focus on nonergodicity scenarios , 
 and ways to identify and circumvent them . subsequently , we examine 
 for the first time how the adaptive biasing force algorithm is used 
 in conjunction with geometrical restraints , which have to be enforced 
 in many problems of interest . finally , we discuss some new strategies 
 for combining thermodynamics and kinetics in importance - sampling simulations , 
 before closing with recommendations for  good practices  
 in applying the adaptive biasing force method and an outlook toward 
 further promising statistical mechanical and algorithmic developments . in this section , the essential 
 idea behind the adaptive biasing 
 force algorithm is first explained in terms appealing to physical 
 intuition , followed by the theoretical underpinnings of the method 
 presented in a more formal language . then , the reader is guided through 
 the common expressions for the mean force , and the adaptive algorithm , 
 both of which are at the core of the method . the adaptive 
 biasing force method is aimed at improving the efficiency of molecular 
 dynamics simulations in which the potential energy surface is sampled 
 ineffectively due to free - energy barriers . in practice , these barriers 
 appear as bottlenecks in the dynamics of certain privileged coordinates 
 that describe the transitions between physically important states 
 ( transition coordinates ) . they also cause the system to become trapped 
 in some states for durations exceeding the time scale of the simulation , 
 resulting in incomplete sampling . the free energy along a transition 
 coordinate can be seen as a potential resulting from the average force 
 acting along the coordinate ( i.e. , the negative of the gradient of 
 this potential ) , hence the name potential of mean force . in the formalism of thermodynamic integration , 
 on which the adaptive biasing force is based , the instantaneous force acting 
 along the coordinate may be decomposed into the sum of the average 
 force ( which depends only on the value of the transition coordinate ) 
 and a random force with zero average , reflecting fluctuations of all 
 other degrees of freedom . hence , in a low - dimensional view of the 
 process , the transition coordinate evolves dynamically in its time - independent 
 potential of mean force , and this evolution is driven by the random 
 force . in many instances , the random force can be satisfactorily approximated 
 as diffusive , leading to a simple physical picture in which the system 
 diffuses along the transition coordinate in the potential of mean 
 force . the idea behind the adaptive biasing force algorithm 
 is to preserve 
 most characteristics of this dynamics , including the random fluctuating 
 force , while flattening the potential of mean force to remove free - energy 
 barriers , and thus accelerate transitions between states . this is 
 done adaptively , without any prior information about the potential 
 of mean force . to accomplish this , the instantaneous force acting 
 along the coordinate is calculated , and its running time average is 
 recorded , thus providing an on - the - fly estimate of the derivative 
 of the free energy at each point along the pathway . at the same time , 
 an external biasing force is applied , exactly canceling the current 
 estimate of the average force . over time , as the estimate converges 
 to the average force at equilibrium , the total , biased average force 
 stabilizes at values very close to zero . then , the system experiences 
 a nearly flat potential of mean force and displays accelerated dynamics 
 along the transition coordinate . the fact that the biasing force is 
 exactly equal to the mean force is actually not crucial . what is important 
 is that the biasing force yields sufficiently uniform sampling of 
 the transition coordinate that the remaining barriers can be easily 
 traversed in response to thermal fluctuations . the adaptive biasing force 
 algorithm is not inherently tied to any specific type of dynamics 
 but does rely on sampling of the canonical ensemble . in explicit - solvent 
 simulations , langevin dynamics with sufficiently soft damping and 
 small stochastic forces becomes a mere perturbation of hamiltonian 
 dynamics and may be used as one simple way to achieve canonical sampling . 
 for convenience , but without loss of generality langevin dynamics can be written 
 as1where 
 ( xt ,  pt ) denotes the positions and 
 momenta of the particles at time t , m is the mass tensor , v :   is the potential energy function ,  
 is the friction coefficient , wt is the wiener process that underlies the random force ( white 
 noise ) , and  = kbt , where kb is the boltzmann 
 constant and t is temperature . the dynamics of eq 1 is ergodic ( under mild conditions on v ) with respect to the canonical measure zp1  exp[pmp/2 ]  dp  (dx ) , where (dx ) = zx1  exp[v(x ) ]  dx . ergodicity means 
 that long - time averages converge to canonical averages:2and that the law at time t of the stochastic process converges to the canonical measure 
 in 
 the long - time limit:3where   is 
 the expected value . the first limit 
 in eq 2 is of particular interest for practical 
 applications because it allows for computing canonical averages from 
 trajectory averages . for a typical potential v , the associated boltzmann measure  is multimodal : high - probability 
 regions are separated by low - probability regions . the former correspond , 
 for instance , to the most likely conformations of a biological object , 
 which are typically separated by transition regions of very low probability . 
 for these reasons , estimating averages with respect to the probability 
 measure  is , in general , a difficult task . in particular , the 
 ergodic properties of the dynamics of eq 1 are 
 not sufficient to devise reliable numerical methods , because under 
 these premises , the stochastic process remains trapped in large - probability 
 regions , and , as a consequence , the long - time asymptotic regime t   is very difficult to reach in the ergodic 
 limits in eqs 2 and 3 . the fact that the system remains for a very long time in some region 
 of phase space before hopping to another region is called metastability , 
 and the corresponding states of the system are called metastable . 
 the inability to reach the ergodic limit is often called quasi nonergodicity ; 
 the system appears nonergodic on the time scales of the simulations . a typical example of a metastable state is a local free - energy minimum 
 in the conformational space of a protein . the adaptive biasing 
 force method relies on modifying the potential v in 
 such a way that the energy landscape is flattened along 
 a given transition coordinate  : . here , we restrict ourselves to a one - dimensional 
 transition coordinate , leaving generalization to high - dimensional 
 transition coordinates for the section expressions 
 for the mean force . at[(x ) ] and at is 
 updated in such a way that it converges to the free energy a , defined ( up to an additive constant ) by4where the measure (x)z  ( dx ) is supported by the subset 
 { x , (x ) = z } and 
 is such that (x)z  ( dx ) dz = dx . in practice , the bias is only applied in a window [ zmin , zmax ] as explained 
 in 
 the section justification of a stratification strategy . notice that , by the definition of a , the canonical 
 measure associated with the biased potential v(x ) = c , where c is a constant 
 independent of z. therefore , if the biased potential v(x )  a[(x)]1(x)[zmin , zmax ] is 
 used , the marginal along  is a uniform 
 law over [ zmin , zmax ] . here , 1(x)[zmin , zmax ] is 
 an indicator function , which is one when (x )  
 [ zmin , zmax ] and zero otherwise . let us recall that the marginal law is defined 
 as follows : if x is distributed according to a probability 
 distribution  , then the law of (x ) is called 
 the marginal of  along . if we knew the free energy a , it would be a good idea to use a as a biasing potential , because sampling along  
 would be easier . this is illustrated in simple two - dimensional examples 
 in figures 1 and 2 . notice in particular that the free energy seems to be a good biasing 
 function for efficient sampling of both energetic barriers ( figure 1 ) and entropic barriers ( figure 2 ) . upper row : ( a ) original two - dimensional , double - well potential 
 displayed as level sets ; ( c ) time trajectory of the first coordinate x in the stochastic process , showing oscillations between 
 the two metastable wells . lower row : ( b ) level sets of the same potential 
 biased by the free energy associated with the transition coordinate 
 (x , y ) = x ; ( d ) time trajectory of the transition coordinate x in the adaptively biased dynamics , showing no metastability . upper row : ( a ) the original two - dimensional 
 potential is zero inside 
 the hourglass shape , and +  outside ; ( c ) time trajectory of 
 the first coordinate x in the stochastic process , 
 showing oscillations between the two metastable wells . lower row : 
 ( b ) free energy along the transition coordinate (x , y ) = x , featuring a purely entropic 
 barrier ; ( d ) time trajectory of the transition coordinate x in the free - energy - biased dynamics , showing no metastability . the main ingredient that we now 
 need is an update rule for at such that limtat = a. this is based on the following formula , which defines 
 the mean force , namely , the negative of the derivative 
 of the free energy . more general formulas can be derived , which give 
 rise to many variants of the adaptive biasing force method ( see section expressions for the mean force),5where the instantaneous force is defined by6 equation 5 is a consequence of the definition 
 of the free energy in eq 4 and of the co - area 
 formula ( which is a generalization of the fubini theorem ) ; see for 
 instance ref ( 8) . from eq 5 it follows that a 
 is the conditional average of the instantaneous force , f , with respect to the canonical measure conditioned 
 to a fixed value of the transition coordinate : a(z ) = [ f(x)(x ) = z ] . an important observation is that 
 eq 5 remains 
 true if the potential v is changed to v at : 
 for any function at7 in other words , a biasing potential at depending solely on  leaves averages 
 conditioned 
 by  unchanged . this is the intuition behind the adaptive biasing 
 force dynamics , which can be written as8 indeed , if ( xt ,  pt ) were at equilibrium with respect 
 to the biased canonical measure zp1  exp(pmp/2 ) dpzt1  exp[(v at)(x ) ] dx , then at would be equal to a. in practice , 
 of course , the sampling is not sufficiently fast for the process to 
 be instantaneously at equilibrium with respect to the time - varying 
 biased potential v at , and this is why this heuristic 
 is not sufficient to fully understand the convergence of the adaptive 
 biasing force dynamics ( see the section convergence 
 and error analysis ) . however , this simple reasoning is sufficient 
 to check that if at converges 
 to some limit a , then necessarily , a = a. from a practical 
 viewpoint , the conditional expectation in eq 8 can be computed using two procedures : either time 
 averages over a single long trajectory or averages over many replicas 
 run in parallel . these procedures will be discussed in ample detail 
 in the section multiple - walker strategies below . intuitively , the adaptive biasing force dynamics thus 
 consists 
 of adding a force at[(x)](x ) that exactly compensates the average of the original force , 
 v(x ) , along a given 
 transition coordinate . if  is well - chosen , the hope is to observe 
 a fast convergence ( at least compared to the original dynamics embodied 
 in eq 1 at equilibrium ) of at to the mean force a. given a transition coordinate 
 (x ) , the mean force is a well - defined quantity , 
 yet its expression as an ensemble average of an instantaneous force f is not unique , as we will see below . 
 in adaptive biasing force simulations , the choice of a convenient 
 expression for f is driven by practical 
 considerations , notably ease of implementation and numerical behavior , 
 such as variance . the classic expression for the instantaneous 
 force involves an explicit coordinate transformation  from 
 cartesian to generalized coordinates , which include the transition 
 coordinate . that is ,  : , with 1 =  and 
 components i for i > 1 are generalized coordinates of no particular physical significance , 
 but necessary to the mathematical framework . a valid expression for 
 the instantaneous force is then9which in the physics literature is 
 more commonly 
 written as10where |j| is the determinant 
 of the jacobian matrix ( ij1)(i , j ) . from the arbitrary choice of 
 i , i>1 is derived 
 a ( somewhat arbitrary ) function f , whose -restricted ensemble average nevertheless yields the 
 uniquely defined mean force ( eq 5 ) . equation 10 provides an intuitive view that different choices 
 of  and f correspond to 
 different ways of projecting the cartesian forces , v , onto the transition coordinate . the direction along which 
 the forces are projected in this expression is the vector 1 , which we call  inverse 
 gradient  . as the gradient of can be seen as the changes in  corresponding to infinitesimal 
 changes in x , the inverse gradient is the vector along 
 which a change in  is propagated in cartesian coordinates , 
 other generalized coordinates ( i , i>1 ) being constant , hence the dependence of the 
 inverse 
 gradient on the explicit coordinate transformation . the alternate 
 notation for the inverse gradient , x/ , 
 has the drawback of hiding this dependence on the choice of . numerically , eq 10 is impractical for two 
 reasons . one is that defining i , i>1 explicitly can be exceedingly difficult , especially 
 if  is a collective coordinate ( e.g. , the radius of gyration 
 of a set of particles ) . supposing that this step is done , the second 
 difficulty comes with the numerical computation of the jacobian derivative , 
 as it involves second derivatives of  whose 
 analytical derivation and numerical implementation may be cumbersome , 
 again , depending on the nature of the transition coordinate . to circumvent this issue , the original adaptive biasing force method 
 was introduced with an instantaneous force estimator based on a constraint 
 force that is calculated iteratively , but never applied . in the initial implementation of the adaptive 
 biasing force algorithm in the popular 
 molecular dynamics program namd , eq 10 was used because the small set of implemented coordinates 
 made it practical . as this set was greatly extended in the framework 
 of the collective variables module , more 
 versatile expressions of the instantaneous force were required . den otter put forward the idea that the 
 inverse gradient can be 
 replaced with an arbitrary vector field ( satisfying certain requirements ) . in other words , changes in  may be propagated 
 along an arbitrary direction in cartesian coordinates , without explicitly 
 defining a complete set of generalized coordinates . consider a vector transition coordinate ( i ) , in the presence of a set of constraints 
 of the form k(x ) = 
 0 . for each coordinate i , let vi be a vector field   satisfying , for all j and k:1112 the ith partial derivative of the free energy 
 can then be calculated as the conditional average of the following 
 instantaneous force:13of which 
 eq 6 is a special 
 case , limited to a single coordinate and choosing v = 
 /|| , which satisfies the 
 condition ( eq 11 ) above . note that this estimator 
 still requires the calculation of second derivatives , in the form 
 of the divergence of the vector field v , although the 
 relative freedom in choosing v can be taken advantage 
 of to make the divergence calculation practical . the choice v = /|| is always 
 valid , and as such , convenient for theoretical purposes , but certainly 
 not always optimal when implementing specific generalized coordinates . expressions that were chosen in practice for those coordinates implemented 
 in the collective variables module are listed in ref ( 39 ) . described 
 an estimator that does not require second 
 derivatives , but rather first derivatives with respect to time and 
 space , and is valid for multidimensional adaptive biasing force calculations . this estimator resembles a statistical form 
 of newton s equation of motion : instead of relating acceleration 
 to the potential energy gradient , it relates the mean acceleration 
 to the gradient of the free energy . in a considerable simplification , only the first derivative of the force with respect to  needs 
 to be derived . the time derivative is calculated numerically in the 
 same fashion and at the same level of accuracy as time derivatives 
 of other quantities in molecular dynamics . other ways of simplifying 
 calculations of instantaneous forces will be discussed in the context 
 of extended adaptive biasing force simulations , or eabf , in the section the extended adaptive biasing force method . the final , essential ingredient 
 of the theory underlying the adaptive force method is an algorithm 
 for deriving the current estimate of the average force as simulations 
 progress . in its generic , one - dimensional implementation , the transition 
 coordinate ,  , connecting two end points , is divided into m equally sized bins of width  in which forces are accrued in the course of the simulation . in 
 a nave approach , the approximation to the average force , f(nstep , k ) in bin k after nstep molecular dynamics steps is just the simple , unweighted 
 average of all force samples in this bin14provided that the bin has already been visited 
 at least once . nstepk is the number of samples 
 accrued in bin k after nstep steps and fk abbreviates the th 
 force sample in this bin . however , when only a few samples are available 
 in a given bin k , the running average might be a 
 poor estimate of the actual average force in this bin . large fluctuations in the running estimate of the average force 
 are undesirable , as they may drive the system away from equilibrium , 
 thus slowing the convergence of the algorithm and reducing the efficiency 
 of the method . to control these effects , a procedure is needed to 
 reduce variations in early estimates of f(nstep , k ) . a number of schemes can be applied for this purpose . in current 
 implementations , the biasing force in bin k at time t is applied in full only if the number of samples ntk is above a threshold , nfull . it is ramped up smoothly as ntk varies from 0 to nfull . in one implementation , the 
 ramp is linear and the force is proportional to ntk / nfull ; in another , the biasing force is zero for ntk < nfull/2 and is ramped linearly 
 above that value , proportionally to 2ntk / nfull  1 . both 
 implementations have proven to be efficient in a number of applications , 
 but other , more advanced schemes are possible . so far , there have 
 been no systematic studies on the efficiency of different adaptive 
 algorithms , but it is anticipated that it may be strongly system - dependent . for a sufficiently large nstep , f(nstep , k ) approaches the correct average force in each 
 bin . then , the free - energy difference , a , between the end point states can be estimated simply 
 by way of summing the force estimates in individual bins.15if the average force is 
 a rapidly changing 
 function of  , a more sophisticated integration algorithm may 
 be required . it might be also possible to develop binless integration 
 algorithms , similar to those proposed for some other free - energy calculation 
 methods . a common trait of importance - sampling 
 algorithms is the discretization 
 of the transition coordinate ,  , in bins of width  in which statistical information is accrued in the course of the 
 simulation . in the umbrella - sampling scheme , for instance , a histogram is constructed , corresponding 
 to the biased probability of occurrence of the molecular assembly 
 of interest at the different values of the transition coordinate . 
 in the adaptive biasing force algorithm , bins are utilized to store instantaneous values of the thermodynamic 
 force that acts along the transition coordinate . as has been observed 
 in practice previously for diffusive dynamics , the instantaneous force 
 in any given bin obeys a normal distribution . at thermodynamic equilibrium , by definition , its average is exactly 
 equal to a(z ) , that 
 is , the gradient of the free energy along the transition coordinate . 
 insisting upon being at thermodynamic equilibrium is pivotal here , as application of a poorly estimated time - dependent 
 bias , i.e. , from a distribution out of equilibrium , will not yield 
 a hamiltonian bereft of a mean force acting along the transition coordinate . the adaptive algorithms described above contain two adjustable 
 parameters : bin width ,  , and the number 
 of samples , ninit , below which r(nstepk ) is not equal to 1 or , equivalently , 
 averaging does not follow eq 14 . the choice 
 of these parameters should be done with some care . at constant  , large values of ninit yield better estimates of the average force 
 once the number of samples collected in a given bin reaches this threshold 
 value and conventional averaging begins , at the price of delayed application 
 of the full averaging and slow , initial progress along the transition 
 coordinate . small values of ninit , in 
 turn , tend to drive the system out of equilibrium . typically , setting ninit in the range between 200 and 500 appears 
 to be a good compromise that allows for avoiding both types of problems . at constant ninit , large values of this may have adverse effects on the 
 accuracy of integration in eq 15 . on the other 
 hand , small values of  require longer 
 simulation times to collect sufficient force statistics in every bin . 
 if the transition coordinate is a distance in an atomistic system , 
 the choice of  as 0.1 or 0.2  
 usually represents a satisfactory trade - off . in the last four decades a 
 number of strategies have been developed 
 for computing free - energy changes as a function of geometrical transformations , 
 each endowed with advantages , as well as limitations . to a large extent , umbrella sampling , 
 whether in its original formulation or variants thereof , remains one of the most widely 
 utilized routes to address rare events in molecular simulations . in 
 its original form , umbrella sampling referred to incorporating an 
 external biasing potential in the simulation , i.e. , an umbrella , ideally the negative of the free energy , which 
 would yield a broad , if not uniform exploration of the transition 
 coordinate . in other words , in ideal circumstances , the system would 
 evolve in the collective - variable space on a flat free - energy hyperplane , 
 as is also the case for the adaptive biasing force method . under most 
 circumstances , however , the form of the optimal umbrella is unknown . thus , for any qualitatively new problem , the end - user must resort 
 to an educated guess regarding the shape of the biasing umbrella potential , 
 usually on the basis of prior knowledge of this and related problems . this may constitute a daunting task . poorly 
 predicted biasing potentials yield nonuniform probability distributions 
 across the transition coordinate . this decreases the efficiency of 
 umbrella sampling , which , in extreme cases , may reduce rather than 
 improve the efficiency compared to results with unbiased calculations . this common shortcoming led to the development of an adaptive variant 
 of the umbrella - sampling algorithm , wherein 
 the initial guess of the biasing potential is progressively refined 
 in light of a series of short simulations . adaptive umbrella 
 sampling is one member of a broader family of techniques called adaptive 
 biasing potential methods . local elevation , conformational flooding , or its more recent avatar , metadynamics , adaptive biasing molecular dynamics , and the wang the idea is to penalize the 
 already visited states by changing the potential v to v at , at(z ) being related to the occupation time of the 
 value z of the transition coordinate up to time t. the longer the time spent in a bin { x , (x )  [ z0 , z1 ] } , the larger the biasing potential at(z ) , z  [ z0 , z1 ] . the biasing potential at is built as a sum of gaussian kernels that are periodically 
 added to the hamiltonian along the  variable . this pushes the 
 system away from states that have already been visited and , by doing 
 so , improves sampling . it ought to be noted that the potential at , rather than its derivative , 
 is computed in these two cases , hence , the name adaptive biasing potential 
 methods . one important downside of the adaptive biasing force algorithm , 
 compared to the class of adaptive biasing potential methods , lies 
 in its inability to handle discrete transition coordinates , for instance , 
 coordination numbers . this drawback can be understood by considering 
 that the free energy is now a map from integers to reals , and , thus , 
 has no derivative and , hence , no mean force . from a mathematical 
 viewpoint , the adaptive biasing force method , 
 just like adaptive biasing potential methods , is an adaptive importance - sampling 
 procedure . there is , however , a salient difference between these two 
 techniques . in the latter , the potential of mean force or , equivalently , 
 the corresponding probability distribution along the transition coordinate 
 is being adapted . in contrast , the former relies on biasing the force , 
 i.e. , the gradient of the potential . this difference is more important 
 than it might appear at first sight , as potentials and probability 
 distributions are global properties whereas gradients are defined 
 locally . in terms of probability distributions , it means that the 
 count of samples in the neighborhood of a given value of the transition 
 coordinate is insufficient to estimate probability . knowledge of the 
 underlying probability distribution over a much broader range of  
 is required . this may considerably impede efficient adaptation . in 
 contrast , all that is needed to estimate the gradient is the knowledge 
 of local behavior of the potential of mean force . thus , in many 
 instances , adaptation proceeds markedly faster . using a common metaphor , the difference between the adaptive biasing potential and adaptive 
 biasing force methods can be compared to inundating the valleys of 
 the free - energy landscape as opposed to plowing over its barriers 
 to yield an approximately flat terrain , conducive to unhampered diffusion . for example , in metadynamics and its ancestors , 
 the widths and weights of gaussian functions and the frequency with 
 which the biasing potential is updated have to be carefully chosen , 
 which often requires considerable experience . in the adaptive biasing force method , estimating the biasing gradient happens automatically 
 by way of a simple algorithm , described in the previous subsection . another technical concern about adaptive methods is to ensure that 
 adaptation vanishes once at approaches the converged free energy . there are a number of ways to fulfill this condition more - or - less 
 automatically , but adaptive biasing force and 
 adaptive biasing potential techniques remain intrinsically different 
 from this point of view . in the adaptive biasing force algorithm , 
 if the correct free energy is given as an initial guess ( namely if v is replaced by v a in eq 8) , then the biasing 
 force will not be updated ( at in eq 8 will be constant over time ) . moreover , if the derivative 
 of the biasing potential is needed ( for example to bias the langevin 
 dynamics as in eq 8) , the advantage of the adaptive 
 biasing force algorithm is that at is directly 
 computed , whereas in adaptive biasing potential algorithms , one needs 
 to differentiate the evaluated biasing potential at , which may lead to very noisy results 
 because at is estimated 
 along a stochastic trajectory . because the basic quantity calculated , 
 and subsequently integrated , 
 in the adaptive biasing force method is the force , this approach belongs 
 to the thermodynamic integration class of methods . however , in contrast 
 to conventional implementations of thermodynamic integration and its 
 generalizations , such as the blue - moon ensemble approach , the adaptive biasing force algorithm does not 
 rely on constrained molecular dynamics , but instead is based on unconstrained 
 simulations ; i.e. , the free - energy difference is not determined at 
 discrete values of the transition coordinate through solving constrained 
 equations of motion . sampling of the transition pathway proceeds in 
 a continuous , unhampered fashion , guided by the diffusion properties 
 of the system of interest , obviating the need for re - equilibration 
 at fixed , predefined values of the transition coordinate , even in 
 stratified simulations . as will be discussed further in this paper it is also 
 of interest to compare the adaptive biasing force method 
 with reconstructions of free - energy landscapes from nonequilibrium 
 trajectories that represent repeated pulling experiments . the latter are based on the groundbreaking jarzynski identity , or its extension to bidirectional transformations , combined with steered molecular dynamics . even 
 though both approaches involve molecular dynamics trajectories that 
 are initially away from equilibrium , there is a fundamental distinction 
 between them . in the adaptive force method , only the average , or systematic 
 force is removed to erase the ruggedness of the free - energy landscape , 
 preserving the random force responsible for diffusion . once a good 
 estimate of the average force becomes available , the equilibrium behavior 
 of the system is restored . in pulling experiments , the transformation 
 always proceeds away from equilibrium at constant velocity and , therefore , 
 the instantaneous force acting along the transition coordinate is 
 nil . the random force actually appears in the formalism after averaging 
 over the ensemble of pulling experiments . moreover , achieving convergence 
 in calculations based on jarzynski s identity usually requires 
 large numbers of independent realizations , which comes at a significant computational cost . taken together , 
 when a geometrical transformation can be undertaken at equilibrium , 
 it is not clear whether there is any practical advantage of handling 
 the problem at hand by means of nonequilibrium work experiments rather 
 than the adaptive biasing force method . an important advantage 
 of gradient - based methods is 
 the possibility of formally decomposing the free - energy change into 
 physically meaningful contributions , thereby helping 
 to dissect qualitatively the nature of the intermolecular interactions 
 at play . it is worth noting that different energy terms contribute 
 to the mean force both explicitly through force terms and implicitly 
 through the boltzmann weights in the canonical average ; contributions 
 can only be separated numerically at the former level , not at the 
 latter . decomposition of the free energy is generally handled a posteriori 
 through computing the thermodynamic force between , for instance , groups 
 of atoms of interest . this force is then projected onto the transition 
 coordinate determined for each stored configuration , prior to the 
 construction of a histogram from which the average force is inferred . integration of the latter yields the desired contribution to the total 
 free - energy change across the entire transition pathway . among 
 many options for reconstructing free - energy landscapes along a transition coordinate , which one is the best ? considering that the 
 efficiency of different methods strongly depends on their implementation 
 in software packages and , very likely , on a system of interest , attempts 
 to answer this question appear somewhat misguided and unproductive . that said , one aspect of the adaptive biasing force algorithm pleading 
 in its favor is its simplicity . how the algorithm operates 
 is physically intuitive , requiring , in principle , 
 very little prior knowledge of the free - energy landscape , or input 
 from the end - user , even for qualitatively new problems . central to geometric transformations is the concept 
 of a transition 
 coordinate . in this section , this concept is illuminated in the context 
 of free - energy calculations aimed at tackling rare events . specifically , 
 we will discuss how the transition coordinate is explored with the 
 adaptive biasing force algorithm and delve into the practical aspects 
 of defining this coordinate . stratification , a common technique for improving the efficiency of free - energy calculations 
 by partitioning the reaction pathway into ranges of the transition 
 coordinate , will be discussed with the focus on the justification 
 for and limitations of this strategy . for any transition 
 from a macrostate , a , to another macrostate , b , of the same system 
 there exists an exact one - dimensional transition 
 coordinate : the committor probability . in most cases , this coordinate 
 is difficult to calculate and usually offers very limited insight 
 into the nature of the process of interest . for these reasons , it 
 is often more useful to employ a transition coordinate that is only 
 an approximation to the committor probability but is physically more 
 meaningful and easier to handle . sometimes it might be helpful to 
 extend the reduced representation of the transition to a transition 
 coordinate that extends beyond one dimension . not only physical significance 
 but also efficiency of sampling , given a transition coordinate , is 
 of concern . as 
 in any enhanced sampling method based on a reduced representation , 
 adaptive biasing force sampling relies on stochastic exploration along 
 the transition coordinate ,  , enhanced by the adaptive bias , 
 combined with equilibration of other , orthogonal degrees of freedom . equilibration in the orthogonal space is critical in two respects : 
 it affects the mobility along  ( see figure 1b for a diffusive example ) , and it determines the rate of 
 convergence of a(z ) = f(x) , which is an average over the orthogonal degrees of freedom . in 
 the ideal situation of time scale separation , all slow degrees of 
 freedom are captured by the transition coordinates , and relaxation 
 in the orthogonal space is comparatively fast . in other words , the 
 adaptive biasing force algorithm removes metastability along the transition 
 coordinates , provided that other degrees of freedom are not metastable . complex systems such as biological macromolecules , however , possess 
 many slow , coupled degrees of freedom , making time scale separation 
 difficult or impossible to achieve . fortunately , empirical results 
 suggest that time scale separation 
 is not an absolute requirement of adaptive biasing force sampling . 
 one reason behind it might be that enhanced diffusion in the transition 
 coordinate space reduces metastability in the orthogonal space , by 
 letting the dynamics sidestep orthogonal barriers 
 rather than cross them , making some  multichannel cases 
 ( see the section hidden barriers and other challenges 
 to obtaining accurate results ) tractable with standard adaptive 
 biasing force simulations . more encouraging still , convergence in 
 such multichannel cases can be markedly accelerated by multiple - walker 
 formulations of the adaptive biasing force algorithm ( see the section multiple - walker strategies ) . yet as will be extensively discussed further in this paper , the end - to - end 
 distance of the -helical deca - alanine peptide provides a good 
 example of an inadequate transition coordinate . depending on the range 
 of values , the coordinate exhibits completely different behavior . 
 for values corresponding to the stable -helix ( 14  ) and 
 larger , separation of time scales is obeyed and the adaptive biasing 
 force converges well . in contrast , smaller 
 values of the end - to - end distance correspond to a rich set of metastable , 
 compact states that are not resolved by the transition coordinate . as a result , adaptive biasing force dynamics 
 becomes trapped in these states , and the free - energy estimator does 
 not converge in accessible simulation times . attempts to resolve these metastable states by two- and three - dimensional 
 coordinates give improved results , allowing the exploration of all 
 metastable basins , yet those basins are not all resolved even in three 
 dimensions . these difficulties 
 are not due to specific deficiencies of the adaptive biasing force 
 method , but rather , to shortcomings of the reduced representation , 
 which would constitute an obstacle to any sampling method . in practice , finding an effective reduced representation is still a process largely 
 guided by physical intuition about the process of interest , as well 
 as trial - and - error . more systematic and robust approaches for this 
 dimension reduction step are an area of active research . a limiting factor is often the availability 
 of usable numerical implementations of the generalized coordinates 
 of choice . the collective variables module is an attempt to overcome this limitation , by providing a rich and 
 flexible toolbox to define many types of coordinates , in particular 
 those useful for the description of biological macromolecules . in 
 this module , once an intuitive understanding of the relevant 
 generalized coordinate has been obtained , some technical choices remain 
 to be made to express this coordinate as a function  of atomic 
 cartesian coordinates . though these decisions may seem ancillary , 
 they have a strong influence on the accuracy , convergence , and computational 
 performance of the adaptive biasing force algorithm . when objects 
 of interest are composed of many atoms , there are 
 often several nearly equivalent ways to define the transition coordinate . 
 in sufficiently long simulations , different definitions produce nearly 
 identical potentials of mean force , but the efficiency might vary 
 considerably . for example , the distance between two proteins could 
 be defined by selecting one central atom in each protein , or by selecting 
 the centers of mass of large groups of atoms . the largest contribution 
 to the instantaneous force on each atom in a molecule is due to rapidly 
 oscillating bonded interactions ; more generally , in all applications , 
 forces on the particles will contain some background noise . if many 
 atoms contribute to the projected force ( e.g. , the first right - hand - side 
 term in eq 13 ) , contributions from those noisy 
 terms average out , which lowers the variance of the instantaneous 
 force estimator . in the initial stage of an adaptive biasing force simulation , nonequilibrium 
 effects occur if the biasing force applied to some degrees of freedom 
 varies faster than coupled degrees of freedom can relax . collective variables 
 that involve many cartesian coordinates with smoothly varying contributions 
 to the gradient ( hence , to the biasing force vector ) . at equilibrium , 
 a smooth motion may be described by a nonsmooth variable and this 
 may make the convergence of the adaptive bias more difficult . in the 
 deca - alanine stretching toy example , the geometric process of interest 
 involves all atoms in the peptide , yet our classic approach uses the 
 end - to - end distance as a biasing coordinate , with the implicit assumption 
 that biasing forces exerted on the terminal atoms propagate , and that 
 the entire peptide relaxes rapidly , so that the biased trajectory 
 remains close to equilibrium . a more robust approach is to replace 
 that coordinate with the radius of gyration of the peptide . the gradient 
 of the radius of gyration has components on each atom proportional 
 to its distance from the center of the group , so that atoms close 
 to the center also experience moderate biasing forces and do not lag 
 behind the terminal atoms when the biasing force varies rapidly over 
 time . a more elaborate discussion of this problem can be found in 
 ref ( 42 ) . in 
 some cases , however , coordinates involving large collection 
 of atoms will have less resolution than more local ones . a biophysical 
 example is permeation through an interface , such as a lipid bilayer . the common choice of transition coordinate is the distance of the 
 permeant molecule to the center of the bilayer center , projected onto 
 the bilayer normal . in such a case , the physically relevant phenomenon 
 is interaction of the permeant with the membrane surface , on a local 
 scale . if the bilayer patch is large enough , it will experience fluctuations 
 away from planarity , thus the local position of the interface will 
 fluctuate with respect to the bilayer center . in turn , this will cause 
 spurious fluctuations in the transition coordinate . a comparable situation 
 arose in a study of glycerol permeation through the water channel 
 protein glpf . interaction of the permeant 
 with the protein depended on distances to neighboring pore - lining 
 residues , which fluctuated with respect to the bulk of the protein . 
 therefore , the global coordinate measured between the protein and 
 glycerol molecule had insufficient resolution on a local scale , and 
 a more local coordinate had to be defined to resolve the structure 
 of the free - energy profile in the constricted region of the selectivity 
 filter . depending on system 
 size and implementation 
 details , choices of coordinates may impact performance noticeably . a common application case is a biomolecular simulation performed with 
 the namd package , with the adaptive biasing 
 force algorithm implemented in the collective 
 variables module . namd is highly parallelized 
 and can simulate large systems on supercomputers with nearly linear 
 scaling . the current implementation of the adaptive biasing force 
 algorithm is not parallelized and runs on one node , leading to two 
 potential bottlenecks : ( 1 ) poor serial performance on the master node : 
 for the most expensive variables ( e.g. , those involving sums on atom 
 pairs ) , the bias calculation on the master node might take longer 
 than the force calculation on other nodes . ( 2 ) scaling may suffer 
 even for computationally simple coordinates , if too many atom coordinates 
 and forces have to be communicated across nodes , increasing latency . this second case may affect any highly parallel application with coordinates 
 defined on many atoms . in practice , one often has to find an acceptable 
 trade - off between variance and performance by selecting a reasonable 
 number of cartesian coordinates that are most representative of the 
 quantity of interest . to describe conformational fluctuations of a 
 protein , for example , the root - mean - square deviation of  carbon 
 coordinates is often a good compromise . to increase 
 the efficiency of exploring the transition coordinate in adaptive 
 biasing force or umbrella sampling , it is common to break down the transition path 
 into a series of sequential strata or windows . this 
 idea arises from the intuition that the 
 time to convergence grows as the square of the range of the transition 
 coordinate . convergence 
 of the free energy over the 
 entire range   is achieved 
 after t0 . let us now divide the transition 
 path into n nonoverlapping windows of lengths , for which convergence 
 is attained after t1 , ... , tn. as shown in the supporting information , t0 > iti. we illustrate this result 
 in a simple example of a tagged water 
 molecule diffusing in a bulk environment over a stretch of 20  . the transition coordinate is the projection of the distance separating 
 the centers of mass of the tagged water molecule and the simulation 
 cell along a given direction of cartesian space . translational invariance 
 due to the isotropic nature of the liquid imposes the condition that 
 the free - energy change along the transition coordinate be zero . for 
 this system , the potential of mean force was determined in a single , 
 20  long stratum , two 10  strata , four 5  strata , 
 and eight 2.5  strata . the simulations continued until the root - mean - square 
 deviation between the computed potential of mean force and the reference 
 zero free - energy profile was less than 0.1 kcal / mol . as can 
 be observed in figure 3 , t0 , the time necessary to attain convergence 
 within this preset tolerance without stratification is on the order 
 of 100 ns . when the transition coordinate is divided into two strata , 
 convergence in each 10  stratum is reached in approximately 
 40 ns , i.e. , in 80 ns over the full 20  range . the effect of 
 stratification increases further , as the reaction coordinate is decomposed 
 into four and eight windows . in these cases , convergence is achieved 
 in approximately 6 and 2 ns per window , respectively , which corresponds 
 to 24 and 16 ns for the complete 20  range . stratification strategies 
 for the translationally invariant toy 
 model of a tagged water molecule diffusing in a bulk aqueous medium . the transition path spans 20  and is handled in a single window 
 ( a ) , in two 10  windows ( b ) , in four 5  windows ( c ) , and 
 in eight 2.5  windows ( d ) . for each stratification strategy , a potential of mean force calculation is carried out until the root - mean - square 
 deviation with respect to the accurate zero free - energy profile is 
 less than 0.1 kcal / mol . both theoretical considerations and a simple example given 
 above 
 appear to suggest that extensive stratification should be always preferred . first , simulations in 
 each window require initial equilibration , which may erase benefits 
 gained from stratification into many windows . perhaps more importantly , 
 extensive stratification may impede ergodic sampling of the phase 
 space . this behavior , shared with the standard thermodynamic integration , 
 will be discussed in section addressing nonergodicity 
 scenarios . in contrast to umbrella sampling and its adaptive 
 variants , the 
 adaptive biasing force algorithm does not require that consecutive 
 windows of a dissected transition coordinate overlap . provided that 
 convergence has been achieved in each window , the gradient , a( ) , can be reconstructed merely by joining the gradients 
 from individual windows at the boundaries . this improves efficiency , 
 as the requirement for overlap between windows may frequently add 
 as much as 50% to the total simulation time . if the gradient between 
 consecutive windows is not continuous to within statistical error , 
 this is usually a sign of difficulties with ergodic sampling . again , 
 this problem will be considered in section addressing 
 nonergodicity scenarios . in this section , we 
 examine the convergence properties of the adaptive biasing force algorithm 
 and the reliability of the computed free - energy estimates . first , 
 the reader is invited to follow a demonstration that the numerical 
 scheme formally converges . then , we discuss how statistical errors 
 associated with the reconstructed free - energy landscapes can be measured 
 and managed . the aim of this section is to explain convergence properties of 
 the adaptive biasing force algorithm . for more details , we 
 consider the overdamped langevin dynamics,16where xt is defined in the n - dimensional 
 torus   ( namely , in [ 0 , 1 ] with periodic boundary 
 conditions ) and (x1, ... we direct the reader to refs ( 67 ) and ( 68 ) for extensions to more 
 general situations of the results presented below . starting 
 from eq 16 and using the above choice 
 of the transition coordinate , the adaptive biasing force dynamics 
 can be represented as17where xt1 denotes the 
 first coordinate of the vector xt , e1 is the vector with coordinates ( 1 , 0 , ... , 0 ) , 
 and 1v denotes the partial derivative 
 of v(x1, ... ,xn ) with respect to x1 . as explained in the section theoretical backdrop , it is clear at least formally that the 
 only possible stationary 
 state for at is the free 
 energy ( up to an additive constant ) . indeed , if at converges to some stationary state a , then the law of xt converges to z1  exp({v(x ) a[(x ) ] } )  dx , which implies that ( 1v(xt)xt1=x1 ) is a(x1 ) . 
 yet , this does not provide a proof of the convergence of at to a , and it does not explain 
 why the adaptive biasing force dynamics of eq 17 indeed converges faster to equilibrium than the original , unbiased 
 dynamics of eq 16 . one way to understand for the original dynamics 
 of eq 16 , the density  satisfies the 
 fokker planck equation:18for the adaptive biasing force dynamics , 
 the 
 density  satisfies19 it ought to be noted that eq 19 is a 
 nonlinear 
 partial differential equation ( pde ) , which makes the study of its 
 long - time behavior much more complicated than for the linear fokker using appropriate mathematical 
 tools ( namely , entropy techniques , 
 see the supporting information ) , one can 
 show that if the transition coordinate is well chosen , the convergence 
 to equilibrium for the adaptive biasing force dynamics of eq 19 is much faster than for the original unbiased dynamics 
 of eq 18 . roughly speaking , the assumption on 
 the transition coordinate is that the canonical measure z  exp[v(x ) ] dx is  more multimodal  
 than the conditional measures at a fixed value x1 of the transition coordinate this 
 is typically the case for the simple 
 illustrative two - dimensional potentials in figures 1 and 2 if (x1,x2 ) = x1 . this can be fully quantified , and it actually gives 
 a way to measure the quality of the transition coordinate . in 
 the analysis outlined above , it is assumed that the conditional 
 expectation appearing in the adaptive biasing force dynamics is computed 
 exactly . this analysis is therefore well adapted to discretizations 
 that involve many replicas in parallel , which indeed converge to the 
 adaptive biasing force dynamics with the exact conditional expectation ; 
 see ref ( 69 ) . analysis 
 of the adaptive algorithms ( adaptive biasing force or adaptive biasing 
 potential ) with estimates of the conditional expectations based on 
 trajectory averages along a single path are much more complicated . see refs ( 70 ) and ( 71 ) for preliminary results 
 for the wang just like any experimental 
 measurement , free - energy calculations , of either alchemical or geometrical 
 nature , ought to be reported with the associated error bars . in the 
 absence of an error estimate , a free - energy difference is generally 
 of limited utility , making direct comparison with experiment difficult 
 and speculative . although much effort has been devoted in recent years 
 to the characterization of errors that are associated with free - energy 
 calculations , estimating the reliability of such calculations remains an intricate 
 task . this explains why it is not unusual that calculated free energies 
 are still being published without error estimates . ideally , 
 any free - energy difference should be determined from a series of n independent simulations . if this were , indeed , the case , 
 the best possible estimate of the target free - energy difference would 
 be the expected value over the n simulations , i.e. , , where  denotes the estimate of 
 the exact quantity , a , inferred from one individual 
 free - energy calculation . then , the 
 associated mean - square error can be written as20the first term 
 of eq 20 ,  = , is the variance , 
 or the precision of the free - energy calculation . in other words , 
 it is a measure of its statistical error . the second term , , is the square 
 of the bias of the free - energy 
 estimator , i.e. , the square of the difference between the expected 
 value of the estimator and the actual free - energy change , a . the bias , also referred to as the accuracy of the free - energy calculation , is a measure 
 of systematic error of the latter . how to estimate statistical 
 error for the adaptive force method 
 will be discussed in the next section . one important source 
 of bias arises from incomplete sampling in finite - length simulations 
 due to quasi nonergodic behavior of the system . although this behavior 
 can not be quantified , in many instances it can be detected . then , 
 a number of remedial tools are at our disposal . how to recognize and 
 remedy problems with quasi nonergodicity will be discussed in section addressing nonergodicity scenarios . other common 
 sources of bias are inaccurate treatment of intermolecular interactions 
 and algorithmic artifacts arising primarily from imprecise numerical 
 integration of the equations of motion . these contributions to bias , 
 which are common to all simulation methods , will not be discussed 
 here . the adaptive biasing 
 force method is typically applied to nanoscale molecular processes 
 under physiological conditions , a realm dominated by thermal noise . thus , estimating free - energy differences using adaptive biasing force 
 requires careful consideration of statistical error . within the adaptive biasing force framework , free - energy 
 differences are determined by integrating the estimated mean force 
 of the system exerted along the transition coordinate . namely , a free - energy 
 difference on the interval za    zb can be expressed as21where fz is the 
 average of the instantaneous force on the transition coordinate  
 at the position (x ) = z. thus , to determine the statistical error of the free - energy differences , 
 we must delve into the statistics of the mean system force . we assume 
 that the transition coordinate ,  , is discretized and instantaneous 
 forces calculated during the course of the simulation are collected 
 in appropriate bins along . as derived in the supporting information , one way to estimate the error of the 
 mean force in bin i is given by22where f(xt ) = f(xt )  
 fi is the random component of the instantaneous force , ni is the number of samples 
 accrued in bin i , t is the 
 time step of the simulation , and i and f2i are the autocorrelation time and variance of f(xt ) in 
 bin i. given a reliable estimate of the error 
 of the mean force in each 
 bin , we are prepared to analyze how these estimates are propagated 
 to yield free - energy differences . on a discrete grid along the transition 
 coordinate , the integral in eq 21 becomes23where ia and ib are bin indices delimiting the -interval [ za , zb ] . assuming independent behavior in each bin , the error of 
 a sum of mean forces is approximated from from the bienaym 
 formula as equal to the square root of the sum of squares of the errors 
 of these mean forces,24 a notable property of this formula is that the error increases 
 with the size of the interval over which the free - energy difference 
 is calculated . for concreteness , 
 we now consider the error in free - energy differences for the reversible 
 folding of deca - alanine in vacuum . we define the transition coordinate 
 of interest  as the end - to - end distance of the peptide , specifically 
 the distance between the carbonyl carbon atoms of the first and the 
 tenth residue . below we discuss the behavior of the three quantities 
 that enter eq 22 , namely . the number of samples 
 in each bin ni , the standard 
 deviation of the random force ( f2i ) , and the autocorrelation 
 time of this force i . the black curve in figure 4a is the number of samples in each bin of width  = 0.1  for the adaptive biasing force 
 calculated in the range from 12 to 32  ( the red curve will be 
 discussed later in the paper ) . the number of samples , ni , in this range varies from 17  000 
 to 36  000 . thus , as expected from the adaptive biasing force 
 algorithm , the number of samples in each bin approaches uniformity . 
 in this case , nonuniformity of sampling exceeds only slightly a factor 
 of 2 , which corresponds to the variations of the biased free energy 
 not exceeding 0.5 kcal / mol . for comparison , the unbiased free energy 
 changes in the same range by approximately 30 kcal / mol . in figure 4b , we show the distribution of 
 instantaneous forces in four different bins . the forces in each bin 
 are approximately normally distributed , with similar standard deviations 
 of about 20 kcal mol  . this point is underscored in figure 4c , 
 which illustrates that , for deca - alanine in vacuum with langevin dynamics 
 emulating buffeting of the molecule by solvent , the standard deviation 
 of the force acting along the transition coordinate changes very modestly , 
 even though the peptide explores structures that are as disparate 
 as can be imagined . the peptide courses through different compact 
 forms for  < 10  , remains mostly -helical for 
 10 <  < 16  , and forms extended structures with 
 diminishing helical fractions as  increases beyond 16  . the approximate uniformity of f2i seen here 
 is also characteristic of many other systems . for example , it has 
 been previously found that the standard 
 deviation of the instantaneous force on the center of mass of a water 
 molecule is about 2 kcal mol  , irrespective of whether the molecule lies in the bulk aqueous phase 
 or in the hydrophobic core of a lipid bilayer . on the other hand , the transfer of a solute across a liquid  vapor interface is 
 an example of a transition for which f2i is expected to vary considerably 
 with . coordinate dependence of sampling and the system force 
 distribution 
 in reversible folding of deca - alanine . ( a ) samples in each bin for 
 a 10 ns adaptive biasing force calculation on the domain   
 [ 4 , 32 ]  ( black curve ) and a 100 ns calculation on the domain 
   [ 4 , 32 ]  ( red curve ) . note the logarithmic 
 scale on the vertical axis . ( b ) distribution of the -component 
 of the instantaneous system force for different ranges of . 
 ( c ) standard deviation of the -component of the instantaneous 
 system force as a function of . note that f2i is considerable . because this term enters prominently 
 the expression for statistical error in eq 22 , there are significant merits in reducing the dispersion of instantaneous 
 forces . first , because the expression for the ensemble average of instantaneous 
 force is not unique we can , in principle , choose one that reduces 
 f2i . second , variation of forces can be also reduced by a thoughtful 
 choice of the transition coordinate . the third variable in eq 22 , the correlation 
 time of the instantaneous system force i , is also the most difficult to calculate . the sampling in 
 a single bin is rarely sufficient to obtain a converged autocorrelation 
 function of the system force ; thus , in figure 5a we plot the autocorrelation function averaged over 40 bins along 
 different regions of . the correlation time for each region , 
  , is determined by fitting an unscaled exponential function 
 e to the positive values 
 of the autocorrelation function . coordinate 
 dependence of correlations in the system force in reversible 
 folding of deca - alanine . ( a ) autocorrelation functions of the random 
 component of the instantaneous system force for different ranges of 
 . the functions are normalized by the variance f2i ( b ) correlation time for -ranges 
 in panel a. using the calculated 
 values of ni , f2i and i , we now estimate the uncertainties 
 of free - energy differences between deca - alanine structures with different 
 end - to - end distances . the mean system force , with uncertainties calculated 
 by eq 22 , is shown in figure 6a . because only free - energy differences between two points 
 along  , rather than free energies at single points , have a 
 clear physical meaning , we focus on errors associated with these differences . 
 to compute the error in the free - energy difference at points za and zb , we must accumulate the uncertainties of 
 the force between these two points , in accord with eq 24 . if stratification is used and points za and zb are in different windows , then the bienaym formula 
 needs to be used across all windows separating these points . propagation 
 of error in the mean force to the error of free - energy 
 differences . ( a ) mean system force on  for a 10 ns adaptive 
 biasing force calculation , with error bars determined according to 
 eq 22 . we show the estimated 
 error 
 of a between the -helical minimum free - energy 
 state ( za = 14 ) and other 
 states within the interval zb  [ 12 , 32 ]  . if we want to calculate , for instance , 
 the difference in free - energy between the minimum and  = 16 
  , we obtain 3.6  1.8 kcal / mol . larger distances in  
 yield larger error : the free - energy difference between the minimum 
 and the plateau at 25.5  is 19  4 kcal / mol . it is 
 expected that , for sufficiently long total simulation time , t , statistical errors of both the average forces and free - energy 
 differences will decay proportionally as t. the same dependence on t should apply to deviations from nonuniform sampling . if force statistics 
 is collected in m bins , then the root - mean - square 
 deviation from the uniform distribution can be defined as the square 
 root of the variance , var(t),25where nstepk and nstep are the number of samples 
 in bin k and the total number of samples at time t , respectively . for large t , the variance 
 is expected 
 to be proportional to 1/t or , equivalently to 1/nstep . in other words , an 
 example of how var(t ) behaves in a typical simulation 
 is shown in figure 7 . errors that clearly deviate from such 
 behavior are a sign of insufficient sampling . if the problem persists 
 with increasing t , and especially if errors exhibit 
 large fluctuations , then most likely problems with quasi nonergodic 
 behavior have been encountered . var(t ) as a function of the number of molecular dynamics 
 steps , nstep . this quantity was calculated 
 for the permeation 
 of k through a transmembrane hexametric channel of a peptaibol , 
 trichotoxin in a window along the normal to the membrane spanning 
 the range between z equal to 15 and 9 
  . note that for nstep between 0.6  10 and 2.3  10 , nstep  var(t ) differs 
 from its average value in this range by no more than 10% . the inset : 
 the number of force samples along z in the window 
 as a function of nstep . a common manifestation of pathological free - energy calculations , 
 in particular those of geometrical nature , is quasi nonergodicity , 
 wherein sampling along the selected transition coordinate appears 
 to be hampered . here , we inspect closely the effects that impede accurate 
 results to be obtained from free - energy methods , including the adaptive 
 biasing force scheme ( notably hidden barriers in the slow manifolds ) , 
 and discuss how to identify these effects and outline possible remedies , 
 by increasing the dimensionality of the transition coordinate , improved 
 stratification , or sampling aided by multiple - replica strategies . the primary objective of importance - sampling schemes is to facilitate 
 exploration of the transition pathway with a uniform probability . among these schemes , as has been previously emphasized , the adaptive 
 biasing force algorithm uses a local estimate of the gradient , a , acting along the transition coordinate , to erase 
 progressively the original ruggedness of the free - energy landscape . 
 as has already been discussed in section formal 
 convergence of the adaptive biasing force algorithm , this feature 
 is valid from a theoretical standpoint . how true is this in practice ? 
 in most instances , satisfactorily uniform sampling is achieved quite 
 efficiently . occasionally , however , the adaptive biasing force algorithm 
 does not perform as expected . the reminder of this section is devoted 
 to explaining and identifying these special , yet important cases . potential difficulties in applying the adaptive biasing force algorithm 
 are intimately related to the choice of transition coordinate . a basic , 
 yet seldom verified assumption that underlies this choice is the separation 
 between time scales of motions along the transition coordinate and 
 orthogonal degrees of freedom ( see section transition 
 coordinates and rare events ) . for complex , rugged free - energy 
 landscapes , notably those formed by parallel valleys separated by 
 considerable barriers in the direction orthogonal to the transition 
 coordinate ( figure 8) , assuming time scale separation may turn out to be unwarranted . because 
 the adaptive biasing force algorithm exerts no direct action in the 
 orthogonal space , it will not improve sampling at constant . 
 returning to the foundational expression for the adaptive biasing 
 force method , eq 8 , which relates the gradient 
 of the free energy to an ensemble average at constant value of the 
 transition coordinate , the inability to cross hidden barriers in the 
 orthogonal space is tantamount to incomplete ensemble averages and , 
 hence , poor estimates of free - energy changes . common free - energy landscape 
 featuring parallel valleys , collinear 
 to the transition coordinate , . these valleys are separated 
 by substantial free - energy barriers in the direction  , orthogonal 
 to .  can be interpreted as a slow degree of freedom 
 coupled to  and hampering progression along the latter direction . excessively stratified reaction pathways preclude spontaneous crossing 
 of high barriers , typically a1 . wider windows should allow diffusion toward values of  , where 
 the barrier separating valleys in the direction  is smaller , 
 typically a2 . their presence is a guaranteed 
 sign of flawed free - energy calculations , but their absence is not 
 a sufficient condition to ensure that such calculations converged 
 to the correct value . as has been discussed previously , sampling along 
 the transition coordinate in well - behaved simulations should approach 
 uniformity with time , and the statistical error associated with the 
 biasing force or , equivalently , the free - energy differences should 
 decrease in a predictable fashion . if this is not the case , difficulties 
 in equilibrating the system along orthogonal degrees of freedom are , 
 most likely , at play . to illustrate this point , we return to the toy 
 model of deca - alanine reversibly folding in vacuum and to figure 4a . in the range of  between 12 and 32  , in contrast , 
 in the [ 0 , 12 ]  range , sampling remains quite nonuniform , even 
 after 100 ns ( red curve ) . as has already been pointed out in the section transition coordinates and rare events , there 
 are a number of metastable states in this region , all corresponding 
 to similar values of . difficulties in properly averaging over 
 these states markedly impedes equilibration along . in the 
 context of stratified simulations , the presence of hidden barriers 
 along degrees of freedom orthogonal to  often leads to discontinuous 
 biasing forces between adjacent windows . another strategy for 
 exposing apparent nonergodic behavior is to 
 carry out bidirectional calculations , i.e. , initiate the adaptive 
 biasing force simulation from both end points along the transition 
 coordinate . just as in free - energy perturbation calculations , the resulting hysteresis is a good indicator 
 ( although not necessarily a measure ) of error . if the hysteresis markedly 
 exceeds statistical error , the calculated free - energy values are , 
 most likely , poorly converged . in such a case , simply combining data 
 from both directions is not likely to improve accuracy significantly , 
 as their proper weighting remains unknown . when quasi nonergodicity 
 scenarios are encountered , additional 
 simulation strategies , such as multiple walkers or multidimensional 
 transition coordinates , should be brought to bear . formally , the adaptive biasing 
 force algorithm does not prescribe 
 a particular window size nor does it require that consecutive windows 
 along a stratified transition coordinate overlap . provided that convergence 
 has been achieved in each window , the gradient , a( ) , can be reconstructed merely by joining the gradients from 
 individual windows at the boundaries . as we have already argued , from 
 the efficiency point of view it might appear that using small windows 
 is always beneficial . one 
 concern about extensive stratification is that efficient sampling 
 of rugged landscape along orthogonal degrees of freedom may require 
 temporary excursions beyond the window s boundaries , where 
 the barrier separating adjacent valleys is smaller . the barrier between the minima 
 along the orthogonal degree of freedom ,  , may be difficult 
 to cross at  in the range of [ 0 , 0.2 ] , but not in the [ 0.8 , 
 1 ] range . in this case , stratification will create kinetic traps , 
 most likely reducing rather than improving efficiency . a similar 
 problem was observed in simulations aimed at determining the potential 
 of mean force for a small , proline - rich peptide , p41 , bound to the 
 sh3 domain of abl kinase using the root - mean - square deviation ( rmsd ) 
 with respect to the native conformation as the transition coordinate . 
 there , a second minimum was only discovered through creation of a 
 window that did not encompass the native state ( see figure 2a in ref ( 89 ) ) . this second minimum 
 represents a shift in register of the peptide in its binding site , 
 and although it would occur spontaneously given sufficiently long 
 time , its probability can be enhanced by first driving the peptide 
 to higher values of rmsd , which disrupts some of the bonds characteristic 
 of the native state , and then letting it come back . because 
 it is often not possible to predict orthogonal barriers a 
 priori , determining an appropriate windowing scheme typically 
 requires an adaptive procedure . the basic criterion here is the continuity 
 of the biasing force across consecutive windows . if one is concerned 
 about quasi nonergodicity , a good strategy is to start with large 
 windows . if the continuity of forces appears to be satisfactory , one 
 might attempt to improve efficiency through further stratification . the advantage of this strategy is that the approximation to the biasing 
 force acquired from the large - window simulation can be used in smaller 
 windows . if no windowing scheme yields continuous forces , then other 
 strategies , described in the next section , should be employed . simulations involving multiple 
 replicas are perhaps the most powerful strategies to accelerate ergodic 
 sampling along degrees of freedom orthogonal to the transition coordinate . 
 beyond the  embarrassingly parallel  strategy of running 
 independent simulations to obtain more sampling than is possible with 
 a single simulation in the same real time , a number of schemes have 
 been devised that significantly enhance sampling , at the cost of transferring 
 information between replicas . in a prototype example of multiple - replica 
 algorithms , originating with monte carlo simulations , replicas are run at different temperatures and exchanges of system 
 configurations are attempted periodically between replicas , so as 
 to maintain a canonical ensemble at each temperature . the advantage of this method is that replicas at higher 
 temperatures cross energetic barriers more quickly and can pass the 
 resulting configurations to replicas at lower temperatures , preventing 
 the latter from remaining in metastable states . exchange of replicas 
 with different umbrella - sampling potentials , known as hamiltonian 
 exchange , has also been widely successful . with 
 the adaptive biasing force method , each replica multiple - walker strategies range from simply running similar independent 
 adaptive biasing force calculations in parallel , to more complex ones , 
 involving communication between replicas . here we consider two communication strategies that can be used in 
 concert . in the first strategy , which we refer to as shared adaptive 
 biasing force , the instantaneous forces sampled from each walker are 
 collected in a single shared buffer as the simulations progress simultaneously . 
 in the current implementation , the shared buffer is merely conceptual : 
 each walker retains its own buffer , which is synchronized with all 
 the others at fixed intervals . regardless , shared adaptive biasing 
 force can result in significantly faster exploration of the transition 
 coordinate and improved convergence of the free energy , as compared 
 to the case for independent walkers . a second strategy , 
 complementary to the first , is the application of so - called walker 
 selection rules . these selection rules eliminate replicas with values 
 of the transition coordinate that are already well sampled , while 
 duplicating replicas in relatively unexplored regions , enforcing more 
 uniform sampling . selection rules may be implemented using a 
 so - called resampling 
 procedure and weights that are associated with each replica . there are many ways to choose the weights and 
 to implement this idea in practice . the basic requirement is that 
 the selection mechanism automatically vanishes at equilibrium , namely 
 when the biasing force is the mean force . here , we calculate the weight 
 of each walker in a somewhat simpler way than in previous studies , while obtaining similar results . the weight assigned to each walker 
 is the inverse of the number of samples accrued in a neighborhood 
 of bins near the walker . specifically,26where j is the bin occupied 
 by walker i at the time of application of the selection 
 rules , nj(t ) nj(tlast ) is the number of samples that have been 
 accrued in bin j since the last selection and resampling 
 step , and h defines the number of bins surrounding j that are included in the sum . the walkers are resampled 
 on the basis of these weights , and the selection mechanism is switched 
 off when the smallest and largest ni(t ) differ by less than 20% , i.e. , when [ max(1/wti ) we consider 
 reversible unfolding of deca - alanine as a model system for testing 
 different multiple - walker strategies . in figure 9a we compare the results of adaptive biasing force calculations of 
 8 ns total simulated time to a converged free - energy profile . a single , 
 8-ns long simulation gave the poorest results , whereas 16 short , independent 
 calculations perform somewhat better , even though the region 31.432 
  has not been sampled at all . in contrast , shared adaptive biasing 
 force yields complete and more uniform sampling , and a smoother , more 
 reliable free - energy profile . the best free - energy profile , almost 
 indistinguishable from the converged one , was obtained from shared 
 adaptive biasing force combined with walker selection . note that , 
 as expected , the free energies converge to the same values for sufficiently 
 long simulations , irrespective of the multiple - walker strategy , as 
 demonstrated in figure 9b . an implementation 
 of the shared adaptive biasing force algorithm with the selection 
 rules used here is expected to be available in future releases of 
 the molecular dynamics program namd . ( a ) potentials of mean force obtained for different multiple - walker 
 strategies and total simulation times of 8 ns . for reference , the 
 black curve shows the result for a total simulation time of 128 ns . single  refers to a single long 8 ns simulation , and 
  independent  denotes the result of combining force 
 samples for 16 walkers after they ran independently for 0.5 ns . for 
 the curves marked  shared  and  selection  , 
 force samples were synchronized among walkers every 20 ps . also included the application of walker selection , as described in 
 the text , on the same interval . ( b ) potentials of mean force obtained 
 for different multiple - walker strategies and total simulation times 
 of 128 ns . there are several 
 reasons to perform 
 adaptive biasing force calculations with a transition coordinate of 
 dimension greater than one . first , we may be simply interested in 
 how the free energy varies as a function of more than one coordinate . for example , to study the interaction between two molecules , it may 
 be of interest to obtain the free energy as a function of both their 
 distance and relative orientation . as described in section hidden barriers and other challenges to obtaining accurate results , sampling along a single collective variable may be hampered by barriers 
 in orthogonal dimensions . if one can identify these orthogonal dimensions , 
 the application of adaptive biasing force along them will remove the 
 barriers and improve sampling . by itself , obtaining a multidimensional 
 free - energy landscape requires more sampling than is needed to calculate 
 a one - dimensional profile . for example , if one uses n bins along collective variable  and m bins 
 along collective variable  , each component of the gradient 
 over both collective variables , a( , 
  ) , requires mn bins . with one sample obtained 
 per time step , the time required to generate a statistically significant 
 number of samples in all bins will increase approximately proportionally . 
 despite the additional , substantial computational effort , adding a 
 second or third dimension may improve efficiency , as it accelerates 
 convergence and increases uniformity in sampling by allowing the system 
 to more rapidly cross between multiple parallel valleys subsequently , one can recover the potential 
 of mean force along a single dimension through integrating out the 
 additional dimensions . however , when 
 one wishes to include a large number of dimensions , the generalized 
 adaptive biasing force algorithm ( see section the 
 generalized adaptive biasing force algorithm ) may be more appropriate . if the transition coordinate involves angular degrees of freedom 
 that span their full range ( e.g. , [ 0 , 2 ] ) , an additional complication 
 develops . the exact forces and free energies are periodic in these 
 variables , but this is not necessarily true for these quantities burdened 
 with statistical errors . one way to restore the required periodicity 
 is to approximate the free energy as a function of angular variables 
 with spline functions , coefficients of which are computed to distribute 
 errors smoothly across the whole hyperspace spanned by these variables . how to do this has been described by darve et al . as an example 
 application of multidimensional adaptive biasing force , we consider 
 the free - energy landscape of deca - alanine as a function of the root - mean - square 
 displacement ( rmsd ) from three reference structures . figure 10a shows these reference structures : an -helix , 
 a 310-helix , and a particular compact conformation that 
 we refer to as the  conformation due to its visual similarity 
 to this greek letter . in agreement with the one - dimensional calculations , 
 extensively discussed above , the free - energy minimum occurs when the 
 structure is close to the -helix , i.e. , when rmsd = 0.3  . in figure 10b the three - dimensional 
 potential of mean force the violet region on the left 
 corresponds to the -helix , whereas the region on the right 
 is closer to the  conformation than to the other two structures . the minimum in the latter region is only 1.2 kcal / mol higher than 
 the minimum for -helix and lies at ( rmsd rmsd310 rmsd ) = ( 4.3 5.0 2.7 )  . in figure 10c , the multidimensional free - energy 
 landscape thus yields more insight into the correspondence between 
 the free energy and molecular conformations , revealing , in particular , 
 a minimum for a compact structure that was not identified when the 
 one - dimensional transition coordinate was used . the root - mean - square deviation 
 of selected atoms from their positions in each of the three reference 
 structures defines each of the three transition coordinates . ( b ) free - energy 
 isosurfaces as a function of the three transition coordinates . the 
 violet , pink , and gray surfaces contain all points for which the free 
 energy is less than 5 , 10 , and 20 kcal / mol , respectively . in all cases the tick marks represent a root - mean - square deviation of 1  . 
 the minimum value of the free energy , which occurs at ( rmsd rmsd310 rmsd ) = ( 0.3 2.6 2.7 )  , 
 is defined to be zero . ( c ) another view of the surface shown in panel 
 b , with representative structures shown for the two energetically 
 favorable regions . previously , 
 we considered one particular approach to multidimensional adaptive 
 biasing force . here first , it 
 is easy to check that if the dimension of the transition coordinate 
 is larger than one , the biasing force can not be in general derived 
 from a gradient ( it is not conservative ) . on the other hand , this is of course true of the expected long - time 
 limit . a natural idea is , therefore , to project , using , for example , 
 the classical helmholtz hodge projection : change the biasing 
 force ft to t where t = arg min  ft  l . because this is consistent with the expected stationary 
 state , this should not alter the convergence properties . moreover , 
 the interest of this projection is that the variance of the force 
 is reduced , because the nonconservative part is set to zero . second , it would be interesting to develop efficient 
 adaptive biasing force techniques for higher dimensional transition 
 coordinates . the standard implementation of the adaptive biasing force 
 algorithm is currently limited to a dimension up to 3 , because it 
 relies on a cartesian grid of the transition coordinate values , whose 
 complexity is exponential with respect to the dimension of the transition 
 coordinate . it would not be very useful , however , to develop a technique 
 that yields a flat energy landscape along the transition coordinate 
 if the dimension is high , as mapping a high - dimensional cube would 
 be quite time - consuming . one alternative idea , explored in ref ( 101 ) , is to use a bias of 
 the form i=1matii , where 
 ( 1 , ... , m ) denotes 
 an m - dimensional transition coordinate . under appropriate 
 assumptions , one can show the convergence of the method , and preliminary 
 numerical results are encouraging ; see ref ( 101 ) . another route would be to consider a bias 
 of the form i=1matii , or 
 even following greedy algorithms and tensor product approaches used 
 in nonlinear approximation theory ( see , for instance , ref ( 102 ) ) , k1i=1mati , ki , where the functions ( a , ... , a)k1 would be 
 iteratively computed . the objective of this section is to assess the 
 influence of geometrical 
 restraints on the free - energy landscape and how such restraints ought 
 to be treated in the context of adaptive biasing force simulations . 
 toward this end , the nature of the forces at play , either thermodynamic 
 forces or forces arising from external harmonic potentials , will be 
 clarified , and the current strategy for handling the latter , using 
 an extended - lagrangian formalism , will be outlined . one subtlety of the adaptive 
 biasing force algorithm that is often 
 overlooked is its dependence on the measure of the thermodynamic force , 
 which is not necessarily synonymous with the total force acting on 
 the transition coordinate . for most biomolecular simulations , for example , hydrogen - bond lengths are often constrained in a simulation 
 via the rattle algorithm . if only the 
 heavy atom of the bonded pair is involved in the transition coordinate , 
 then the adaptive biasing force algorithm will include the constraint 
 force emanating from the hydrogen in addition to the thermodynamic 
 force , thus contaminating its estimate . a straightforward solution 
 is to include both the hydrogen and its parent atom in the collective 
 variable(s ) defining the transition coordinate , causing the constraint 
 forces to cancel each other . in other words , if the constraint / restraint 
 force is zero in the collective variable s center of mass , 
 it will not contribute to the measured potential of mean force . alternatively , 
 if restraints nonuniformly affect atoms in the collective variable(s ) 
 but are not accounted for in the thermodynamic force used by the adaptive 
 biasing force algorithm , then convergence is impossible ; the adaptive 
 biasing force scheme can not remove forces that it does not measure . though it may seem apparent that one will always want to calculate 
 what the potential of mean force would be in the absence of artificial 
 restraints , there are cases in which externally imposed restraints 
 are meant to be included . ligand 
 binding free energies , which utilizes a staged procedure involving 
 a series of geometrical restraints . by design 
 of the procedure , these restraints must involve the same atoms also 
 being biased at each stage and their contributions are individually 
 determined and tabulated at the end . in this case , to ensure that 
 all necessary forces are included , a procedure such as extended adaptive 
 biasing force is required . as mentioned 
 above , it may be cumbersome to write an analytical expression for 
 the instantaneous force in dimensions larger than one or for complicated 
 transition coordinates . one possible way around this issue is to extend 
 configurational space with a fictitious degree of freedom  
 and define an extended potential,27where k > 0 is a ( large ) 
 force constant that couples  to the transition coordinate , 
 and proceed utilizing the extended transition coordinate , (x, ) =  , instead of . 
 it can be checked that the free energy associated with this extended 
 system is the convolution of the original free energy with a gaussian 
 kernel , a( ) = dz k(z )  a(z ) , where 
 k(z ) = z  exp[(/2k)|z| ] is the gaussian kernel 
 with variance k. the constant k should thus be chosen sufficiently large to ensure that a is a good approximation of a. this enables quick convergence to equilibrium . the adaptive 
 biasing force algorithm can be applied to the extended system , the potential energy being v and the transition coordinate being . notice that the instantaneous force on the extended degree of freedom , 
 corresponding to free energy a , is 
 trivial to compute , because it is equal to the harmonic spring force . the original free energy a can then be recovered 
 by using either deconvolution procedures or simple unbiasing techniques based on formulas such as eq 34 . provided that a good estimate of the biased 
 marginal of the transition 
 coordinate , (z ) , has been collected , 
 the following unbiased estimator of the free - energy derivative can 
 be used,28where z is the conditional average value of  for 
 (x ) = z. here , the biased marginal 
 of  is used as a correction to the inaccurate biasing force . it should be noted that this is valid in a more general context , including 
 the absence of biases ( removing the second term of the right - hand 
 side ) , in which case one recovers the trivial histogram - based estimator 
 of the potential of mean force . although the adaptive biasing force algorithm 
 is often considered 
 merely a free - energy calculation technique , it is , at its core , an 
 importance - sampling scheme that can have applications beyond obtaining 
 free energies . one such application is to calculate kinetic parameters , 
 such as the diffusivity along a transition coordinate , where an importance - sampling 
 scheme is often essential to obtain reliable estimates in all regions 
 of the coordinate , particularly near all - important free - energy barriers . together , the diffusivity and free energy as functions of the transition 
 coordinate can be used to construct a kinetic model of a process of 
 interest , which yields insight beyond the static picture of molecular 
 phenomena given by free - energy calculations . in this section , we describe 
 a recently developed scheme leveraging adaptive biasing force to compute 
 simultaneously the free energy and diffusivity along a transition 
 coordinate . this scheme is just an example , as a number of other approaches 
 in which the adaptive biasing force is used to improve kinetic descriptions 
 of different systems are being investigated . also note that , although 
 this scheme was designed with the adaptive biasing force method in 
 mind , it can be straightforwardly applied to other importance - sampling 
 techniques . determining kinetic parameters can be more subtle 
 than mapping 
 the free - energy landscape , because kinetic descriptions are usually 
 approximations and can not be exactly derived from statistical mechanics . a commonly invoked diffusive model is overdamped 
 langevin dynamics , in which the following equation of motion is assumed for a set of m collective variables z = ( z , z , ... , z),29where ti is the time derivative 
 of the ith collective variable at time t , dij(zt ) is the ( i , j ) component 
 of the diffusivity tensor for the configuration of the collective 
 variables at time t , fj(zt , t ) is the total force on collective variable j for the configuration at time t , j = /z , and i(t ) is a stochastic variable with i(t ) j(t) = 2dij(zt ) (tt ) . 
 the total force on variable j , fj(zt , t ) , is the sum of the biasing force , fjbias(t ) , and the system force , 
 which can be expressed as the negative of the gradient of the potential 
 of mean force fj(z ) = ja(z ) . here , we will focus on diffusive 
 models along a single collective 
 variable , but note that multidimensional transitions have been analyzed 
 by similar approaches . in one dimension , the overdamped langevin model has two free parameters , 
 the free - energy profile a(z ) and 
 the position - dependent diffusivity d(z ) . the adaptive biasing force method , as well as other techniques , 
 provides a route to a(z ) , whereas 
 a number of methods have been applied to determining position - dependent 
 diffusivity . however , some of the most basic methods are not compatible with 
 nonuniform free - energy landscapes , and 
 others are based on the assumption that the free energy can be approximated 
 as a harmonic well . in methods specifically designed for compatibility with the adaptive 
 biasing force method , it is assumed that the dependence of diffusivity 
 on position is weak . to parametrize a 
 diffusive model , a(z ) and d(z ) should be determined consistently , 
 because , for example , coarsening a(z ) results in a reduction of the associated d(z ) . furthermore , to save computational 
 resources , it would be desirable to obtain both the free energy and 
 the diffusivity in the same calculation . however , in many methods 
 cited above , equilibrium or steady - state statistics is assumed , which 
 is likely to yield erroneous results with time - dependent biasing forces . solutions to the problem of calculating diffusivity in simulations 
 with time - dependent biases take advantage of statistical tools such 
 as maximum likelihood or bayesian inference . conceptually , these methods rely on optimizing 
 the parameters of the diffusive model such that the observed trajectory 
 has the greatest likelihood of occurring . the bayesian inference 
 scheme described below begins with assuming a particular dynamical model for the collective 
 variable ( or collective variables ) of interest . in practice , we represent 
 the functions dij(z ) and fi(z ) by piecewise cubic interpolation from a discrete grid in m - dimensional z - space . we seek the optimal parameters comprising the values of the 
 functions dij(z ) and fi(z ) at each grid node that best correspond to the simulated trajectory , 
 denoted by t. for simplicity , we represent these optimal 
 parameters as h*. we consider the trajectory as a set 
 of discrete hops of duration t . given trial 
 parameters h , as well as the biasing forces f(t ) , we compute for each hop p ( { z+1 , t+1}{z , t } , f(t ) , h0 ) , the conditional 
 probability of arriving at transition - coordinate configuration z+1 at time t+1 , given that the system occupied the configuration z at time t. 
 the probability of the complete trajectory given the parameters is 
 the product of the probabilities at each step30this equation yields the 
 probability of the 
 trajectory given an assumed set of parameters . using bayes 
 theorem , we can infer the probability 
 of the parameters given the trajectory:31 an 
 advantage of the bayesian approach 
 is that it permits inclusion of any prior knowledge about the form 
 of the parameters in a consistent way by defining the prior probability pprior(h ) of the parameters . for example , one can assume scale invariance of the function values or smoothness of the functions . finding the optimal parameters 
 then becomes a problem of finding the set of parameters h * that maximizes posterior probability p(h*|t ) , which can be carried out by generating a markov 
 chain of states hk using 
 the metropolis  hastings algorithm . in the one - dimensional 
 case ( or for independent motions along multiple dimensions ) , we can 
 simplify and discretize eq 29 to yield32where gt is a random 
 variable with a standard normal distribution . 
 from the properties of gt , it is evident that eq 30 becomes33below , 
 using this formula and eq 31 , we are able to 
 reconstruct d(z ) and f(z ) for deca - alanine , 
 yielding a complete diffusive model . indeed , we may even set f = 0 and 
 construct the diffusive model from an equilibrium simulation . however , 
 the accuracy of the results of a bayesian scheme are wholly dependent 
 on the quality of sampling near each point along z , necessitating importance sampling for rugged free - energy landscapes . 
 in the supporting information , we discuss 
 assessment of the precision and reliability of the results of the 
 bayesian scheme , which includes estimating the statistical error , 
 as well as checking the consistency of the diffusive model with itself 
 and the output of adaptive biasing force . the z - dependent diffusivity 
 for reversible unfolding of deca - alanine 
 is shown in figure 11b . the variation with z could not have been inferred in any simple way from knowledge 
 of the system force or free - energy profiles . diffusion along z , which can be thought of as the rate at which the end - to - end 
 distance randomly changes , is highest near z = 16 
  , corrsponding to slightly unfolded -helical structures . the ensemble of diverse compact structures in the range z  [ 5 , 10 ]  seems to be associated with the lowest diffusivity 
 values , whereas a secondary minimum appears near z = 22  . given these data , one might hypothesize that the diffusivity 
 is inversely related to the conformational degeneracy at z. note that figure 11b shows considerable 
 dependence of the calculated diffusivity on t , the time over which the trajectory is discretized in eq 33 , which implies that the motion on the observed 
 times is not well modeled by overdamped langevin dynamics , and that 
 a more sophisticated model of motion along z may 
 be needed . ( a ) comparison of the system force as determined 
 by adaptive biasing force and that calculated from the bayesian scheme 
 on the same trajectory . ( b ) diffusivity as a function of the end - to - end 
 distance of deca - alanine for different observation time intervals 
 t . free - energy calculations have become 
 standard tools of statistical 
 mechanics applied to a wide variety of problems in chemistry and biology . this has been possible due to significant theoretical advances in 
 this area , their highly efficient implementations in software packages 
 developed for modern , parallel computers , and remarkable improvements 
 in computational power available for free - energy calculations . yet , 
 these calculations are still not sufficiently mature to be carried 
 out without careful supervision of the end user . this may be fortunate , 
 as there is no substitute for physical insight of the researcher . this also implies that a number of good practices should be followed in free - energy calculations in general , and in 
 applying the adaptive biasing force in particular . these good practices 
 are usually quite simple and involve careful design of simulations , 
 monitoring their progress and postprocessing analysis . in favorable 
 circumstances , they will simply increase confidence that the calculated 
 free energies are reliable . in less favorable cases , they are even 
 more important , as they allow for identifying and correcting shortcomings 
 of the calculations that might adversely affect the results and their 
 interpretation . below , we recapitulate these good practices , which 
 have been already discussed in a considerable detail in the preceding 
 sections.(1)careful 
 choice of the transition coordinate 
 is a key step toward successful free - energy calculations . in making 
 this choice , it is desirable not only to capture the physical nature 
 of a problem of interest but also to ensure small variance and smoothness 
 of the biasing force . ignoring the latter 
 issues may adversely affect efficiency and even correctness of the 
 calculation . reducing the variance should be balanced with the cost 
 of calculating the instantaneous value of the transition coordinate , 
 as usually the former will decrease and the latter will increase with 
 the number of atoms involved . this balance will depend on implementation 
 details and , therefore , some familiarity with the underlying code 
 might be required.(2)a suitable stratification strategy 
 should be planned in advance . if there are concerns about possible 
 quasi nonergodic behavior of the system , and in particular the existence 
 of parallel channels between the end point states , it is recommended 
 to use large windows . if these concerns prove to be unjustified , it 
 is always possible to switch to smaller windows without losing information 
 about the biasing force that has already been accrued . in other instances it is usually more efficient to stratify the transition coordinate 
 into smaller windows throughout a simulation.(3)free - energy calculations should always 
 be accompanied by estimates of statistical errors . for the adaptive 
 biasing force algorithm , a formula to do so exists and should be applied 
 whenever possible . without error estimates the reliability of the 
 calculated free - energy values is questionable and the ability to compare 
 them with experimental measurements is seriously hampered.(4)in stratified simulations , 
 lack of 
 continuity in the biasing force across consecutive windows that clearly 
 exceeds statistical errors is a sure sign of problematic free - energy 
 calculations . this issue should not be ignored , and any attempt to 
 circumvent it by way of some sort of averaging is unlikely to succeed . 
 instead , remedial steps aimed at improving ergodic behavior of the 
 system , such as applying the multiple - walker strategy , should be taken.(5)it is usually very useful 
 to monitor 
 the behavior of the total , biased force , as it should converge to 
 zero in long simulations . moreover , the rate of convergence at sufficiently 
 long times t should become proportional to 1/t . if the biased force along the transition coordinate or 
 within a window exhibits large deviations from nonuniformity , which 
 clearly does not decrease with time , or the convergence rate is erratic , 
 we have , again , a likely indication of quasi nonergodicity in orthogonal 
 degrees of freedom . then , there is no basis for assuming that a reliable 
 free - energy dependence on the transition coordinate can be extracted 
 from the calculated forces . instead , in such circumstances , techniques 
 for removing quasi nonergodicities along orthogonal degrees of freedom 
 should be brought to bear . careful 
 choice of the transition coordinate 
 is a key step toward successful free - energy calculations . in making 
 this choice , it is desirable not only to capture the physical nature 
 of a problem of interest but also to ensure small variance and smoothness 
 of the biasing force . ignoring the latter 
 issues may adversely affect efficiency and even correctness of the 
 calculation . reducing the variance should be balanced with the cost 
 of calculating the instantaneous value of the transition coordinate , 
 as usually the former will decrease and the latter will increase with 
 the number of atoms involved . this balance will depend on implementation 
 details and , therefore , some familiarity with the underlying code 
 might be required . if there are concerns about possible 
 quasi nonergodic behavior of the system , and in particular the existence 
 of parallel channels between the end point states , it is recommended 
 to use large windows . if these concerns prove to be unjustified , it 
 is always possible to switch to smaller windows without losing information 
 about the biasing force that has already been accrued . in other instances it is usually more efficient to stratify the transition coordinate 
 into smaller windows throughout a simulation . free - energy calculations should always 
 be accompanied by estimates of statistical errors . for the adaptive 
 biasing force algorithm , a formula to do so exists and should be applied 
 whenever possible . without error estimates the reliability of the 
 calculated free - energy values is questionable and the ability to compare 
 them with experimental measurements is seriously hampered . in stratified simulations , lack of 
 continuity in the biasing force across consecutive windows that clearly 
 exceeds statistical errors is a sure sign of problematic free - energy 
 calculations . this issue should not be ignored , and any attempt to 
 circumvent it by way of some sort of averaging is unlikely to succeed . instead , remedial steps aimed at improving ergodic behavior of the 
 system , such as applying the multiple - walker strategy , should be taken . it is usually very useful 
 to monitor 
 the behavior of the total , biased force , as it should converge to 
 zero in long simulations . moreover , the rate of convergence at sufficiently 
 long times t should become proportional to 1/t . if the biased force along the transition coordinate or 
 within a window exhibits large deviations from nonuniformity , which 
 clearly does not decrease with time , or the convergence rate is erratic , 
 we have , again , a likely indication of quasi nonergodicity in orthogonal 
 degrees of freedom . then , there is no basis for assuming that a reliable 
 free - energy dependence on the transition coordinate can be extracted 
 from the calculated forces . instead , in such circumstances , techniques 
 for removing quasi nonergodicities along orthogonal degrees of freedom 
 should be brought to bear . following 
 these simple good practices guarantees the improved quality 
 of free - energy calculations carried out by way of the adaptive biasing 
 force method but does not ensure that all problems encountered in 
 such calculations , especially those related to quasi nonergodicity , 
 have been identified . to this end , the toughest challenges to practitioners 
 of the adaptive biasing force algorithm and related methods are linked 
 to multiple slow degrees of freedom leading to orthogonal barriers . one way to address these problems is through a multiple - walker strategy , 
 which has proven effective in a test system that exhibits metastability 
 in the orthogonal space , a significant challenge to the classic adaptive 
 force method . along the same lines , the 
 basic method can be integrated with other enhanced sampling techniques 
 that involve multiple - copy schemes , such as parallel tempering and 
 hamiltonian exchange . another promising 
 research direction is to increase the dimensionality of the bias , 
 as is done in an implicit way in the generalized adaptive biasing 
 force scheme . real - life  
 applications would , however , require an accessible and well - documented 
 implementation , now available as part of the independent collective 
 variables module . although better 
 sampling of orthogonal degrees of freedom is the 
 most promising direction for increasing the efficiency of the adapting 
 biasing force method , there is also room for improvement in achieving 
 convergence along the transition coordinate . for example , it would 
 be of interest to develop adaptive algorithms for the early stages 
 of simulations that would converge faster than the currently used 
 step or ramp functions . one possibility along these lines is to employ 
 a kernel function along . in other words , adaptation in a given 
 bin would depend not only on samples accrued in this bin but also 
 on samples in neighboring bins . another avenue to improve convergence along the transition 
 coordinate , which has not been explored so far , is to exploit the 
 freedom in defining the ensemble average of instantaneous force and 
 identify choices that reduce the variance . for example , the error in integrating force to calculate free energy 
 that is due to binning of force values could be reduced by using improved , 
 smooth interpolation schemes . further , we observe that the best accuracy 
 is obtained not when the free energy is flat but rather when the statistical 
 error of the average force is the same everywhere along the transition 
 coordinate . to achieve this , the number of samples , nk , in bin k should be 
 such that k/(nk ) rather than nk is constant for all bins , 
 where k is the standard deviation 
 of estimated average force in bin k. this is realized 
 if the biasing potential due to the average force used in the standard 
 adaptive biasing force method is supplemented by the term 2kbt ln(k/0 ) , where 0 is the standard 
 deviation at a reference point . in practice , k and 0 are estimated from their running values 
 or their approximations during the course of a simulation . the additional 
 term becomes important whenever the friction coefficient changes markedly 
 with  or parts of the system undergo large fluctuation , for 
 example because they are near phase transition . a number of 
 extensions to the adaptive biasing force method have 
 not been explored yet . for example , the method could be straightforwardly 
 extended to transformations along a parameter of the hamiltonian if 
 an equation of motion was associated with this parameter , as is done 
 in metadynamics . it would be of interest to check whether such application 
 of the adaptive biasing force method were more efficient and reliable 
 than other , related methods currently used for this purpose , such 
 as conventional thermodynamic integration , metadynamics , or the wang another extension of the method that was outlined in one 
 of the original studies on the adaptive biasing force but has not been pursued so far , is to add 
 to the biasing force another contribution that would depend only on 
  and would favor certain states of particular interest to the 
 user or drive the dynamics in a specified direction . the latter would 
 make the adaptive biasing force method a more efficient and better 
 controlled alternative to steered dynamics . finally , we emphasize 
 that the adaptive biasing force is not just 
 a free - energy calculation technique but can also be seen as a general 
 adaptive biasing scheme ( see , for instance , ref ( 124 ) for applications in bayesian 
 statistics ) . indeed , once a correct sampling of the biased measure zp1  exp(pmp/2 )  dp  zt01  exp{[v at0](x ) }  dx has been obtained ( t0 being a fixed time , and the bias being fixed from time t0 : t  t0 , at = at0 , it is easy to recover canonical averages using standard unbiasing 
 procedures:34 although 
 the adaptive biasing force dynamics has been applied so 
 far almost exclusively to chemical and biological systems , it could 
 also be very useful for sampling a multimodal measure ( namely , a 
 probability measure for which high - probability regions , called  modes  
 in this context , are separated by low - probability regions ) in other 
 fields , e.g. , among others , free - energy computations in material sciences or markov chain monte carlo techniques for 
 bayesian inference . in summary , the 
 combination of simplicity , versatility , and strong mathematical underpinnings 
 makes the adaptive biasing force method an attractive target for a 
 wide variety of extensions in statistical mechanics and beyond .