the rational design of 
 high - affinity ligands , in applications like 
 drug discovery and supramolecular chemistry , depends on an understanding 
 of how interactions involving the ligand , receptor , and solvent contribute 
 to the overall observed affinity .
calorimetric studies , which break 
 down the binding affinity into enthalpic and entropic components , 
 have been helpful in this respect and have led to discussions of entropy 
molecular simulations offer additional insight , as they link atomistic 
 detail and binding thermodynamics , but , although decompositions of 
 computed binding free energies into enthalpic and entropic terms are 
 sometimes performed , this practice is not common , perhaps because these terms are regarded 
 as being more difficult to compute .
importantly , simulations often 
 allow the thermodynamic components themselves to be further subdivided 
 to provide rich detail regarding the driving forces of the binding 
 interaction , as we recently demonstrated for the binding enthalpies 
 of a series of host  guest systems . here , we extend our previous work to describe what may be 
 termed 
 computational calorimetry , in which both a binding free energy and 
 a binding enthalpy are obtained self - consistently from the same set 
 of simulations .
we focus on host  guest systems ( figure 1 ) for which ample experimental 
 data are available and extensive sampling can be achieved via gpu - enabled simulations . for these systems ,
just a few days of computing 
 on commodity hardware yields numerical precision as good as or better 
 than the corresponding reported experimental uncertainties .
this precision 
 is enabled , in part , by enhancements to the attach  pull 
release 
 ( apr ) free energy framework , which include the addition of orientational 
 and conformational restraints that facilitate efficient simulations . 

as the centerpiece of this study , we validate the self - consistency 
 of this approach by demonstrating excellent agreement between binding 
 enthalpies computed directly , via end - state mean potential energies , and indirectly , 
 via a vant hoff fit of the binding free energy .
furthermore , we compare the numerical 
 efficiency of these direct versus vant hoff approaches , determine 
 whether hydrogen mass repartitioning ( hmr ) simulations can be used to speed the calculations without impacting their fidelity , 
 and compare the results from two widely used approaches to estimate 
 the statistical uncertainty of mean quantities from simulations with 
 correlated data sets : statistical inefficiency analysis and blocking analysis .
structures 
 of the host ( left ) and guest ( right ) molecules studied 
 in this work .
the protonation state used in this study is shown for 
 each guest and reflects the dominant species under experimental conditions . 

binding thermodynamic values were calculated for cucurbituril ( cb7 ) 
 with all guests except hex and for -cyclodextrin ( cd ) 
 with hex .
the extensive sampling and uncertainty analysis we perform , coupled 
 with the relative simplicity of host guest systems , rules out sampling 
 error as a major contributor to the deviation .
rather , to build rational 
 explanations of binding , further testing and improvement of simulation 
 force fields will be necessary
. the computational calorimetry framework 
 we describe will be helpful in performing such validation and development 
 studies .
the model host  guest 
 systems we chose for this work , cucurbituril ( cb7 ) and -cyclodextrin 
 ( cd ) , exemplify key features of protein  ligand binding , 
 but they are far more computationally tractable due to their small 
 size .
cb7 systems feature high binding affinities and thus require a rigorous approach 
 to deal with the steep energy landscape encountered when removing 
 the guest .
additionally , binding enthalpies have been previously calculated for these systems using the same force field 
 parameters but with a different simulation box setup , thus providing 
 a convenient comparison to study the influence of the choice in setup , 
 as detailed below .
although their 
 binding affinities tend to be more modest than those of cb7 systems , 
 their flexibility and large number of hydrogen - bond donors and acceptors 
 offer a potentially more realistic representation of protein 
a complicating issue with cyclodextrins is that many guest 
 molecules can bind in two orientations that do not exchange on the 
 microsecond time scale of the present simulations .
for example , hexanoate 
 can bind to cd with the carboxylate group near the primary 
 hydroxyl face ( primary binding mode ) or at the secondary hydroxyl 
 face ( secondary binding mode ) .
thus , binding calculations must be 
 performed for both orientations and combined to yield a single value 
 ( see section 2.2.3 ) . finally , we study an extremely simple model of binding , an ion - pair , 
 to validate our computational procedures at high precision .
each set of simulations is named according to the 
 binding pair ; the subscript temp indicates a temperature series used 
 for vant hoff binding enthalpy calculations , and the subscript 
 hmr indicates calculations using hydrogen mass repartitioning .
refer 
 to figure 1 and table 1 for illustrative 
 and technical summaries , respectively , of the simulation sets discussed 
 below .
simulation setup types are as following : 
 i = rectangular box with a 3 anchor particle restraint setup ; ii = 
 truncated octahedron box with 1 distance restraint ; iii = rectangular 
 box with a two anchor particle restraint setup .
aggregate simulation time over the 
 entire simulation set is indicated for direct calculations of the 
 binding free energy and binding enthalpy .
note that the binding free 
 energy total time includes the binding enthalpy simulation time since 
 the latter can be included with the former .
free energies and enthalpies were 
 determined for the guest in two orientations inside the cd 
 host , requiring twice as many simulations .
additionally , the total 
 time listed for the h includes the two free 
 energy calculations required to appropriately weight the contributions 
 from the two binding orientations .
calculations were performed only 
 on the primary orientation , as these simulations served as a consistency 
 check .
all cb7 simulation sets used the same 
 parameters as those described 
 by fenley et al . in order to maintain 
 consistency with that publication .
briefly , parameters for the cb7 
 host and the guests were assigned by the program antechamber : partial charges used the am1-bcc method , whereas bonded and nonbonded parameters came from gaff .
all simulations ( both bound and unbound ) were 
 solvated in an orthorhombic box of 2210 tip3p waters with dimensions 
 of approximately 36  36  52  , with the 
 exception of the larger b11 guest , for which 2500 waters were used 
 in a box of dimensions 36  36  60  . 

the simulation with the mvn guest also included one molecule of protonated 
 tris buffer in order to match experiment . included in each simulation 
 box
were three noninteracting anchor particles , which , in conjunction 
 with several restraints , aligned the host 
hydrogen mass repartitioning was tested 
 for the association of cb7 with b2 . for these cb7-b2hmr simulations ,
the parameter / topology file was modified with the parmed.py 
 program ( distributed with amber ) to reallocate the solute atom masses 
 in accordance with previously described hydrogen mass repartitioning 
 schemes .
cyclodextrin parameters were taken from czard 
 et al . , and hexanoate parameters used 
 resp partial charges from the r.e.d . 
 server program and other bonded / nonbonded parameters 
 from gaff . in an attempt to closely match 
 typical experiments , which include 1020 mm phosphate buffer , 
 we included two molecules each of h2po4 and hpo3 and neutralized 
 with an appropriate number of sodium ions .
to our knowledge , there 
 are no established parameters for inorganic phosphate ions in the 
 amber community , so we began with parameters for existing bio - organic 
 phosphates .
partial charges for each species were determined with 
 a resp fit while lennard - jones terms and bonded parameters were taken 
 from gaff . because we observed nonphysical aggregation of phosphate 
 and sodium ions in the initial simulations , we increased the lennard - jones 
 radius parameters for the phosphate oxygen atoms from their default 
 values by 0.29 and 0.33  for h2po4 and hpo3 , respectively , 
 in accordance with conclusions from steinbrecher et al .
this change eliminated the aggregation , but 
 we do not regard these parameters as optimal , and we note the caution 
 from steinbrecher et al . in regards to the uncertainty of their suitability . 

the entire system was solvated with 2210 tip3p waters in an orthorhombic 
 box whose dimensions were approximately 36  36  52  .
as for the cb7 simulations , anchor particles were used to 
 orient the system along the long axis of the box .
calculations were 
 done with ( cd - hexhmr ) and without ( cd - hextemp , cd - hexalt ) hydrogen mass repartitioning 
 of the solute . as a toy system 
 that would allow unambiguous numerical convergence , and thus rigorously 
 test our computational setups
we used a single k and cl ion pair , with cation parameters from qvist and anion parameters from dang , but we modified the ion parameters to create a diatomic 
 system with favorable enough affinity to bind together during simulation . 

in particular , although this combination is known to produce unphysical 
 ion pairing , the binding was not sufficiently 
 tight to mimic a host  guest - like system , so we increased the 
 binding affinity by artificially changing the charges on each ion 
 from 1.0 to 1.3 .
the k - clsimple simulation 
 set , which did not include orienting anchor particles ( see section 2.5 ) , was solvated 
 with 2000 tip3p waters in an isotropic truncated octahedron .
in contrast , 
 the k - cltemp simulations , which did use orienting anchor 
 particles , were solvated with 1000 tip3p waters in an orthorhombic 
 box with dimensions of approximately 25  25  50  , where the long dimension of the box corresponds to the guest 
 pulling axis .
our general approach to free energy calculations extends the theoretical 
 framework described by velez - vega and gilson , in which the binding free energy is computed in terms of the sum 
 of the work , w , required to attach restraints to 
 the complex , pull the ligand from the binding site and then release 
 the attached restraints ( see schematic summary in the left - hand panel 
 of figure 2):1here , rather 
 than performing a single dynamic 
 pulling calculation , we break the calculation into a series of independent 
 simulation windows .
the collection of independent simulations can 
 still be viewed as spanning an attach  pull  release ( apr ) 
 process , although , due to their independent nature , the process could 
 be equally viewed as unrelease 
has two advantages : 
 ( 1 ) because each simulation is independently built and equilibrated , 
 little or no memory effect is encountered from one point on the path 
 to the next ; ( 2 ) the entire calculation can be distributed across 
 a large , heterogeneous computational cluster , in order to improve 
 turnaround . in the following subsections ,
we provide more detail on 
 the binding free energy calculations performed in this work , starting 
 with a thermodynamic integration ( ti ) approach for calculating the 
 binding free energy from the apr simulation data .
release 
 ( apr ) binding 
 free energy calculation ( left ) and the direct binding enthalpy calculation 
 ( right ) .
the data shown was obtained from the primary orientation 
 of the cd - hextemp simulation set at 300 k. error 
 bars are 95% ci . following the approach detailed by velez - vega 
 et al .
, we first define the chemical 
 potential of the 
 restrained host  guest system in the following manner2where   1/rt , r is the gas constant , t is the 
 absolute temperature , c is the standard concentration 
 ( 1 m or 1 particle/1660  ) , and e(x ) is the potential energy as a function of atomic 
 coordinates , x , defined in the reference frame of 
 the host .
the equation includes two nonphysical summation terms that 
 define harmonic restraint potentials on the rotational and translational 
 coordinates , or pose , of the guest relative to the host , xnp , and on a set of conformational coordinates internal to the 
 host and/or the guest , xnc .
the np / nc number of harmonic restraints 
 are described with force constants , kp / kc ( which implicitly include a factor 
 of 1/2 , according to the convention in amber ) , and target values , x0,np / x0,nc , along with scaling 
 coefficients p/c , which are varied 
 from 0 to 1 in order to turn the restraints on and off in the course 
 of the free energy calculations .
in addition , e(x ) implicitly includes a set of energy terms that restrain 
 the overall translational and rotational coordinates of the host within 
 the lab frame , without influencing its conformational distribution . 

most simply , a single atom of the host may be restrained 
 to a lab - frame position .
more generally , it may be convenient to keep 
 the host positioned in and aligned with the long axis of a rectangular 
 simulation box via restraints tied to a set of artificial particles 
 anchored in the lab frame ( black restraints in the figure s1a example ) . in the following sections
, we will 
 refer to five host  guest states , which depend on the status 
 of the variable restraints : bound , attached , pulled , pulled - free , 
 and unbound .
the bound state has the guest bound to the host but with no added 
 restraints .
the attached state is the same , except that restraints 
 have been turned on . in the pulled state , the restraints are still 
 on , but the equilibrium value of the restraint that sets the host 
guest 
 distance has been adjusted to keep the guest far from the host . in 
 the pulled - free state ,
the conformational restraints have been turned 
 off , but the pose restraints are still in force . in the unbound state , 
 the pose restraints are removed and the guest is considered to be 
 at the standard state concentration .
the bound state does not include rotational or translational restraints on the guest or conformational 
 restraints on either the host or the guest .
therefore , using the definition 
 in eq 2 , the bound state 
 corresponds to p = 0 and c = 0 . 
 with p = 1 and c = 1 ,
the chemical 
 potential is defined as the attached state , as it now includes restraints 
 necessary to extract the guest ( depicted in gray in the figure s1b example ) .
the force of attaching these 
 restraints is -dependent and can be found by taking the partial 
 derivative of hg with respect to p and c ( intermediate steps given in supporting information).3note that the angle 
 brackets 
 indicate the mean value from the simulation at a given value of p/c .
also , the scaling of p/c should be performed simultaneously for the above 
 equation to be valid . to compute the work
, we estimate the integral 
 over the force curve ( figure 2 , left)4from a series of simulations 
 with discrete  steps from 0 to 1 , as detailed in the appendix .
the geometry of 
 the host molecule and our restraint setup ( figure s1b ) allows us to vary the target value of a single distance 
 restraint parameter , designated x0,1p in our general equations below , 
 in order to pull the guest away from the host while the other restraints 
 remain constant . at this point ,
all of the restraints are fully turned 
 on ( p = c = 1 ) , and we compute 
 the work of increasing x0,1p from its initial to its final value , 
 which is the pulled state ( figures 3 and s1c ) .
first , we write the mean force as the partial derivative of the 
 chemical potential with respect to the varying parameter , x0,1p ( see intermediate steps in supporting information):5note that , 
 during the pulling 
 phase , usually x0,1p > x1p , so the 
 present expression 
 will generally yield a positive work for pulling the guest from the 
 host , as expected for favorable binding . to find the work , we integrate 
 the force over the change in equilibrium restraint length from the 
 initial attached distance to a final pulled distance at which the 
 host and guest can be considered to be noninteracting.6as in the attachment phase , 
 in practice the integration is estimated from simulations at discrete 
 steps in x0,1p .
restraint scheme demonstrating the use of anchor 
 particles for 
 attach  pull  release ( apr ) free energy calculations . 

p13 are anchor particles , h13 are host atoms , and 
 g12 are guest atoms .
the ed1 , ea12 , and et13 labels indicate distance , angle , and 
 torsion restraints , respectively , which modify only the host translational 
 and rotational degrees of freedom and are included implicitly in the 
 potential energy function , e(x ) , 
 as described in the main text .
these restraints are held constant 
 throughout all simulation windows and do not perturb the host s 
 conformational degrees of freedom .
indicated in gray with labels x1p , x2p , x3p , and xnc are restraints 
 that are attached over a series of simulation windows and are subsequently 
 used to pull the guest out of the host cavity . in this example , x1p is a distance restraint , x2p and x3p are angle restraints , 
 and xnc represents host conformational restraints 
 ( of any harmonic form ) . in parentheses ,
the corresponding spherical 
 coordinate or euler angle is indicated for each translational / rotational 
 pose restraint , as defined in the main text .
see figure s1 for a full illustration of the apr process . in the release 
 phase ,
all of the restraints added during the attachment phase must 
 be removed either by numerical integration or analytically .
those 
 restraints that perturbed the conformational distribution of the system , 
 either host or guest , are removed in precisely the same manner as 
 the attachment phase , except that now the guest is held far from the 
 host . in this case , only c is varied , while the 
 scaling factor for the translational and rotational pose of the guest 
 relative to the host , p , is held constant . the -dependent 
 force curve is integrated and multiplied by 1 to find the 
 work of releasing the nc conformational 
 restraints.7 after removal of the conformational 
 restraints , the system is in the pulled - free state ( figure s1d ) .
note that , in cases where conformation - perturbing 
 restraints are not used , there is no distinction between the pulled 
 and pulled - free states .
finally , since the guest s rotational 
 and translational 
 degrees of freedom are still restrained relative to the frame of reference 
 of the host , we calculate the work to release these restraints and 
 place the guest at standard concentration ( the unbound state , depicted 
 in figure s1e ) .
this work is the difference 
 between the chemical potential of the pulled - free state and the sum 
 of the chemical potentials of the separate host and guest at standard 
 concentration:8here , xh and xg represent the coordinates of the host and 
 guest respectively , so that x = ( xh , xg ) .
the cancellation of 
 the energy potential term is dependent on the assumption that the 
 host and guest are far enough apart to be noninteracting so that e(x ) = e(xh ) + e(xg ) .
std term , it is helpful to reconsider the guest s translational 
 degrees of freedom in the spherical coordinate system ( r ,  ,  ) and its orientation in terms of euler angles 
 ( a , b , c ) . here , 
 we rewrite eq 8 with 
 appropriate jacobian terms and , in accordance with our example figures , 
 identify the two translational restraints as r and 
  and one rotational restraint as b:9this equation can be evaluated semianalytically , 
 without the use of molecular simulations . note that it may readily 
 be generalized to a calculation in which restraints are also applied 
 to  , a , and/or c. one may use the same set of dynamics trajectories to compute wattach , wpull , and wrelease  conf values by the mbar estimator 
 approach .
this is done 
 by postprocessing the simulation trajectories to obtain the restraint 
 potential energy of each simulation frame in the restraint potential 
 of each simulation window .
the freely available pymbar implementation 
 of the mbar estimator can then be used 
 to determine the relative free energy along the entire free energy 
 path . for
some host  guest complexes , the guest can adopt multiple , 
 distinct orientations that do not readily convert without unbinding . 

for example , an asymmetric guest may bind to -cyclodextrin 
 head - in or head - out . in such cases ,
the overall binding free energy 
 can be computed from nb separate binding 
 free energies of the nb noninterconverting 
 configurations:10note that eq 10 places greater weight on results of lower free energy , 
 so errors in the component free energies propagate nonlinearly .
therefore , 
 we estimate both the mean and the uncertainty of gbind , all by numerically sampling from a gaussian distribution based on the 
 values of gbind , n and their respective 
 uncertainties ( see appendix ) . the 
 enthalpy calculations in this article are performed in one of two 
 ways , termed the direct and vant hoff methods , as now detailed .
the binding enthalpy 
 is the difference between the partial molar enthalpies of the complex 
 and the free molecules . in the direct approach , this quantity is obtained 
 from differences in mean potential energy between simulations of the 
 bound and free states .
( note that the pressure  volume contribution 
 to the binding enthalpy is negligible for binding in solution . ) previously , 
 fenley et al . implemented this method 
 by a multi - box approach , which uses four separate simulations of the 
 free host , the free guest , the host 
guest complex , and any 
 additional pure solvent needed to exactly balance the stoichiometry 
 of the bound and unbound simulations . here , we instead use a single - box 
 approach , which involves simply subtracting the mean energy of the 
 host  guest system in the pulled - free state from that of the 
 bound complex ( figure 2 , right)11where u is defined as the 
 potential energy without any restraining potentials included .
it should 
 be emphasized that , in these bound and pulled - free states , the restraining 
 potentials have no influence on the internal coordinates of the molecules . 
 in the pulled - free state
, they only hold the host and guest far apart , 
 along an axis aligned with the long axis of the rectangular simulation 
 box .
the two states considered here correspond to simulations in the 
 first ( figure s1a ) and last ( figure s1d ) simulation windows of the apr binding 
 free energy approach .
as shown in results , 
 this single - box method gives excellent agreement with the prior multi - box 
 approach .
when nb noninterchangeable 
 binding orientations exist , as for the present cd systems , 
 where nb = 2 , the binding enthalpy of 
 the various orientations must be weighted by their respective free 
 energies to yield the overall binding enthalpy:12again , this combination step is carried out 
 using resampling from gaussian distributions based on the mean and 
 sem of each input value .
this approach 
 involves using the vant hoff equation or related forms to 
 extract the binding enthalpy and entropy from calculations of the 
 binding free energy ( described above ) at multiple temperatures .
if 
 one makes the approximation that enthalpy and entropy are invariant 
 with temperature , then both values can be obtained by fitting to the 
 familiar integrated form of the equation13via 
 linear regression analysis on a plot of 

ln(keq ) versus 1/t . here , 
 however , we allowed for temperature dependence by explicitly including 
 a fitted change in the heat capacity , cp , using the following equation , where the 
 thermodynamic parameters were fitted with nonlinear curve fitting 
 tools:14here 
 the ref subscripts indicate the reference 
 temperature for which the parameters are fitted .
we used simulations 
 to compute g(ti ) at either five or seven different temperatures , ti , at uniform temperature intervals 
 and with a central temperature ( i = 3 or i = 4 , respectively ) of 300 k ; chose tref = 300 k ; and used the values of g(ti ) at all 
 five or seven temperatures to fit the values of h(tref ) , s(tref ) , and cp(tref ) via nonlinear 
 fitting .
it is worth remarking that this approach yielded a value 
 of h(tref ) that 
 was statistically indistinguishable from that obtained by using eq 13 as mentioned above . 
 on the other hand , a pure finite difference method , in which the entropy 
 or enthalpy is computed from just two free energy values at nearby 
 temperatures , results in very large uncertainties , as may be appreciated 
 by inspection of the data in the top panels of figure 7 .
thus , numerical uncertainty is substantially 
 decreased by using free energies computed at multiple temperatures , 
 as done here .
we use numerical molecular simulations to estimate 
 the mean quantities 
 in the equations above , so the reported binding enthalpies and free 
 energies are estimates associated with statistical uncertainties , 
 in the sense that replicating the simulations with slightly different 
 starting conditions or random number seeds would yield somewhat different 
 means , given finite sampling .
it is possible to estimate the uncertainty 
 of the mean without actually running time - consuming replicates so 
 long as time correlation within the simulation data is accounted for 
 properly .
first , the blocking method iteratively 
 averages a time course data series ( e.g. , the potential energy , in 
 an enthalpy calculation ) into successively larger blocks and calculates 
 the sem of the block averages at each block size . the resulting plot 
 of block size versus sem will , ideally , display a plateau feature 
 in the curve .
the sem value of the plateau region is an error estimate 
 that accounts for correlation in the data .
however , it is not straightforward 
 to automate detection of a plateau , and , in some cases , there is no 
 clear plateau , as detailed in section 3.5 .
therefore , in order to automate the usage 
 of this method , we conservatively chose to use the maximum sem value 
 obtained by the blocking analysis , even if it comes from the noisy 
 region of very large block sizes at the far right of the graph rather 
 than from a plateau region .
second , the statistical inefficiency ( stin ) method involves estimating the autocorrelation 
 function of the data series .
the effective number of uncorrelated 
 data points is then found by dividing the total number of correlated 
 data points by the stin value . the usual formula for the sem of an 
 uncorrelated data series can then be used .
following what has been 
 presented in the literature , we couple the ti free energy calculations 
 with the blocking method and the mbar calculations with the stin method 
 ( although any combination could be used ) .
further description 
 of the error propagation in each value calculated is given in the appendix . in order to perform 
 the binding free energy calculations described in this work
, restraints 
 must be used to guide the host and guest apart along a free energy 
 path . in one case ( k - clsimple )
, we used a single distance 
 restraint between the ions , primarily as a reference benchmark for 
 more complicated but efficient restraint schemes that control the 
 position and orientation of the pulling axis .
we discuss the restraint 
 setup details for each simulation set ( table 1 ) in the following text .
many of the simulations we present 
 include 2 or 3 noninteracting 
 anchor particles , which are positionally restrained to fixed cartesian 
 coordinates .
the anchor particles have zero charge , zero lennard - jones 
 radius and well - depth , and a mass of 220 da .
to these anchor particles 
 we attach harmonic distance , angle , and torsion restraints that position 
 and orient the host and guest .
figures 3 , s1 , and s2 diagram a typical anchor particle restraint setup for a 
 host  guest system and an ion - pair toy model . for the ion - pair 
 simulations ,
we use two restraints per ion : a distance restraint and 
 an angle restraint set to 180. this combination effectively 
 locks all three translational degrees of freedom .
six restraints are 
 required to control the three translational and three rotational degrees 
 of freedom of the polyatomic hosts .
these include one distance restraint , 
 two angle restraints , and three torsional restraints ; similar restraint 
 sequences have been described previously .
it is critical to note that , in order for the binding free energy 
 results to be valid , these six restraints must be designed so that 
 they do not perturb the conformational distribution of the host molecule ; 
 see , e.g. , figure s1a .
these host restraints 
 are present during the entire apr process , from the beginning to the 
 end ; thus , they are not imposed during the attach phase or detached 
 during the release phase . in all cases with anchor particles , except 
 the cd - hexhmr simulation set ,
the host translational 
 and rotational restraints had the following force constants : distance 
 restraints used a 10.0 kcal / mol force constant , 
 whereas angle and torsional restraints used 100 kcal / molrad .
a weaker distance restraint force constant of 5.0 kcal / mol was tested for the cd - hexhmr simulation 
 set , whereas the other force constants were identical , but no difference 
 was observed .
the target equilibrium values of these host restraints 
 were set equal to the initial values of the input structure , which 
 was manually oriented in ucsf chimera . in all simulations involving cb7 and cd
, restraints were added 
 during the attachment phase to adjust and restrain the host conformation 
 in a manner that would facilitate sampling during the pulling phase . 

this approach is conceptually similar to the confine - and - release method 
 that has been described for alchemical free energy calculations of 
 protein  ligand complexes .
the 
 cb7 host , which is a very stiff molecule with narrow entry portals , 
 poses a high energetic barrier to removing guests from its cavity . 

we used 14 distance restraints , seven spanning each of the cavity 
 portals , to enlarge the cavity opening and reduce the energy barrier 
 for guest removal ( figure s3 ) .
the work 
 of adding these restraints is relatively easy to converge , and opening 
 the portal in this manner avoids sampling problems that would otherwise 
 occur during the pulling phase , as the guest pops like a cork from 
 the restrictive portal .
in contrast to the 
 cb7 molecule , cd is very flexible and can get dragged into 
 distorted conformational substates for many nanoseconds as it clings 
 to the guest during the pulling - phase simulations ; the associated 
 fluctuations and conformational trapping can make convergence difficult . 
 to alleviate this problem , 14 torsional restraints ( two per glucopyranoside 
 monomer )
were added to the -cyclodextrin backbone during the 
 attachment phase of the free energy calculation .
the effect of these 
 restraints is to maintain the symmetrical , canonical shape of the 
 cd molecule .
the work of adding these restraints is easy to 
 converge , as the host s conformation does not fluctuate much 
 in the presence of the bound guest .
overall , the present restraint 
 schemes significantly improve the reproducibility of the calculated 
 binding free energies for the systems we have tested .
additionally , 
 the work of releasing such conformational restraints must be included 
 in the release phase , with the guest now far from the binding site . 

we accomplished this by performing the exact same simulations as the 
 attachment phase except the guest is either held at a constant remote 
 distance ( in the case of the cb7-b2temp , cb7-b2hmr , and cd - hexalt simulation sets ) or removed entirely 
 from the simulation ( all other cb7 and cd simulation sets ) . 
 in the latter cases , we added a single simulation window of the guest 
 in the pulled - free state with host conformational restraints released 
 ( figure s1d ) in order to compute the direct 
 enthalpy .
control 
 of the necessary degrees of freedom was achieved with just two or 
 three restraints for the ion - pair and host  guest systems , respectively , 
 due to convenient properties of 180 angle restraints .
we note 
 here that it was necessary to include constant flatwell restraints 
 during the attachment phase for the weaker affinity systems in order 
 to prevent occasional dissociation of the bound complex .
an arbitrary 
 choice must be made about where to place the boundary of the bound 
 state in these cases .
the details of the guest restraints and 
 chosen number of windows for each calculation are given in the supporting information .
some of the details vary 
 from one simulation set to another as we optimized our protocol during 
 the course of this investigation .
however , none of the changes were 
 made due to concerns about the quality of a previously collected simulation 
 set ; rather , we optimized our overly conservative initial approach 
 for the sake of greater computational efficiency .
all simulations 
 were performed using the pmemd , pmemd.mpi , and pmemd.cuda 
 programs from amber 14 .
production simulations 
 were run in the npt ensemble , with temperature control using a langevin 
 thermostat with collision frequency 1.0 
 ps and pressure control provided by the monte 
 carlo barostat for all simulation sets 
 except the cb7-b2temp , which used the berendsen barostat to check consistency with previous work ( note that we have not observed differences 
 between the two barostats ) .
direct space nonbonded interactions were 
 truncated with a 9.0  cutoff , whereas long - range electrostatics 
 were handled with the pme method using default 
 amber settings .
shake constraints were applied 
 to bonds involving hydrogen , and the simulation time step was set 
 to 2 fs , with the exception of hmr simulations , which used a 4 fs 
 time step .
anchor particles , when present , were maintained at the 
 desired coordinates using 100 kcal / mol positional 
 restraints , except for the cb7-all8 and cd - hexhmr sets , which used 50 kcal / mol .
the equilibration 
 process consisted of the following steps : ( 1 ) minimization of the 
 initial system with positional restraints maintained on the anchor 
 particles , host , and guest , ( 2 ) 50 ps of npt simulation at 50 k with 
 the same positional restraints as the minimization , ( 3 ) 1.0 ns of 
 npt simulation with smooth heating from 50 k to the desired final 
 temperature without positional restraints on the host or guest , and 
 ( 4 ) 2.0 ns of npt equilibration at the final temperature without positional 
 restraints .
two strategies were employed to determine the length of 
 the production simulations . for the cb7-b2temp set , we 
 chose a fixed simulation length for each window involved in the apr 
 path .
for all other simulation sets , we chose a variable simulation 
 length for each window that depended on either reaching a specified 
 error value in the forces or reaching a maximum simulation time , whichever 
 came first .
the simulation windows for the bound and pulled - free states 
 were then extended to at least 1 s when the direct enthalpy 
 calculation was required . for the hmr simulations , which have previously 
 been reported ,
amber s parmed.py program was used to edit 
 the standard parameter / topology file in order to increase hydrogen 
 masses to 3.024 da and reduce the mass of neighboring heavy atoms 
 accordingly .
computational 
 calorimetry simulations were performed on nine host  guest systems 
 for which experimental data are available ( table 2 ) , using the attach  pull  release 
 ( apr ) method for the free energy and the direct , single - box method 
 for the binding enthalpy . the calculated binding thermodynamics display 
 a strong tendency to overestimate the measured favorability of the 
 binding free energy , the binding enthalpy , and the binding entropy 
 ( table 2 and figure 4 ) .
the rmsd values 
 of these thermodynamic quantities , relative to experiment , are all 
 greater than 4 kcal / mol , with the largest rmsd , 8 kcal / mol , 
 observed for the binding free energies due to overestimation of both 
 the enthalpic and entropic contributions . on the other hand ,
the computed 
 binding free energy and enthalpy correlate well with experiment , with r values of 0.92  0.09 and 0.87  
 0.17 , which suggests that there may be a systematic cause for the 
 large absolute deviations from experiment .
comparison of the experimental 
 and calculated binding values for 
 the systems in table 2 .
a1 , a2 , b2 , b5 , and b11 were at 298.15 
 k , g8 and g9 were at 300 k , and mvn was at 298 k. data for mvn is an averaged obtained 
 from the published value and two additional measurements via personal 
 communication .
all simulations were performed at 300 k. the calculated 
 binding enthalpy was via the direct method . calculated ts was determined by subtracting 
 the binding enthalpy from the binding free energy , and the uncertainty 
 was determined by addition of the respective errors in quadrature . 

uncertainty values for the thermodynamic data given as 95% ci calculated 
 with the ti / blocking approach .
the uncertainty in the error statistics 
 are also 95% ci , determined by resampling the data with replacement .
the statistical uncertainties 
 of both the computed and the measured 
 thermodynamic quantities are on the order of 0.5 kcal / mol ( table 2 ) . as a consequence , 
 the deviations of the calculations from experiment can not be attributed 
 to inadequate sampling during the molecular dynamics simulations . 

this view is supported by the highly stable convergence graphs shown 
 in figure 5 ; indeed , 
 the binding free energies and enthalpies obtained with 100% of the 
 simulation data in all cases remain within the confidence intervals 
 obtained with only 10% of the data .
finally , inspection of the simulation 
 trajectories suggests that good conformational sampling has been achieved . 

this is demonstrated in figure 6 , where we overlay the time - averaged structure for the bound 
 state of every host 
guest system reported in table 2 . due to averaging of the guest s 
 orientation within the host , the extensive sampling results in a stick - like 
 appearance of the guests inside the host ; the averaging also gives 
 the host a regular , symmetric appearance . in the case of cb7 with 
 the aromatic guests ,
this differs substantially from the elliptical 
 distortion characteristic of a typical instantaneous bound conformation . 

likewise , the time - averaged cd structure displays a uniform 
 conformation for each glucopyranoside monomer , which contrasts with 
 a typical instantaneous binding pose in which each monomer flexibly 
 adopts its own conformation , producing an asymmetrical appearance . 

the thorough convergence and small statistical uncertainties of these 
 results highlight the feasibility of performing numerically precise 
 computational calorimetry studies on host  guest systems .
convergence 
 plot of the binding free energy and binding enthalpy 
 for the systems in table 2 .
cb7 was 
 the host molecule for all guests except hex , which bound to cd . 

convergence is indicated by percentage of total data due to the varying 
 trajectory length of the simulation windows considered in the gbind or hbind calculations .
. time - average structures for the bound state of the host  guest 
 systems presented in table 2 .
( right ) superposition of the cd - hextemp host  guest 
 system at 300 k in both the primary and secondary orientations .
one potential concern with the 
 present single - box direct approach 
 to computing binding enthalpies is that the guest in its pulled - free 
 state might not really be far enough from the host to be considered 
 noninteracting .
we addressed this issue by comparing the present results 
 with those from our prior multi - box direct enthalpy calculations , in which the noninteracting state is represented 
 by two separate simulations , one for the host and one for the guest , 
 so that host 
figure s4 confirms that the present 
 single - box results are indistinguishable from those of the multi - box 
 approach and thus supports our choice of separation distance as being 
 adequate . additionally , to confirm that the calculated thermodynamic 
 values did not depend on our specific choice of atoms for the restraint 
 scheme
, we performed an additional set of calculations , labeled cd - hexalt , in which an alternate set of atoms was chosen as the restraint 
 points .
the results showed excellent agreement with the cd - hextemp values ( see table s2 ) .
the binding 
 enthalpy can be computed directly , i.e. , from potential 
 energies , or via the vant hoff equation , which exploits the 
 temperature dependence of the binding free energy .
since we can now calculate both binding free energies and 
 direct binding enthalpies at high precision , we have the opportunity 
 to test for self - consistency by comparing binding enthalpies obtained 
 directly with those obtained from binding free energy calculations 
 at several temperatures .
we explored this in the cb7-b2temp and cd - hextemp simulation sets , wherein binding 
 free energies were calculated at seven and five temperatures , respectively , 
 centered at 300 k ( figure 7 , top ) .
the consistency between the direct 
 and vant hoff binding enthalpy results , shown in table 3 , supports the validity 
 of both the binding free energies and the direct enthalpies .
it is 
 also of interest that convergence plots of the binding enthalpies 
 from both approaches ( figure 7 , bottom ) show that the vant hoff approach is much 
 more sensitive to extending the simulation time and hence is significantly 
 more difficult to converge than the direct approach .
also , for both 
 systems tested , the final statistical uncertainty is much larger for 
 the vant hoff approach ( table 3 ) and comes at a much higher computational cost in 
 terms of simulation time ( tables s1 and s2 ) .
demonstration of the temperature dependence of the binding free 
 energy ( top ) and the convergence of the vant hoff and direct 
 binding enthalpy calculations at 300 k ( bottom ) .
data is shown for 
 the cb7-b2temp ( left ) and cd - hextemp ( right ) 
 simulation sets .
uncertainty values are 95% ci . careful inspection of table 3 shows that the uncertainty values produced 
 by the blocking 
 method are consistently larger than those from the stin method .
although 
 our automated implementation of the blocking method is expected to 
 produce slightly larger values than an investigator - performed analysis 
 in which the sem would be assigned based on a plateau in the sem as 
 a function of block size , the large difference between the two uncertainty 
 estimates was surprising .
for example , when comparing the direct 
 and vant hoff binding enthalpies of cb7-b2temp , 
 a two - sample t - test for equal means (  = 0.05 ) 
 would not indicate a significant difference based on the blocking 
 uncertainty values , but it would indicate a significant difference 
 when using the stin uncertainties . this concern led us to develop 
 a toy model that would allow an even more precise validation of our 
 computational calorimetry methods ( see section 3.3 ) and a deeper investigation into the uncertainty 
 estimation procedures ( see section 3.5 ) .
although 
 the agreement between the direct and vant hoff binding enthalpies 
 for the two host  guest systems was strongly encouraging , we 
 wanted to test the comparison to even finer precision in order to 
 detect any potential subtle pathologies . using a modified k - cl ion 
 pair , which was tuned to exhibit favorable binding during simulation , 
 we performed several further simulation studies to validate the present 
 approach , as follows .
first , we established that the use of 
 noninteracting anchor particles to align the pulling coordinate did 
 not influence the results .
the k - clsimple simulation set , 
 which used a single distance restraint in a large water box , produced 
 binding values statistically indistinguishable from the anchor particle 
 aligned k - cltemp simulation set ( see table 3 for results , and see figure s2 for the anchor particle approach with ions ) .
additionally , 
 comparison between the k - clsimple and k - cltemp results suggests that use of a smaller rectangular simulation box 
 aligned to the pulling axis does not introduce error into the binding 
 results ( table 3 ) . 

the binding free energies and direct enthalpies showed excellent agreement 
 to high precision , with confidence intervals of less than 0.15 and 
 0.40 kcal / mol , respectively .
our modified k - cl system could be considered 
 particularly challenging due to the rather large charges that we placed 
 on the ions ( 1.3 ) .
thus , it is reassuring that the optimized 
 rectangular box ( with alignment via anchor particles ) did not incorrectly 
 treat the electrostatics of the system .
finally , we also observed 
 excellent agreement between the direct 
 and vant hoff binding enthalpy for this simple test case ( table 3 ) .
the small uncertainties 
 in the vant hoff binding enthalpies , relative to those for 
 the full host  guest systems , can be attributed to the better 
 behaved temperature dependence of the binding free energy ( figure s5 , top ) , and the convergence behavior 
 of both binding enthalpy methods is correspondingly smooth ( figure s5 , bottom ) .
interestingly , we noticed 
 that , although the mbar / stin uncertainties 
 are consistently smaller than the ti / blocking uncertainties for the 
 k - cl data , the ratio between the two values is smaller than the ratio 
 observed for the host 
although a slight overestimation of the 
 error is expected by our implementation of the blocking method , the 
 variability of the ratio between the two uncertainty estimates suggests 
 the difference might be influenced by the nature of the system being 
 simulated . here
, the simplicity of the k - cl systems possibly minimizes 
 the apparent difference in the uncertainty .
taken together , 
 the results of the k - cl simulations strongly suggest 
 that the theory and implementation of the apr computational calorimetry 
 framework is valid and numerically sound . as a consequence
, it seems 
 likely that the modest deviations between the direct and vant 
 hoff binding enthalpies for the host  guest systems arise from 
 statistical , rather than methodological , error .
curiously , however , 
 there is still a systematic tendency for the blocking estimate of 
 error to exceed the stin estimate .
the hmr technique is appealing 
 because it allows one to double the simulation time for a given elapsed 
 wall - clock time , by going from a 2 fs to a 4 fs md time step without 
 inducing numerical instabilities .
it 
 should be noted that , because we are working within the classical 
 approximation of statistical thermodynamics , the quantities that we 
 aim to compute are independent of the atomic masses .
therefore , any 
 error introduced by hmr would derive from changes in the numerical 
 properties of the simulations .
in particular , increasing the time 
 step can increase error by increasing the mean spatial excursion the 
 atoms make in one time step ; conversely , increasing the masses of 
 the fastest atoms , i.e. , the hydrogens , reduces their mean velocities 
 and hence their mean excursions , for a given temperature , and thus 
 compensates , at least partially , for the increased time step .
the 
 net effect could be to influence the values that the computed quantities 
 approach , in the limit of a long simulation .
in addition , changing 
 the masses and time step influences the dynamics of the simulations and hence can influence the rate of convergence . 

therefore , we carried out an empirical evaluation of the feasibility 
 of speeding thermodynamic calculations , while maintaining reliable 
 results , by performing several binding free energy and binding enthalpy 
 calculations with hmr and a 4 fs time step
. the results for the cb7-b2hmr and cd - hexhmr simulations sets , given 
 in table 3 , indicate 
 generally excellent agreement between the standard approach and hmr . 

( we did not carry out vant hoff binding enthalpy measurements , 
 given the excellent agreement between hmr and the traditional approach 
 for the free energy calculations at 300 k. ) the one potential outlier 
 is the comparison of the direct binding enthalpy for cd - hex . 

although a t - test (  = 0.05 ) would not indicate 
 a significant difference between the binding enthalpy means for cd - hextemp and cd - hexhmr when using blocking uncertainties , 
 the stin uncertainties are smaller and thus indicate that the difference 
 in the means is potentially significant . in the following section
, 
 we investigate this difference further , show that the difference is 
 due to random statistical chance , and offer evidence that the stin 
 uncertainty estimate can sometimes be deceptively small .
as 
 noted by others , the speedup from 
 hmr does not always double the effective sampling rate due to alterations 
 in the viscosity of the system .
for example , with a traditional 2 fs time step , 
 the blocking uncertainty values in the cd - hextemp direct binding enthalpy for the primary and secondary guest orientations 
 are 0.56 and 0.52 kcal / mol , respectively ( table s2 ) . with hmr and the same simulation time ( ns ) , and double 
 the number of data points due to reducing the interval between data 
 writing by one - half
( note , however , that the real time expenditure of gpu 
 hours for hmr was roughly 60% of the traditional approach . )
this suggests 
 that hmr may influence the effective sampling rate to some degree , 
 but the result could be due to random chance and thus requires further 
 study to make a definitive conclusion .
the same trend is not always 
 evident for similar comparisons elsewhere in our data , possibly because 
 the effect is clouded by our evolving methodology for performing the 
 calculations .
overall , hmr does seem to offer an effective increase 
 in sampling 
 rate without impacting the fidelity of the calculations , but the effective 
 speed up is not always straightforward to determine . as discussed in 
 prior subsections , the two approaches we employ to 
 estimate uncertainty , blocking analysis and stin , can give very similar 
 estimates for the sem of a time - averaged quantity obtained from md 
 simulations , especially for simple systems , such as our k - cl model 
 ( figure 8 , top ) .
however , 
 we noticed that the two approaches can disagree substantially for 
 more complicated systems , such as a host  guest complex ( figure 8 , bottom ) . in the 
 latter case
, it appears that the stin approach may register relatively 
 fast fluctuations while missing the consequences of slower processes , 
 leading to underestimates of the sem .
typical blocking curves for the potential 
 energy of the bound state 
 for a k - cltemp simulation ( top ) and a cd - hextemp simulation ( bottom ) .
error bars are an estimate of the 
 uncertainty in the sem ; see flyvbjerg et al . for further details .
the sem predicted by the statistical inefficiency 
 method ( stin ) is indicated with a red line . in the present context , it was of particular concern that 
 comparison 
 of the direct binding enthalpy values obtained from the cd - hextemp and cd - hexhmr simulation sets ( table 3 ) shows a discrepancy 
 in the uncertainty values , which could impact the conclusions drawn 
 about the difference in the means , as noted above .
we traced the primary 
 difference between the two binding enthalpy values to a single source : 
 the average potential energy value of the bound state in the secondary 
 orientation .
there were no structural indications to suggest that 
 the sampling had been biased , so we reran one replicate of this particular 
 simulation for each simulation set .
the mean of the cd - hexhmr replicate was less than 0.1 kcal / mol from the original 
 value , whereas the cd - hextemp replicate differed 
 from the original by 0.6 kcal / mol , and four additional cd - hextemp replicates differed from the original by 0.30.7 
 kcal / mol .
thus , it appeared that the original cd - hextemp mean value was simply a low probability , but reasonable , value from 
 the tail of the expected normal distribution of mean values .
using 
 the cd - hextemp replicates , we compared the observed 
 sem value with the average blocking sem value and average stin sem 
 value .
the observed sem for the six calculations ( i.e. , the standard 
 deviation of the means for the one original value and five additional 
 replicates ) was 0.24 kcal / mol , the average blocking sem was 0.20 kcal / mol , 
 and the average stin sem was 0.14 kcal / mol .
the results show closer 
 agreement between the blocking sem and the observed sem and suggest 
 that there is no statistical difference between the cd - hextemp and cd - hexhmr direct binding enthalpy 
 values ( according to the t - test ) , whereas the stin 
 uncertainty values would indicate a statistical difference . in order to test more systematically
whether the stin method underestimates 
 the uncertainty in some cases , we developed a python script ( figure s6 ) that generates artificial simulation 
 data containing a short - correlation / large - variance process , mixed 
 with several long - correlation / small - variance processes .
the resulting 
 time - series data ( figure s7 , top ) yield 
 blocking analysis graphs ( figure s7 , bottom ) 
 similar to those from real host  guest data ( figure 8 , bottom ) . by using this method 
 to generate 500 distinct , artificial simulation data sets , we were 
 able to calculate the actual sem as a function of 
 simulation time and use this to test the ability of the blocking and 
 stin methods to estimate this observed sem . in order to generate a 
 summary comparison , we averaged the results from blocking analysis 
 and stin over the same 500 artificial time series ( figure 9 ) . early in these runs ( < 10% 
 of the simulation length )
however , for the remainder of the data collection , the blocking 
 method closely tracks ( and slightly overestimates by 1015% ) 
 the observed sem , whereas the stin method underestimates the sem by 
 roughly 50% .
comparison of the observed sem ( black ) with the average 
 stin sem 
 ( red ) and average blocking sem ( blue ) for a data set consisting of 
 500 artificial simulations , each containing 500  000 correlated 
 data points .
a possible explanation for the underestimate by the stin 
 method 
 can be traced to its approach for computing the integrated correlation 
 time . as detailed previously , the tail 
 of the correlation function is generally dominated by noise , which , 
 if included in the integration , tends to make the computed correlation 
 time inaccurate . following prior practice
, we addressed this issue by truncating the correlation function when 
 it first crosses zero . to check
whether this approach misses useful 
 long - range correlation signals , we performed several tests on the 
 artificial data sets in which the correlation function was not truncated . 
 in each case , however , the estimated uncertainty was even smaller 
 than that obtained with the truncation approach , which itself underestimated 
 the observed uncertainty . therefore , regardless of whether or not 
 truncation is used , the stin method appears to be unreliable for the 
 data sets that we studied in this work .
we suspect that better reliability 
 of the blocking method can be attributed to two features : ( 1 ) the 
 blocking method seems to be sensitive to long - correlation / small - variance 
 signals even when they are obscured by short - correlation / large - variance 
 signals , whereas the stin method seems to recognize only the latter 
 in the same situation ; ( 2 ) our automated implementation does not search 
 for a plateau , as in the canonical approach , but reports the largest 
 estimated sem from the entire blocking curve .
the second feature is 
 likely responsible for occasional overestimation of the sem , but it 
 seems to be the simplest and most conservative method for automation .
examples 
 of binding affinity calculations are numerous in the literature since 
 they are typically the most sought after value in lead optimization 
 studies .
our computational calorimetry approach provides additional 
 insight simply by extending the end point simulation windows of the 
 apr binding free energy calculation and extracting the binding enthalpy 
 from the difference in mean potential energies .
as detailed previously , binding enthalpies obtained by the direct method 
 may also be easily decomposed to provide information on the force 
 field terms and/or structural elements that drive binding . as a larger 
 set of host 
guest analyses is collected , a more detailed understanding 
 of the role each chemical moiety plays can be established .
comparison 
 of the calculated and experimental binding free energies , enthalpies , 
 and entropies shows good correlations but large systematic deviations . 

given the correlations , one might wonder whether the systematic deviations 
 should be regarded as problematic since , at least in the field of 
 computational ligand design , predicting a correct ranking of ligand 
 affinities is nearly as useful as predicting absolute affinities . 

however , it is difficult to escape the conclusion that errors of up 
 to 10 kcal / mol point to significant problems with the force 
 field that require corrective action , especially given the high precision 
 of the results .
moreover , the types of errors observed here could 
 generate serious practical problems in the more complex setting of 
 protein  ligand binding .
in particular , there is no reason to 
 expect that chemically varied ligand moieties interacting with varied 
 subpockets of a protein binding site will be uniformly scaled by whatever 
 factor is generating the systematic bias seen here . as a consequence , 

there is no assurance of obtaining a good correlation with experiment 
 in such complex systems .
ligand 
 systems , including deep binding cavities with the potential to generate 
 structured water , steric barriers to entry and exit , conformational 
 fluctuations with long time scale correlations , and multiple binding 
 poses .
we would argue , therefore , that careful force field optimization 
 to improve the agreement of absolute binding thermodynamics for host  guest 
 systems , in combination with more traditional target data , is a useful 
 path to improve the reliability of protein  ligand modeling 
 for drug discovery .
the demonstration of good agreement 
 between the direct and vant 
 hoff binding enthalpies for host  guest systems , and of even 
 greater precision for the ion model systems , offers more than a consistency 
 check . given that robust agreement between the direct and vant 
 hoff methods has seemed to be elusive for not only computational investigations but also experimental studies , there is value in a clear demonstration that rigorous consistency 
 can , in fact , be achieved .
second , consistency between the direct 
 and vant hoff binding enthalpy calculations arguably provides 
 strong support for the validity of the binding free energy method 
 employed .
this suggests that , as an alternative to assessing error 
 by performing multiple free energy simulations which close a thermodynamic 
 cycle , as is frequently done for alchemical transformations , one could instead compare direct and vant 
 hoff binding enthalpies along a single path .
whether this approach 
 is advantageous will depend on the specifics of the system and the 
 free energy calculation .
the present study also sheds light 
 on the relative merits of the direct versus vant hoff binding 
 enthalpy methods .
our data clearly indicate that , for small host  guest 
 systems , the direct binding enthalpy requires less net simulation 
 time to reach a given level of statistical uncertainty .
in addition , 
 the direct approach is far simpler to implement than the vant 
 hoff approach , at least when the system does not have multiple noninterconverting 
 configurations
. the vant hoff approach may still be numerically 
 favorable for larger systems , where the potential energy fluctuations 
 are so large that reducing the uncertainty of the mean to a useful 
 level requires extremely long simulations . on the other hand , the 
 challenge of achieving binding free energies that are sufficiently 
 tightly converged to yield meaningful temperature derivatives is not 
 trivial ; and the mere technical simplicity of the direct enthalpy 
 approach , which does not require a pathway or lambda windows , is itself 
 a powerful argument in its favor . with over 438 s of total sampling 
 across 4434 independent 
 simulations ,
this included generating the initial coordinates for each 
 simulation along the apr path , solvating the systems with precisely 
 the desired number of water molecules ( a feature that is not currently 
 available in amber s leap program ) , determining the correct 
 target values for restraints , and writing all of the input files necessary 
 to perform the simulations . additionally , our heterogeneous set of 
 computational resources ( desktop machines , local clusters , and national 
 clusters ) motivated us to develop execution / queuing scripts that simultaneously 
 optimized usage of each resource while maintaining compatibility with 
 the generalized input scripts .
of particular importance was implementation 
 of an automated scheme to determine whether an individual free energy 
 simulation window had reached a predetermined uncertainty threshold 
 in the force uncertainty .
our initial simulation set ( cb7-b2temp ) did not use this feature and , because all windows used the same 
 amount of sampling time , likely oversampled many of the easy - to - converge 
 windows near the pulled - free state .
implementation of our adaptive 
 scheme , which involved iteratively running a short 5 ns simulation 
 followed by checking the uncertainty value of the mean force , improved 
 the sampling efficiency of our other simulation sets .
we also note 
 here that the free energy simulation times listed in tables s1s3 all include the entire direct binding 
 enthalpy simulations in which the two end point simulations were extended 
 to at least 1 s each .
if one were interested in only computing 
 the binding free energy , then this extra sampling at the end points 
 would not be required , so the aggregate simulation time for a free 
 energy calculation would be lower .
nonetheless , in the cases presented 
 here , removing the additional sampling from the binding free energy 
 simulations would not tilt the balance in favor of the vant 
 hoff binding enthalpy approach .
for example , table 1 shows that the extended end point sampling 
 in the cb7-b2temp simulation set contributed only , at most , 
 2.7% of the total sampling time for the vant hoff calculations .
this study was feasible in large part to the availability of gpus , 
 both the commodity and research - grade varieties .
typical simulation 
 speeds for the systems studied here ranged between 150 and 200 ns / day 
 without hmr and up to 350 ns / day with hmr , depending on which generation 
 of hardware was used . with a small cluster of gpus
, both binding enthalpy 
 and free energy values can be computed within a matter of days to 
 highly precise values , with 95% confidence intervals in the range 
 of 0.5 kcal / mol .
as noted in the methods section , each host 
 studied here posed a distinct challenge to the calculation of well - converged 
 binding free energies . for rigid cb7 ,
the constrictive portal can 
 generate a cork popping problem , whereas flexible cd can catch 
 on the exiting guest and be pulled into long - lasting distorted conformations . 

we solved these problems by imposing restraints during the attach 
 phase and lifting them during the release phase , but a variety of 
 enhanced sampling technologies might also have been brought to bear .
indeed , given broad current 
 interest in binding free energy calculations , the computational tractability 
 of these host  guest systems , and the fact that each host examined 
 here exemplifies a different type of sampling challenges , these could 
 be valuable test cases for testing and improving enhanced sampling 
 algorithms .
broader 
 use of the methods 
 laid out here should be facilitated by the fact that our implementation 
 of the apr binding free energy method , and the associated direct calculation 
 of binding enthalpy , is naturally suited to run with nearly all simulation 
 packages that support restraints .
automation of the simulation setup still remains a challenge , as it 
 is difficult to foresee a generalized approach for determining the 
 ideal pulling pathway and implementing it with an appropriate restraint 
 setup .
however , in most cases , automation for a given host or protein 
 will require setup only for an initial guest or ligand , which can 
 then be used as a template for subsequent compounds .
as noted 
 above , we expect a primary application of computational calorimetry 
 to be force field validation and development . in anticipation of this
we anticipate that this approach will be informative since most force 
 fields are not optimized against experimental binding data , even though 
 the calculation of binding affinities , in the context of computer - aided 
 drug design , is one of their most important and prevalent uses .
moreover , 
 the fact that a relatively advanced force field can yield serious 
 errors in binding thermodynamics , as clearly demonstrated here , means 
 that the experimental data traditionally used to adjust force field 
 parameters , such as the hydration free energies of small molecules 
 and the properties of neat liquids , do not suffice when one s 
 goal is to compute the thermodynamics of noncovalent association . 

we anticipate that the use of binding data as a force field optimization 
 target , in combination with more traditional targets such as ab initio quantum data , neat liquid properties , and hydration 
 free energies , will be a powerful new tool to improved simulation 
 accuracy .
 we 
 present a strategy for carrying out high - precision calculations 
 of binding free energy and binding enthalpy values from molecular 
 dynamics simulations with explicit solvent . 
 the approach is used to 
 calculate the thermodynamic profiles for binding of nine small molecule 
 guests to either the cucurbit[7]uril ( cb7 ) or -cyclodextrin 
 ( cd ) host . for these systems 
 , calculations using commodity 
 hardware can yield binding free energy and binding enthalpy values 
 with a precision of 0.5 kcal / mol ( 95% ci ) in a matter of days . 
 
 crucially , the self - consistency of the approach is established by 
 calculating the binding enthalpy directly , via end point potential 
 energy calculations , and indirectly , via the temperature dependence 
 of the binding free energy , i.e. , by the vant hoff equation . 
 
 excellent agreement between the direct and vant hoff methods 
 is demonstrated for both host  guest systems and an ion - pair 
 model system for which particularly well - converged results are attainable . 
 additionally , we find that hydrogen mass repartitioning allows marked 
 acceleration of the calculations with no discernible cost in precision 
 or accuracy . finally , we provide guidance for accurately assessing 
 numerical uncertainty of the results in settings where complex correlations 
 in the time series can pose challenges to statistical analysis 
 . the 
 routine nature and high precision of these binding calculations opens 
 the possibility of including measured binding thermodynamics as target 
 data in force field optimization so that simulations may be used to 
 reliably interpret experimental data and guide molecular design . 
Introduction
Methods
Results
Discussion