a valuable theoretical insight on mechanisms 
 and rates of chemical 
 reactions can be obtained from ab initio qm / mm free - energy surfaces . one of the main challenges in computing the free - energy surfaces is 
 high computational cost due to multiple evaluations of electronic 
 structures during configurational sampling . to overcome this limitation , a class of computational methods has emerged , which shares the idea 
 of using a coarse - physics reference potential ( rp ) to compute the minimum free - energy reaction path , to which 
 a correction is added to approximate the ab initio qm / mm free - energy 
 surface and(or ) the corresponding activation barrier . while the computational 
 cost can be reduced , by up to two orders with one - dimensional ( 1d ) reaction coordinate , obtaining reliable 
 estimates for the ab initio qm / mm activation free - energy with these 
 approaches perhaps , 
 one of the main reasons is approximating the fine - physics reaction 
 path with the coarse - physics counterpart . the coarse - physics free - energy 
 surface can be substantially different from the ab initio qm / mm free - energy 
 surface . a possible solution to this problem is to reweight the distribution 
 of the reaction coordinate , using , for instance , the umbrella sampling or some alternative approach . however , it was found that obtaining convergent 
 free energies is problematic due to a poor overlap between the rp 
 and the ab initio qm / mm target potential ( tp ) , what leads to the regime 
 where the umbrella sampling becomes inefficient . the second reason is obtaining convergent estimates for 
 the free - energy perturbation ( fep ) , which is often used to estimate 
 the ab initio qm / mm free - energy correction to the coarse - physics reaction 
 path . a popular solution to latter problem involves fixing the reacting 
 fragments in a configuration from a predetermined reaction path , which 
 essentially boils down to evaluation of a higher level - of - theory solvation 
 free energy for a given solute configuration . an alternative approach involves refinement 
 of the rp , which , however , can be a relatively 
 time - consuming operation , since , in general , it is not a trivial automatic 
 procedure . the third approach involves performing a limited sampling 
 with the tp , thus improving the fep estimates 
 by the multistep transformation ( it is important not to confuse the 
 efficiency of a method with the efficiency of the multistep perturbation , 
 involving sampling from , at least , both the fine - physics and the coarse - physics 
 ensembles ) . two latter strategies have been gradually introduced and 
 realized in the paradynamics model which relies 
 on refined empirical valence bond rp and 
 evaluates the ab initio qm / mm correction using a two - step linear response 
 approximation ( lra ) . described in 
 this work is a protocol which is based on the dual 
 sampling ( from both ensembles ) and on the thermodynamic cycle of the 
 paradynamics approach but is novel conceptually and algorithmically . 
 the main idea of this method is to exactly compute the ab initio qm / mm 
 free - energy surface locally , without relying on ( and even computing ) 
 the reference free - energy surface . in other words , the coarse - physics 
 solution is used as the initial guess for the fine - physics solution . here these solutions are not the actual free - energy surfaces but the 
 free - energy penalties of introducing the bias at reactants and at 
 the transition state to the coarse - physics potential and to the fine - physics 
 potential . once the free - energy penalty of biasing the ab initio qm / mm 
 is known , the umbrella sampling is used to provide a reliable estimate 
 of the free - energy surface locally using the biased potential . while 
 computing potential of mean force ( pmf ) with the umbrella sampling 
 method , one has to calculate these penalties with the fine - physics 
 potential in a multistep transformation . here instead , it is computed 
 using a thermodynamic cycle with the coarse - physics potential by introducing 
 the bias to the coarse - physics potential and then switching the force 
 law to the fine - physics potential . another modification is that the 
 actual fine - physic reaction path is located using a low - accuracy fine - physics 
 free - energy surface obtained from the coarse - physics sampling . the paper is organized as follows : in methods , the techniques used in this work for computing the free change 
 of altering the potential are reviewed : the exponential configurational 
 average , the linear response approximation , 
 and the thermodynamic integration . furthermore , the equivalence of 
 the lra method to the tdi method for considered applications is demonstrated 
 analytically and relation to the free - energy interpolation is shown . 
 next , the umbrella sampling method and the weighted - histogram analysis 
 method ( wham ) are reviewed in the section 
 dedicated to computing pmf . finally , these methods are applied to 
 derive the main equations of the presented approach , which is also 
 shown to yield other formulations of the rp - based strategies as private 
 cases while approximating the umbrella sampling equation . the algorithm 
 for solving is formulated as an iterative strategy for improving accuracy 
 of the target free - energy surface at selected regions to the level 
 achieved in computing pmf . in examples , the 
 application of the considered method is provided with a comparison 
 to the pmf results . the 
 free - energy change 
 of moving from one potential to another ( altering the force law ) is 
 calculated using the fep method . it involves evaluating 
 the exponential average of the energy gap between two potentials . for instance , the free - energy of moving from a reference state described 
 with a coarse - physics rp to the state defined by a fine - physics tp , eref  etgt , is1here  = kt ( boltzmann 
 constant multiplied by absolute temperature ) ; 
 e(x ) = etgt  eref is the energy 
 gap ; x is the coordinate vector ; and n is the number of particles the configurational average can 
 be computed over time or over ensemble using molecular dynamics or 
 monte carlo simulations , respectively:2here and further on  ... ei indicates that 
 sampling is performed with the potential ei , while obtaining the helmholtz free energy for the 
 canonical ( nvt ) ensemble . this formulation can be easily generalized 
 to the gibbs free energy ( npt ) and other thermodynamic ensembles . the free - energy change can be computed for the reverse process 
 by taking the configurational average with the tp:3 ideally ( in the limit of infinite sampling 
 or infinitesimal perturbation):4 however , in practice the forward fep and the backward fep 
 of eqs 2 and 3 can show 
 a hysteresis , 
 and the corresponding average is a way to obtain a more reliable estimate . the perturbation 
 in eqs 2 or 3 can be expanded 
 in the power series and truncated to 
 the linear free - energy expansion . while such an expansion can be done 
 for both eqs 2 and 3 , 
 what would lead to the linear response approximation ( lra ) , the computational 
 power of lra was revealed from inspecting free - energy functions of 
 the energy gap:5 it was demonstrated 
 in the context of the various free - energy calculations and for the rp approach in particular that a more reliable estimate of the corresponding free - energy change 
 is obtained from averaging on both end - point potentials ( rp and tp ) . similar conclusion was made in derivation of the acceptance ratio 
 method by bennett , who also noted that this can be seen from the gibbs bogolyubov 
 inequality:6 to further improve the free - energy 
 convergence , the perturbation is computed in a multistep manner by 
 creating a series of n intermediate mapping potentials 
 between the end points potentials eref and etgt:7 this is achieved by changing the perturbation parameter m from 0 to 1 with an increment of . the multistep implementation of eq 7 involves sequentially applying the fep scheme for 
 adjacent simulation windows , so that the average of the forward fep 
 and the backward fep is given by8 alternatively , the same 
 free - energy change can also be estimated by the thermodynamic integration 
 ( tdi ) approach , which , however , requires analytical partial derivative 
 of the free energy with respect to some perturbation parameter .9here , the integral in eq 9 is computed numerically using a trapezoidal numerical integration . substitution of the corresponding derivative of eq 7 with respect to the perturbation parameter to eq 9 yields10 finally , the multistep fep of eq 7 can be computed with the multistep lra approach , which is just the generalization of the lra 
 to the multistep transformation:11note that em+1 em = e , and therefore , 
 the multistep lra and the multistep tdi eqs 10 and 11 are , in fact , identical . in the 
 limit of a big n , all the approaches give the same 
 result ; however , to minimize the computational cost , a small value 
 of n is desired . the questions of statistical efficiency 
 of analyzing data in free - energy calculations with aforementioned 
 and the acceptance - ratio based schemes ( as well as its relation with 
 aforementioned methods ) are discussed in great detail elsewhere . in ref ( 16 ) , it is 
 mathematically shown that in the case of overlapping distribution 
 ( histograms ) of the energy gap the acceptance ratio method ( which 
 is not covered here ) is statistically very efficient , while in the 
 case of a poor overlap between two energy gap distributions , the free - energy 
 interpolation ( the curve - fitting method ) has an advantage . in appendix 1 , the relation of the lra scheme to the 
 curve - fitting method is analytically demonstrated . thus , this relation 
 provides theoretical support for using the lra to compute the free - energy 
 change while moving from rp to tp ( when one expects a poor overlap 
 of distributions ) and additionally provides a possible recipe for 
 improving the lra estimate . the choice of the lra and multistep lra 
 techniques was already validated by comparing to the average fep with 
 the multistep transformation in the study 
 which also showed that the most important part in computing the free - energy 
 perturbation is to include sampling from both ensembles . in order to efficiently 
 sample the elevated regions of the reaction 
 free - energy surface , the potential is modified by introducing a harmonic 
 bias , centered at a particular value of reaction coordinate:12here  designates a chosen 
 reaction 
 coordinate ; m0 is the center of the harmonic bias , and the 
 original potential is eqm . in order to sample along the whole range of the reaction coordinate 
 values , a set of harmonic biases is created by changing the mapping 
 ( coupling ) parameter m incrementally 
 from 0 to 1:13here lb0 and ub0 correspond 
 to the values of the selected  at the lower and the upper boundaries . since the convergence of fep depends on the energy gap between the 
 two potentials , implementation of eq 7 might 
 give a faster convergence when there is a large energy gap between 
 the two potentials . but the multistep 
 approach of eq 13 provides more control in sampling 
 the free - energy surface along the rc . when reaction coordinate 
 is represented by a two - dimensional vector , with 1 and 2 components , the corresponding mapping potentials 
 include two harmonic biases , em = em(1,2):14 the choice of the rc is problem - dependent , 
 common examples include 
 the energy gap rc ( for empirical valence bond and similar potentials ) and nuclear 
 rc ( for molecular - orbital - based potentials ) . the reference state for the 
 fep methods is defined as the original 
 potential , while the target state can be defined as the biased potential . this formulation is used for the mapping potentials em given by eq 12 in eqs 8 , 10 , and 11 to evaluate the relevant free - energy changes of 
 introducing the bias . while applying fep and lra approaches to estimate the free - energy perturbation is similar to the previously considered 
 case , the tdi approach requires calculation of a new partial derivative , 
 since the perturbation parameter now is given by eq 13 . the corresponding derivative for eq 9 ( with the fixed force constant ) is given by15by expanding the lra eq 11 using eq 12 and eq 13 one arrives at16 in appendix 2 , it is shown that with the perturbation parameter 
 given by eq 13 , the multistep lra formulation 
 is equivalent to 
 the result of numerical thermodynamic integration using a trapezoidal 
 rule with the derivative given by eq 15 . as was already mentioned , the multistep formulation of eq 13 provides an effective way of improving the efficiency 
 of sampling along the reaction coordinate , which is crucial in pmf 
 calculations . the efficiency of sampling can be monitored by creating 
 a distribution function ( histograms ) . for the best sampling efficiency , 
 one should try to achieve a uniform distribution , which can be done 
 by adjusting the force constants . however , increasing the force constant 
 will require increasing the number of simulation windows to ensure 
 the even sampling distribution ( histogram overlap ) along the whole 
 range of the reaction coordinate . special care should be taken when 
 the initial system configuration is far from the minimum on the mapping 
 potential , since it might create strong biasing forces . one way to 
 solve this is to gradually increase the force constant while monitoring , 
 in parallel , the deviation of the reaction coordinate from the center 
 of bias . the ultimate 
 goal in the free - energy 
 simulation , however , is to obtain the probability distribution ( and 
 the corresponding pmf ) with unbiased potentials . therefore , it is 
 necessary to recover the corresponding boltzmann probability distribution 
 from the sampling obtained with mapping potentials . the target pmf 
 obtained with a tp is given by17where the probability 
 distribution 
 function of the reaction coordinate is18 in order to improve the sampling efficiency , 
 a bias is introduced to the original potential ( e.g. , see eqs 12 and 14 ) . the unbiased probability 
 distribution of eq 18 is related to the configurational 
 averages computed with the modified ( biased ) target potential given 
 by eqs 12 or 14 , which 
 is further called tm , 
 through the umbrella sampling formula:19 rewriting this equation in 
 terms of the free - energy changes gives 
 two terms:20 the first right - hand side ( rhs ) term ( the denominator of eq 19 ) is the free - energy penalty for introducing a bias 
 ( cf . eq 2 ) , which can be calculated by any fep 
 method with a multistep scheme . the second rhs term is the numerator 
 of eq 19 , and it arises from the probability 
 distribution of the reaction coordinate . combining the result 
 of several mapping potentials when calculating 
 the pmf also , the free - energy functions can be obtained by calculating 
 the weight - average among different mapping potentials , which is called 
 the weighted pmf ( w - pmf):21 within the umbrella sampling formulation of eq 19 , the bias to the original potentials can also be represented 
 by gaussians , evb mapping potential , etc . once the free - energy levels , corresponding to denominators of eq 19 , are known , one can immediately evaluate multiple 
 distributions [ e.g. , two - dimensional ( 2d ) projections from 1d mapping ] . another tool , employed to 
 combine results of simulations is weighted histograms analysis method , 
 wham:2223 wham equations provide both the pmf 
 and the free - energy penalties of introducing the bias . the methodology described 
 in the previous section is a conventional way of calculating the pmf . however , mapping the reaction free - energy surface with a fine - physics 
 tp is extremely computationally expensive , particularly in a many 
 dimensional reaction coordinate space , since one has to perform a 
 long configurational sampling with multiple mapping potentials . while 
 using the coarse - physics sampling significantly reduces the computational 
 cost , the main question becomes how to use this sampling to obtain 
 reliable estimates for the fine - physics model . a straightforward 
 implementation of the rp approach involves performing sampling exclusively 
 with the rp , followed by applying a reweighting strategy such as umbrella 
 sampling or wham - based reweighting strategy . from sampling performed with a biased reference 
 potential given 
 by eq 12 ( which is further called rm ) , the fine - physics distributions are 
 recovered with eq 19 , which becomes24with the use of the free - energy 
 representation 
 eq 20 , the denominator of eq 24 can be represented as two free - energy terms:25 but 
 f(etgt  eref ) does not depend on the bias m0 and the dependence on the reaction coordinate is integrated out , 
 it is a constant free - energy shift , which can be subtracted from the 
 free - energy shifts for all mapping potentials rm , whereas the second term is accurately evaluated 
 from sampling with the rp using the multistep fep . thus , evaluation 
 of the numerator in this approach is the most 
 problematic step due to a poor overlap in the configurational space 
 ( the regime in which the umbrella sampling is inefficient ) . in other 
 words , configurations generated with the rp are , mostly , unlikely 
 configurations for the tp . the corresponding sampling is not representative 
 for the tp , and the pmf convergence is slow , thus providing a low - accuracy 
 estimate . an alternative 
 to computing 
 the low - accuracy target free - energy surface with the reweighting approach 
 of eq 24 is based on computing the vertical 
 free - energy perturbation from the high - accuracy coarse - physics pmf . the free - energy perturbation is calculated with the rp averages only , 
 that is using eq 1 or its linear expansion . this method of computing the perturbation is very similar to a number 
 of other reference - potential based strategies . then the denominator of eq 19 which is the free - energy of introducing the bias 
 to the target potential is computed via the thermodynamic cycle f(etgt  tm ) = f(etgt  eref ) + f(eref  rm ) + f(rm  tm ) , where the first two terms 
 are discussed in the previous section and the last term is computed 
 using the linear expansion , thus27 while the convergence of eq 27 is limited 
 by convergence of the denominator , the accuracy of evaluating the 
 activation free - energy is also limited by the assumption of eq 26 , which with 1d representation of the reaction coordinate 
 is additionally equivalent to the assumption of the same reaction 
 paths . in certain cases , the assumption of similarity 
 of reaction paths can be justified , but it is hard to quantify the 
 error : it might work in some cases and fail in others . the issue is 
 addressed using a parametric or functional fitting of the rp to the 
 tp , what was proposed and implemented in the paradynamics approach 
 for empirical valence bond method ( for 
 semiempirical methods , such refinement can be adapted from other works ) . however , performing parametric refinement 
 of the rp is a nontrivial task and requires a good knowledge of the 
 potential . an alternative strategy involves 
 functional refinement of the rp by adding the difference between the 
 rp and the tp along the reaction coordinate approximated by gaussian 
 functions . in recent works , an implicit refinement 
 strategy is used by choosing a suitable semiempirical rp , what requires 
 a detailed knowledge on the parameterization of the potential . the paradynamics approach assumes that after refinement the reaction 
 paths on the reference and on the target free - energy surfaces are 
 very close , and the error is limited by the denominator in eq 27 . to address this issue , the lra approach 
 is applied to compute the free - energy perturbation while switching 
 from the rp biased at reactants and the ts , m0  ref , to the 
 tp with the same bias , assuming that 
 tgt  ref. then by virtue of eq 26 , the 
 free - energy perturbation reproduces the difference between two free - energy 
 functions : f(rm  tm ) = ftgt ( tgt )  fref ( ref ) . in the case when the transition paths ( transition states ) 
 are dissimilar , 
 ref tgt , reliably computing the activation free - energy [ even 
 if the free - energy perturbation f(rm  tm ) converges ] would generally represent significant challenges , 
 since it might not provide the necessary estimate of moving to the 
 fine - physics path . in this work , another strategy for improving 
 reliability and accuracy 
 of the free - energy estimates , alternative to refining the reference 
 potential , is proposed . it is based on computing the low - accuracy 
 target free - energy surface using eqs 24 and 27 ( both schemes are realized simultaneously ) . these 
 surfaces are used only for the purposes of locating the reactants 
 and the transition state for the subsequent limited sampling with 
 the tp . once the reaction path on the target free - energy surface is 
 located , the fine - physics targeted sampling is performed at those 
 regions . this is the main difference between the presented approach 
 and other reference potential - based approaches , which compute the 
 correction to the reference pmf , thus relying on relevance of the 
 coarse - physics reaction path . it is essential that the coarse - physics 
 sampling is not restricted 
 along the narrow coarse - physics reaction path ( in order to cover the 
 fine - physics path ) . this can be generally achieved by calculating 
 the free - energy surface for the reaction coordinate representation , 
 which dimensionality is capable of capturing both paths . one can see 
 the problem of using a 1d reaction coordinate in movies provided as 
 additional external files ( https://www.dropbox.com/sh/53a4coext3i3ywq/0f9jxztrui ) : the coarse - physics sampling is confined to the coarse - physics 
 path and the fine - physics path being poorly sampled . this difference 
 can be even more extreme when the corresponding reaction paths are 
 more associative and dissociative . in that case , using a 2d representation 
 of the reaction coordinate ( and mapping the free - energy surface in 
 two dimensions using eq 14 instead of eq 12 ) is rather necessary and would result in a better 
 sampling along the fine - physics path and a more reliable estimate 
 of the low - accuracy free - energy surface but at a greater computational 
 cost . in the proposed approach , the accuracy in evaluating the numerator 
 of eq 19 ( tm)]tm is not reduced compared to pmf 
 calculations by taking advantage of the limited fine - physics sampling 
 targeted at the selected regions of the low - accuracy free - energy surface 
 ( computing the low - accuracy free - energy target surface allows identifying 
 the local pmf regions for the targeted sampling ) . however , computing 
 the reaction coordinate distribution from this sampling provides the 
 full accuracy pmf regions only locally . to compute the reaction 
 barrier accurately in pmf calculations , the denominator in eq 19 is computed using a multistep fep procedure ( in 
 other words , using many overlapping histograms in wham ) . this perturbation 
 is associated with the free - energy penalty of introducing the bias 
 to the original tp ( see figure 1 ) . alternatively , 
 this penalty is computed with the thermodynamic cycle identical to 
 the paradynamics model:28 calculating 
 the activation free - energy by positioning local pmf 
 regions relative to the coarse - physics free - energy shifts instead 
 of computing full pmf from fine - physics sampling . pmf calculations 
 involve sampling with a series of biased fine - physics potentials t1 .. tm .. along the reaction coordinate ( cyan circles ) . these mapping 
 points are sequentially connected using the multistep free - energy 
 perturbation f(t1  tm ) , what 
 is depicted by cyan arrows . an alternative approach involves computing 
 pmf locally ( blue dots ) using fewer fine - physics potentials ( blue 
 stars ) . the free - energy difference f(t1  tm ) is instead computed from the thermodynamic cycle , shown by 
 blue and red arrows . red arrows show connecting a sequence of biased 
 coarse - physics potentials r1 .. rm .. using the multistep free - energy 
 perturbation f(r1  rm ) . blue 
 arrows show the free - energy perturbation of switching from the coarse - physics 
 potential to the fine - physics potential ( both having the same bias ) 
 at reactants and at the transition state , f(ri  ti ) . this allows for the positioning 
 of two local pmf regions ( blue dots ) and to determine the activation 
 free - energy ( shown with a black arrow ) . the first rhs term of eq 28 does not 
 depend 
 on the reaction coordinate and is just a constant shift ( same at the 
 reactants and at the ts ) , the second term is accurately calculated 
 while mapping the free - energy surface with the rp . while m0  tgt can be different from the bias , close 
 to the ref , the corresponding free - energy penalties 
 f(eref  rm ) are computed and well - converged 
 for all biases from mapping with the rp . the last term can be calculated 
 by any fep methods presented above , but the actual estimate is taken 
 from lra ( which as shown in appendix 1 is 
 closely related to the free - energy interpolation approach which provides 
 a more reliable estimate when two distributions of the energy gap 
 poorly overlap):29 note , however , that if the biases are different in rm and tm ( e.g. , different force constants ) , the equation 
 should be30 in practice , better 
 estimates of the ts and rs regions of the target 
 free - energy surface are obtained with several mapping potentials favoring 
 sampling of the respective regions ( see figure 1 ) , those estimates can be averaged using the procedure described 
 below for even more accurate estimates . the overall efficiency of 
 this approach therefore is determined by convergence of eq 29 . for a rp retaining adequate physical description 
 of changes in electron density of reacting fragments , the corresponding 
 structural changes in the system upon moving to the tp biased at the 
 same region of the free - energy surface are smaller than changes caused 
 by moving the system on the tp from the reactants to the transition 
 state . in other words , the perturbation of eq 30 should converge in less simulation windows than the perturbation 
 associated with changing the reaction coordinate according to eq 13 . furthermore , the method computational cost relative 
 to the conventional pmf approaches will increase as the reaction coordinate 
 dimensionality increases . to further improve accuracy of positioning the local 
 pmf regions , 
 the multistep transformation of eq 7 is used until the left - hand side of eq 29 convergence . in practice , the multistep lra eq 11 is used since it involves only computing the energy gap . 
 the difference between three - steps and two - steps lra estimates are 
 added to the rhs of eq 28 . once f(rm  tm ) is converged by virtue of 
 eqs 28 and 20 the estimate 
 of the activation free - energy barrier is of the same accuracy as computing 
 eq 20 with the pmf approach . iterative accuracy improvement of the fine - physics free - energy 
 surface computed from coarse - physics sampling . the first step involves 
 computing a low - accuracy free - energy surface from the coarse - physics 
 sampling by reweighting and by free - energy perturbation approaches 
 ( red upper block ) . all ensemble averages are computed with the coarse - physics 
 reference potential , < >r . once the surface 
 is 
 constructed , the regions of interest ( reactants and the transition 
 state ) are identified . in the next step , this allows computing the ensemble 
 averages with fine - physics target potential < >t , 
 which are used to calculate a high - accuracy free - energy surface at 
 those regions ( blue medium block ) . after the second step , the accuracy 
 becomes limited by positioning the local regions . in the third step , 
 a limited sampling is performed in the same regions with a linear 
 combination of the two potentials . with the use of the generated data , the free - energy perturbation is further improved by including the 
 corresponding average for the intermediate ensemble . all methods , reviewed above , are implemented in 
 the pd - m program . to demonstrate application of this method , a computational study 
 is performed on an enzymatic reaction , which is shown in figure 3 . the system consists of an enzyme , haloalkane dehalogenase , 
 with a substrate , 1,2-dichloroethane . the reaction scheme , which describes 
 nucleophilic substitution , the sn2 mechanism , is given 
 in figure 4 . the qm region includes atoms shown 
 in figure 4 , where one of the hydrogen is a 
 link - atom . a semiempirical pm6/mm potential 
 was taken as a reference potential and a hybrid density functional 
 theory b3lyp//6 - 31g*/mm as a target potential . 
 to demonstrate that the proposed approach is capable of delivering 
 the activation free - energy estimate of the same accuracy as the full 
 pmf scheme , next , the activation free - energy 
 was computed using the developed algorithm , outlined in figure 2 . the protein 
 backbone 
 is shown in green ; the catalytic residue , aspartate , and the substrate 
 are shown as sticks . the dotted red line depicts the direction of 
 the nucleophilic attack on the substrate , 1 , 2-dichloroethane . the 
 reaction free - energy surface was defined with a 1d reaction coordinate , 
 represented as a difference between the breaking ( c  cl ) and 
 the forming ( c  o ) bond lengths :  = d(c  cl )  d(c the sampling 
 along the reaction coordinate was performed with a harmonically based 
 pm6/mm potential . totally 64 equidistant mapping potential were used with 
 biases centered at 1.125  m0  2.025 , 
 the force constant value was k = 125 kcal/(mol  ) . the electronic embedding qm / mm protocol implemented in molaris - xg with 
 a file - based interface to mopac2012 ( for 
 pm6 ) and to gaussian09 ( for b3lyp ) . the 
 mm part was represented with a scaas sphere 
 of solvent ( protein and water ) described with enzymix force - field of 22  radius . all free energies 
 were computed using an in house pd - m program , which simultaneously 
 process data using all free - energy perturbation and pmf methods reviewed 
 in corresponding methods . the free - energy 
 changes reported below were computed using the corresponding multistep 
 lra method described above . pmfs were computed using the umbrella 
 sampling approach of eq 19 combined with the 
 multistep lra method and compared with the wham estimates . while 
 carrying out configurational sampling with the rp , the energy gap 
 with the tp ( b3lyp//6 - 31g*/mm ) was computed for 13000 configurations 
 each 2 fs . the number of configurations for which the energy gap is 
 evaluated can be significantly reduced in real simulations ( by at 
 least ten times ) , and this number was taken to obtain the methods 
 best possible estimate . from the generated data the low - accuracy target 
 free - energy surfaces were computed using eqs 24 and 27 in combination with the w - pmf , eq 21 . the resulting low - accuracy target pmfs are shown 
 in figure 5 . from these surfaces , the region 
 of interest ( reactants , transition state , and products ) were located 
 ( shown with gray - shaded areas ) . selecting regions for targeted fine - physics 
 sampling from low - accuracy 
 b3lyp / mm free - energy surfaces computed from sampling with pm6/mm reference 
 potential for the benchmark enzymatic reaction . these surfaces are 
 used for detecting reactants and the transition state ( shown by shaded 
 gray areas ) . the first surface ( shown in green ) is computed by reweighting 
 the reaction coordinate distribution from pm6/mm sampling using the 
 umbrella sampling method . another low - accuracy surface ( shown in black ) 
 is calculated by adding the linear expansion of the free - energy perturbation 
 to the pm6/mm high - accuracy pmf . at the identified 
 regions , a limited sampling using the target potential with the same 
 biases ( see table 1 ) was performed for 10000 
 2 fs md steps ; simultaneously the energy gap was evaluated . the reaction 
 coordinate distribution for the obtained sampling is shown at the 
 bottom in figure 6 . the data was used to compute 
 the local pmfs ( using 35 mapping potentials ) . then the 
 corresponding free - energy perturbations were computed at each bias , 
 lra estimates and the corresponding averaged energy gaps are reported 
 in table 1 . also the difference between the 
 linear expansion and the ( 2-step ) lra are given in the last column . improving 
 the accuracy of probability distributions of the reaction 
 coordinate at selected regions of the reaction path by performing 
 local sampling with b3lyp / mm potential for the benchmark system . the 
 upper plot shows local regions of the free - energy surface computed 
 from the targeted sampling : green and magenta triangles show free - energy 
 penalties of incrementally moving along the reaction coordinate relative 
 to the first ( leftmost ) bias at each region . they are generated with 
 multistep linear response approximation ( m - lra ) and with wham , respectively . red and blue lines are the local pmf computed using the weight - averaged 
 umbrella sampling approach ( red ) and wham ( blue ) , respectively . the 
 lower plot shows biased distributions of the reaction coordinate for 
 the data points used in calculations . these estimates added to the corresponding free - energy 
 penalties 
 computed with the reference potential are also shown in figure 7 . by comparing them to the one - step lra estimates , 
 one can see that the two - step lra estimate improved the accuracy of 
 eq 28 . however , the error with respect to the 
 exact solution varies among different lra estimates , therefore , the 
 average estimate is computed as described below . improving accuracy of 
 the free - energy shifts associated with changing 
 the center of bias m0 while moving along the reaction coordinate . red 
 stars correspond to free energies computed from sampling with biased 
 pm6/mm reference potential ( r ) using a multistep free - energy perturbation . the cyan circles represent the exact searched solution and correspond 
 to free energies computed from sampling with a biased b3lyp / mm target 
 potential ( t ) using a multistep free - energy perturbation . green stars 
 correspond to the pm6/mm free - energy shifts to which the linear expansion 
 of the free - energy perturbation to b3lyp / mm potential is added ( at 
 corresponding biases ) . this perturbation is computed as pm6 ensemble 
 average of the energy gap between b3lyp and pm6 , e . these shifts are used in constructing the low - accuracy b3lyp free - energy 
 surface . after performing a limited fine - physics sampling with b3lyp 
 target potential ( t ) , a more - accuracte estimate for these shifts is 
 computed ( purple  ) using two - step linear response approximation . these estimates are used for positioning the local pmf regions ( see 
 the main text ) . next , the estimates are improved further with a three - step 
 perturbation ( blue  ) , which involves computing the average 
 energy gap with intermediate potentials m = 0.5r + 0.5 t . from local ( targeted ) sampling with n biased target 
 potentials at each identified region , the algorithm computes pmf locally 
 and connects these windows using the multistep lra eq 11 f(tm  tm+i ) . on the other hand lra provides an estimate 
 of f(ri  ti ) 
 at each m+i0 . the lra and local free - energy 
 shifts are however related:31 therefore , one can compute the average free - energy 
 perturbation 
 over n simulation windows at each local region:32and use this estimate to position the local 
 pmf region . computed values are reported in table 2 . the 
 first column ( in angstroms ) the second and third columns are free energies 
 of changing the bias , relative to the smallest ( the most negative ) 
 bias in the group . estimates computed with eq 32 were used to 
 position the transition state and products regions relative to the 
 reactants , as shown in figure 8 . increasing 
 accuracy of the b3lyp activation free - energy barrier 
 computed by positioning the local pmf regions from the targeted sampling 
 and its comparison to the full pmf . full b3lyp / mm pmf ( dashed gray 
 line ) is computed from 64 molecular dynamics b3lyp trajectories . black 
 line is a low - accuracy pmf computed from 64 pm6 trajectories using 
 the averaged energy gap with the b3lyp potential . local b3lyp pmfs are computed 
 from this targeted sampling , which includes 5 trajectories at reactants , 
 5 trajectories at the transition state , and 3 trajectories at products 
 ( blue line ) , but the relative positions of these regions is not known . they are found by computing the free - energy perturbations of switching 
 from pm6 reference potential ( r ) to b3lyp target potential ( t ) , shown 
 by black arrows . the local pmf regions , positioned using a two - step 
 linear response approximation , are shown in red . three - step free - energy 
 perturbation slightly further improves the positioning , and the activation 
 free - energy estimate ( shown in green ) . the three - step estimate involves 
 additional targeted sampling at the same regions with 0.5r + 0.5 t 
 potential . in the next step , a limited sampling 
 at the same regions ( and using 
 the same biases as in the previous part ) was performed with a linear 
 combination of the reference and of the target potential , mm = 0.5rm + 0.5tm . from eq 11 , it follows that to increase 
 the accuracy of the free - energy perturbation , one has to compute the 
 energy gap between the target and the reference potential . thus , this 
 step constitutes a three - step perturbation approach (  
 = 0.5 ) , which is estimated using33 now by estimating the average difference 
 between two - step lra and three - step lra , and their average over the 
 simulation potentials at corresponding local regions , a higher accuracy 
 estimate is obtained ( see table 3 ) . the 3-step 
 lra correction for positioning the transition state pmf region ( relative 
 to the reactants ) is only 0.19 kcal / mol . the updated positioning is 
 shown in figure 8 , from which one can see that 
 the positioned local pmfs essentially coincide with the full pmf . while 
 the barrier for the forward reaction in figure 8 is in excellent agreement with the estimate drawn from the 
 conventional pmf , the reverse reaction barrier is about 3 kcal / mol 
 higher than the pmf estimate . this is believed to be the result of 
 an insufficient amount of simulation windows in pmf calculations for 
 the high reverse barrier ( big spacing in eq 13 ) . to understand the methodological advantage of positioning locally 
 computed pmf regions over computing the fine - physics correction to 
 the coarse - physics reaction path , it is instructive to inspect movies 
 provided as additional external files . file movie1.mov shows two high - accuracy 
 free - energy surfaces , pm6 and b3lyp , and the minimum free - energy reaction 
 paths computed from 1d pmf and projected on a 2d representation of 
 the reaction coordinate . the other two files , movie2.mov and 
 movie3.mov , show low - accuracy b3lyp free - energy surfaces computed 
 from pm6 sampling using the umbrella sampling reweighting of eq 24 and one - step lra correction of eq 27 , correspondingly , for 2d representation of the reaction coordinate . 
 from comparison of these surfaces with the high - accuracy b3lyp surface , 
 one can see that the minimum free - energy reaction path on the low - accuracy 
 surfaces is not the one on the high - accuracy surface . file movie4.mov 
 shows high - accuracy local regions of b3lyp computed from limited b3lyp 
 sampling and positioned using 2-step free - energy perturbation . while the goal of this work is 
 to provide a proof of the proposed 
 concept , and not to study dehalogenase , it is important to make sure 
 the calculations are consistent . this system was studied earlier with 
 other reference potential approaches , and the authors 
 of those works came to different conclusions on the applicability 
 ( due to the accuracy limitation ) of the approach . the b3lyp//6 - 31g*/mm 
 free - energy barrier calculated for the catalyzed reaction in protein 
 with full pmf and with the proposed approach are 8.8 and 8.7 kcal / mol , 
 which are drastically lower than the experimental estimate of 15.3 
 kcal / mol . the applied computational protocol 
 was earlier verified by calculating the reference reaction in water , 
 and the corresponding catalytic effect with pm3/mm , pm6/mm , and b3lyp//6 - 31g*/mm 
 models . the pmf protocol employed here 
 fairly closely reproduced the experimental catalytic effect of 11.7 
 kcal / mol using pm6/mm activation free energies computed from pmf in 
 the reference reaction in water ( 19.0 kcal / mol + the solvent cage 
 effect ) and in the protein ( 10 kcal / mol ) . 
 both potentials ( pm6/mm and b3lyp//6 - 31 g * ) are known to underestimate 
 the activation free - energy barrier in the gas phase , with b3lyp//6 - 31 g * 
 reportedly underestimating by as much as 8 kcal / mol compared to the 
 mp2//6 - 311 g * * results . additional source 
 of error can be in the way charges for qm atoms are derived ; for instance , 
 for a similar type of reaction in water between methyl chloride and 
 chloride , using mulliken charges with pm3/mm potential leads to about 
 3 kcal / mol lower barrier than using the charges fitted to the electrostatic 
 potential . while reproducing the absolute 
 energetics for the reaction barriers is challenging , the relative 
 catalytic effect leads to the error cancellation and , therefore , is 
 considered to be a reliable tool for checking the consistency of calculations 
 for enzymatic reactions . while 
 analyzing the computational cost of elevating the pmf regions 
 compared to building the full pmf , one should note that for the second 
 and third steps of the algorithm , sampling can be performed with a 
 single mapping potential at each region ( even though it would decrease 
 reliability of the estimates ) . thus , in the case of the 1d reaction 
 coordinate , the second iteration ( for two - step lra ) requires propagating 
 a minimum of two ( three if the products region is included ) simulations 
 with the tp , the third ( and subsequent ) iterations add one extra simulation 
 for each region . considering that the full pmf for the 1d case required 
 more than 60 simulations , local pmf elevated with the paradynamics 
 lra approach requires 60/3 , 60/6 , and 60/12 times less fine - physics 
 sampling simulations ( trajectories ) for the two- , three- , and five- 
 multistep fep correction , respectively . when the mapping of the reaction 
 free - energy surface is performed with a two - dimensional reaction coordinate , 
 the number of simulation windows with the conventional pmf approach 
 increases quadratically , and the local pmf with positioning using 
 the multistep lra approach becomes almost 2 orders of magnitude less 
 expensive than the full pmf in terms of the number of required trajectories . thus , the proposed scheme will be particularly computationally efficient 
 for calculating the free - energy surfaces with 2d representation of 
 the reaction coordinate ) . the focus 
 of the current contribution is to improve accuracy of 
 calculating the fine - physics activation free - energy from coarse - physics 
 sampling , up to the level achieved in computing pmf from the fine - physics 
 sampling . the proposed solution is to compute and to position the 
 fine - physics local pmfs using the fine - physics sampling targeted to 
 selected regions identified on the low - accuracy free energy surface 
 computed from extensive exploratory coarse - physics sampling . limited sampling with the target potentials has been used earlier 
 in the paradynamics model and with other lra - based 
 models for improving accuracy of the 
 free - energy perturbation . here the sampling is targeted more precisely 
 using the low - accuracy fine - physics free - energy surface and is also 
 used for computing the distribution of the reaction coordinate . ( 1 ) instead of computing the free - energy perturbation 
 correction to the coarse - physics minimum free - energy reaction path , 
 the proposed method computes the fine - physics free - energy surface 
 locally , and positions these segments relative to the free - energy 
 shifts computed with the coarse - physics potential . ( 2 ) the segments 
 of the fine - physics reaction path are found from a low - accuracy fine - physics 
 free - energy surface . this surface is obtained from sampling with the 
 coarse - physics potential using the umbrella sampling reweighting of 
 eq 24 and the free - energy perturbation approximation 
 of eq 27 . ( 3 ) no explicit parametric or functional 
 refinement of the reference potential ( along the reaction path ) to 
 the target potential is performed . ( 4 ) the distribution of the reaction 
 coordinate is computed with the target potential from the local targeted 
 sampling , which is analogous to computing pmf locally . ( 5 ) locally 
 computed pmf regions are positioned using the averaged lra eq 32 to minimize the error , which limits the accuracy 
 of the approach . ( 6 ) the accuracy of positioning of the local pmf 
 regions is further improved by computing three - step free - energy perturbation . ( which also can be done over several simulation windows locally , see 
 eq 33 ) . these advances allow for reaching 
 accuracy of the full pmf , what 
 is demonstrated analytically from the umbrella sampling equation and 
 is showed numerically on a benchmark reaction by computing the b3lyp 
 pmf locally ( using the sampling with the pm6 potential ) . these 
 advances can be viewed as a recipe for computing the fine - physics 
 solution from the initial conditions generated with the coarse - physics 
 model . the searched solution in fact is not the free - energy surface , 
 but the free - energy penalty of introducing a bias to the target potential . 
 unlike in the pmf calculations , where this penalty is computed in 
 a multistep transformation of eq 13 , this approach 
 computes the penalty through the free - energy perturbation by switching 
 the reference potential to the target potential at the same bias location . the convergence of the perturbation will be reached in fewer multistep 
 transformations than in pmf mapping along the reaction coordinate 
 if the corresponding structural changes are less extensive . a particularly 
 common and illustrative case in chemistry is the associative versus 
 concerted versus dissociative pathways . the constraint on the 1d representation 
 of the reaction coordinate ( similar to one used for the numerical 
 example ) can satisfy both the associative and the dissociative pathways . thus , if the reference potential and the target potentials favor two 
 different pathways , the perturbation of eq 30 would involve extensive structural changes while switching from 
 the rp to the target potential at the same bias location . in the case 
 of the 2d representation , this problem will not be encountered . the 
 same argument applies to the structural changes caused by different 
 coordination numbers , if the metal center is included , or by structural 
 deformations by mechanical stimuli unaccounted for in the reaction 
 coordinate and not fixed by the bias . successful and efficient 
 use of sampling with the rp in calculating 
 the activation free - energy involves overcoming three key challenges : 
 ( 1 ) choosing a good rp , ( 2 ) choosing an efficient strategy for sampling 
 rare events , and ( 3 ) choosing a good reaction coordinate and an efficient 
 scheme for relevant free - energy calculations . the choice of the rp is , perhaps , 
 the most important factor , which determines the overall method computational 
 cost and therefore its efficiency . first , the proximity of the rp 
 to the tp , or their overlap , determine how representative ( efficient ) 
 sampling with the rp for the fine - physics minimum free - energy reaction 
 path , that is how often the most probable configurations of the tp 
 are generated while sampling with the rp . this , in turn , influences 
 the decision for the third problem , since convergence of the free - energy 
 change while moving to the tp greatly depends on its overlap with 
 the rp . second , the difference in computational cost between the tp 
 and the rp is what essentially determines efficiency of the rp approach . thus , the rp should be orders of magnitude less expensive to compute , 
 and yet it should capture enough physics to roughly approximate the 
 tp . the 
 choice is system dependent and is affected both by the simulation 
 system and by the tp . the choice of semiempirical potential as a rp 
 is motivated by their abundance , popularity of the free - energy based 
 studies of chemical reactions in the condensed phase with semiempirical 
 potentials , and with recently emerged scc - dftb potentials . while the choice of the rp is problem - dependent , it should 
 be noted that the method described in this work is problem independent , 
 still the strategies considered above can be applied for two models 
 of the same dimensionality . trivially , the number of degrees of freedom 
 for a higher level of theory model is higher , thus the number of degrees 
 of freedom should be normalized to the coarser model . here the number 
 of degrees of freedom is reduced at each molecular dynamics steps 
 to 3n for the expense of the electronic degrees of freedom which are 
 condensed in the born oppenheimer approximation to the atomic 
 degrees of freedom via energy gradients acting on nuclei in a particular 
 configuration . another example of reduction of dimensionality is the 
 centroid molecular dynamics approach 
 which reduces dimensionality of the ring polymer to its central distribution . an alternative approach to the normalization 
 is the centroid - based quantized classical particle approach , which increases dimensionality of a classical 
 molecular dynamics trajectory to the ring - polymer type approach in 
 order to describe the nuclear quantum tunneling effect . at least in 
 principle , similar approaches can be extended to normalize dimensionality 
 of atomistic simulations to supra - atomistic coarse grained models , what allows for application of the methods presented 
 in the paper to a wide spectrum of problems in addition to chemical 
 reactions . this approach can also be used to improve accuracy 
 of the activation 
 free - energy by moving from a lower level of theory ab initio qm ( e.g. , 
 hf , dft ) to a higher level of theory qm [ e.g. , ccs(d ) , hybrid dft ] 
 what is routinely done in the energy minimization studies . furthermore , 
 it can be used to increase the qm region in qm / mm , up to moving to 
 a full ab initio qm description . while the presented method 
 removes the need of having the same 
 minimum free - energy reaction path , and therefore of careful refinement , 
 the reference potential the choice 
 of the sampling strategy to construct the free - energy 
 surface is identical to choosing the rare event sampling method . the 
 proposed approach is currently formulated for a common strategy of 
 mapping the free - energy surface with a grid of harmonically biased 
 potentials and combining data from multiple ensembles using free - energy 
 perturbation and umbrella sampling approaches . one thing that needs 
 to be considered in this aspect is the force constant for the harmonic 
 bias . in general , the problem can be that the force constant used 
 to keep the reference potential might not be optimal for the target 
 potential . this issue can be resolved by adjusting the force constant 
 on the fly to make sure that the deviation from the bias center is 
 small . this would improve the umbrella sampling efficiency by reducing 
 the gap between the mapping and the target potentials . this idea in fact 
 is used in the metadynamics model , where 
 the free - energy surface is not computed from the probability distribution 
 but is approximated by the cumulative bias . the proposed algorithm 
 can be generalized to include other sampling schemes , for instance 
 the metadynamics approach . generally , the choice 
 of the reaction coordinate should ensure that the perturbation depends 
 approximately equally on all degrees of freedom , that is , that other 
 degrees of freedom ( on which the perturbation depends strongly ) are 
 included in the reaction coordinate , which is kept fixed by the same 
 bias in eq 30 . this condition also ensures good 
 convergence of the free - energy interpolation and of lra ( see appendix 1 ) . the presented 
 algorithm for computing the activation free - energy is formulated in 
 order to construct a high - accuracy fine - physics free - energy surface 
 locally at selected regions of the reaction path . it is derived and 
 is shown to reproduce the full pmf computed with fine - physics sampling . these regions are positioned using the multistep lra approach , which 
 is chosen based on the analysis in appendices 1 and 2 . the presented algorithm is 
 conceptually similar to another multistage 
 method of computing the target free - energy 
 surface ( and improving its accuracy ) , which computes a low - accuracy 
 ( high - dimensional ) free - energy surface with the tp using the metadynamics approach and then refines the surface to a high - accuracy 
 ( low - dimensional ) free - energy surface using the metadynamics potential 
 as a bias in the umbrella sampling . in this approach , however , both 
 the metadynamics and the umbrella sampling steps require sampling 
 with the target ( fine - physics ) potential , and the improvement of accuracy 
 is achieved by reducing the dimensionality and by increasing the sampling 
 accuracy ( and computing the reaction coordinate distribution ) with 
 the umbrella sampling . a similar idea was also proposed in the paradynamics - based 
 model in the sense that both the reference 
 and the target low - accuracy free - energy surfaces were computed using 
 the energy minimization approach and fitted with gaussians . the difference 
 was additionally used to refine the original reference potential , 
 which was used in the umbrella sampling method ( thus avoiding sampling 
 with the target potential ) , while computing a higher - accuracy free - energy 
 surface estimate . the ideology of this work is to use a coarse - physics 
 sampling to 
 construct a low - accuracy fine - physics free - energy surface , which is 
 used for targeting the fine - physics sampling , computing the fine - physics 
 free - energy surface regions locally and positioning them relative 
 to the coarse - physics free - energy shifts . the use of the coarse - physics 
 reference potential sampling is thus the key to reducing the computational 
 cost of calculating the fine - physics activation barrier or the reaction 
 path ( in many dimensions ) . obtaining accurate fine - physics activation 
 free energies is achieved by , first , identifying the actual fine - physics 
 minimum free - energy path on a low - accuracy fine - physics free - energy 
 surface ( which is computed from the coarse - physics sampling ) and by , 
 second , refining accuracy of selected regions of the reaction path 
 ( rather than by correcting the coarse - physics path ) using targeted 
 fine - physics sampling at those regions .
proposed 
 in this contribution is a protocol for calculating fine - physics 
 ( e.g. , ab initio qm / mm ) free - energy surfaces at a high level of accuracy 
 locally ( e.g. , only at reactants and at the transition state for computing 
 the activation barrier ) from targeted fine - physics sampling and extensive 
 exploratory coarse - physics sampling . the full free - energy surface 
 is still computed but at a lower level of accuracy from coarse - physics 
 sampling . 
 the method is analytically derived in terms of the umbrella 
 sampling and the free - energy perturbation methods which are combined 
 with the thermodynamic cycle and the targeted sampling strategy of 
 the paradynamics approach . 
 the algorithm starts by computing low - accuracy 
 fine - physics free - energy surfaces from the coarse - physics sampling 
 in order to identify the reaction path and to select regions for targeted 
 sampling . 
 thus , the algorithm does not rely on the coarse - physics 
 minimum free - energy reaction path . 
 next , segments of high - accuracy 
 free - energy surface are computed locally at selected regions from 
 the targeted fine - physics sampling and are positioned relative to 
 the coarse - physics free - energy shifts . 
 the positioning is done by 
 averaging the free - energy perturbations computed with multistep linear 
 response approximation method . 
 this method is analytically shown to 
 provide results of the thermodynamic integration and the free - energy 
 interpolation methods , while being extremely simple in implementation . 
 incorporating the metadynamics sampling to the algorithm is also briefly 
 outlined . 
 the application is demonstrated by calculating the b3lyp//6 - 31g*/mm 
 free - energy barrier for an enzymatic reaction using a semiempirical 
 pm6/mm reference potential . 
 these modifications allow computing the 
 activation free energies at a significantly reduced computational 
 cost but at the same level of accuracy compared to computing full 
 potential of mean force .