equilibrium simulation 
 of molecular systems entails sampling of 
 conformations from the appropriate distribution corresponding to a 
 thermodynamic ensemble , such as the boltzmann distribution for the 
 isothermal canonical ensemble .
given a potential energy function that 
 describes the interatomic interactions , equilibrium sampling is required 
 to compute the thermodynamic observables , such as internal energy 
 and heat capacity .
development of efficient sampling methods is an 
 important and highly active field of research in biomolecular simulation 
 and molecular science in general .
most thermodynamic sampling methods 
 are based on either the molecular 
 dynamics ( md ) or monte carlo ( mc ) approaches .
md simulations are prone 
 to trapping in local minima of the potential energy surface ( pes ) 
 if there exist barriers that are larger than the available thermal 
 energy .
consequently , long simulations may be required for equilibration 
 on the high dimensional energy landscapes corresponding to biomolecules 
 of practical interest .
mc simulation involves random perturbations , 
 or moves , in the conformational space , designed to preserve the canonical 
 distribution .
due to the random perturbations , the system may directly 
 jump between local minima of the pes .
however , due to the bonded topology 
 and compact structures of biomolecules , large jumps in the conformational 
 space typically require cooperative motion of multiple atoms .
designing 
 such cooperative moves , which satisfy detailed balance and are computationally 
 efficient , is the primary difficulty in the mc simulation of biomolecules .
we recently developed a method , namely , 
 kirkwood sampling , for surveying the n = 3 m  6 dimensional conformational space of a molecule 
 containing m atoms .
the key challenge in random conformational 
 sampling of compact biomolecules is to avoid steric clashes .
kirkwood 
 sampling addresses this problem by incorporating correlations among 
 internal coordinates , as captured by the joint probability distribution 
 between them .
the joint distributions may be obtained from relatively 
 short high temperature md or mc simulation trajectories , or even using 
 an informatics approach by exploiting statistics from structural databases 
 such as the protein data bank or the cambridge 
 structural database .
results for small 
 molecules and small peptides suggest that incorporating low order ( pairwise 
 and/or 3-fold ) correlations may be sufficient to greatly reduce the 
 occurrence of steric clashes in comparison with sampling that ignores 
 all correlations .
this is fortuitous , since , in practice , sufficient 
 data is likely to be available for populating only the low order pdf s . 

furthermore , neglect of the higher order correlations leads to greater 
 coverage of the conformational space compared to the data used to 
 populate the joint distribution functions .
kirkwood sampling 
 is a geometrical conformational sampling method , 
 independent of the potential energy surface . in the present contribution , 

we describe two new algorithms , which combine kirkwood sampling with 
 existing monte carlo and replica exchange based approaches .
the idea 
 in each case is to improve convergence by taking advantage of the 
 greater conformational space coverage and random ( or non - markovian ) 
 sampling properties of the kirkwood procedure .
section 2 describes the various 
 aspects of the underlying theory , including the kirkwood sampling 
 algorithm , application of kirkwood sampling for conformational sampling 
 of molecules , and the description of the new equilibrium sampling 
 algorithms introduced in this work .
section 3 defines a model nine - atom molecule with an unbranched chain topology , 
 characterizes its energy landscape , and presents a detailed analysis 
 of the convergence of the new algorithms .
section 3 concludes with results for alanine dipeptide , a popular molecule for benchmarking simulation approaches .
let x1 , ... , xn be n discrete random variables 
 such that 
 the ith variable takes di discrete values , xi  { vi,1 , ... , vi , di}. we denote the kth order probability distribution function ( pdf ) of a set 
 of k such variables ( corresponding to coordinates 
 in the present work ) { x1 , ... , xk } as pk(x1 , ... , xk ) . using lower case letters 
 for specific values for a random variable , pk(x1 ,
... , xk ) , denotes pk(x1 = x1 , ... , xk = xk ) .
kirkwood sampling refers to a 
 family of algorithms for sampling points in the full n - dimensional space consistent with select joint distributions among 
 different subsets of the variables , as described previously .
the doublet level kirkwood sampling employs only 
 the 1-d , or singlet , p1(xi ) , and 2-d , or doublet , p2(xi , xj ) , pdf s . in this work , 

we employ the doublet level sampling algorithm , but extensions using 
 different sets of pdf s are feasible.algorithm 1 ( main text ) presents the pseudo 
 code for generating a point , x = ( x1 , ... , xn ) , in the n - dimensional space using the doublet 
 level sampling algorithm .
the variables are sampled sequentially 
 from their corresponding 
 one - dimensional conditional pdf , p1(2)(xk|x1 , .. , xk1 ) .
the normalization factor , nk , of the conditional pdf in algorithm 1 is obtained numerically by summing the 
 probability for all possible values of xk , the variable to be sampled .
the various kirkwood 
 sampling algorithms differ in the expression for the conditional probability 
 distribution .
the doublet level algorithm samples points from the n - dimensional probability distribution1following previous notation , the superscript  ( 2 )  in the above equations 
 denotes doublet level and will henceforth be dropped .
the sampling 
 probability , pn(2 ) , can be readily computed 
 for any given point .
the computational complexity for doublet 
 level kirkwood sampling ( algorithm 1 ) is o(n ) , since there are n(n  1)/2 doublet pdf s . 

the kirkwood sampling distributions are normalized by construction 
 via normalization of the conditional probability distribution of each 
 variable . in general
, the kirkwood sampling distribution involves 
 the product 
 of all pdf s employed . as a result ,
if certain cells in the 
 pdf s have zero probability , then the corresponding combinations 
 of coordinates will never be generated .
thus , if the input pdf s 
 contain such zero probability cells , or  holes 
, then , 
 strictly speaking , kirkwood sampling will not be ergodic , since it 
 eliminates certain regions of the conformational space . in principle
, 
 this problem can be remedied by filling the holes with a small but 
 finite probability . however , if the pdf s employed are representative 
 of the true distributions , then the conformational regions eliminated 
 by the holes are likely to correspond to high potential energy due 
 to unfavorable interactions such as steric clashes .
for example , the 
 unpopulated regions of the ramachandran plots for the backbone  
 and  torsion angles of a protein correspond to conformations 
 with steric clashes .
hence , one can think of the input pdf s 
 as constraints on the conformational space , and the accessible conformational 
 space shrinks as more pdf s are included .
one approach 
 for obtaining the input pdf s is to populate 
 them using conformations generated through an alternative sampling 
 approach . in the present work
, the input pdf s among the coordinates 
 were populated using a representative set of conformations of the 
 molecule obtained by molecular dynamics simulations . in practice , 
 since the order of the input pdf s is likely to be much smaller 
 than the dimensionality of the system , the conformational space accessible 
 to kirkwood sampling will invariably be much larger than that represented 
 by the conformations used to populate the pdf s .
this effect is illustrated for a three - dimensional 
 system in the supporting information .
note 
 that the representative set of conformations need not be exhaustive , 
 just large enough to ensure the absence of spurious holes . in practice , 
 even for high dimensional systems , sufficient data is likely to be 
 available for avoiding spurious holes in the low order ( 1- and 2-d ) 
 pdf s . in the present 
 context of molecular conformational sampling
, the random variables , xi , denote discretized internal 
 coordinates of a molecule .
following previous work , we employ the 
 bond - angle - torsion ( bat ) internal coordinate system .
the singlet and doublet pdf s are obtained as histograms 
 of the internal coordinate values observed using md .
individual binning 
 of each coordinate provides the n singlet ( or 1-d ) 
 pdf s , and joint binning of all pairs of coordinates provides 
 the n(n 
each coordinate is discretized into b equally spaced bins between 
 the minimum and maximum values observed in md .
the random variables xi now denote a bin number and 
 can take values in the discrete set { 1 , ... , b } , that is , for all 
 variables di = b and vi , j = j. the bin width for discretizing the ith coordinate is set to i = i / b , where i is the range of values observed in the simulation . to visualize 
 a conformation and
for computing energies , the bin numbers for the 
 discrete bat coordinates need to be mapped into real values , and subsequently 
 to cartesian coordinates .
a discrete bat coordinate , xi  { 1 , ... , b } , is mapped to 
 the continuous space value , i  , by 
 picking a random point in the corresponding 
 bin2where i , min is the minimum value for the ith coordinate and r    is 
 a uniformly distributed random number 
 in [ 0 , 1 ] .
the left edge of the bin is given by the first two terms 
 of eq 2 .
the conformational 
 distribution for a molecule in contact with a heat bath is given by 
 the boltzmann distribution .
assuming translational and rotational 
 invariance of the potential energy , the boltzmann distribution in 
 terms of internal coordinates is given by3where 
  = kbt is the inverse temperature 
 and kb is the boltzmann constant .
in eq 3 , e( ) is the potential 
 energy and j( ) is the jacobian for 
 the transformation from internal to cartesian coordinates . to present 
 the equations in the following sections using a more familiar notation
, 
 we define4so that 
 the boltzmann distribution can be 
 written as5 monte carlo 
 ( mc ) sampling is a general method for generating points from a multidimensional 
 probability distribution function .
mc sampling involves generation 
 of trial points or moves , which are then subjected to an acceptance 
 criterion designed to impose detailed balance with respect to the 
 distribution of interest , here the boltzmann distribution .
the first 
 mc algorithm , referred to here as perturbation mc ( pmc ) , is a standard 
 implementation of mc for molecular systems , where a trial conformation , 
 t , is generated by random perturbation of 
 each coordinate of the current conformation , c,6with r being a uniform random 
 number in [ 0.5 , 0.5 ] and s a vector 
 of step sizes for the n coordinates . the trial conformation 
 is accepted or rejected using the metropolis acceptance function7the 
 trial conformation is accepted and added 
 to the markov chain if fpmc is greater 
 than or equal to a random number uniformly generated in [ 0 , 1 ] ; otherwise , 
 the move is rejected and the current conformation is added to the 
 chain . during equilibration , components of the step size vector s can be adjusted to achieve a specified acceptance 
 ratio .
after an optimal s has been found , 
 it must be fixed for production runs to guarantee detailed balance.algorithm 2 gives the 
 pseudo code for perturbation mc .
the second mc algorithm is 
 a variant of standard biased mc ( bmc ) 
 where the trial conformations are drawn 
 from a probability distribution , termed the biasing distribution , bn( ) , which is 
 defined over the full conformational space8 in contrast 
 to pmc , which is sequential , the 
 moves in biased mc are independent of the current conformation .
figure 1 schematically compares how the conformational space 
 is explored by the two schemes . for biased mc
, the acceptance function 
 of pmc is modified to reweight the current and the trial conformation 
 according to the boltzmann distribution9algorithm 3 provides 
 the pseudo code for biased mc .
given a temperature , the acceptance 
 ratio in a biased mc simulation is completely determined by the biasing 
 distribution .
mc moves are more likely to be accepted if the overlap 
 between the biasing and the boltzmann distribution is high .
indeed , 
 all moves will be accepted if the biasing distribution matches the 
 boltzmann distribution .
biased move mc can be trivially parallelized , 
 since the trial conformations are generated independent of the current 
 conformation in the markov chain . as a result , unlike perturbation 
 mc , no equilibration is required for a biased mc simulation . here
, 
 we apply the biased mc algorithm with the trial moves generated by 
 kirkwood sampling .
the key objective of this work is to obtain a boltzmann 
 distributed set of conformations corresponding to a given potential 
 energy function and temperature , using conformations generated by 
 kirkwood sampling .
schematics for the exploration of conformational 
 space by ( a ) perturbation 
 monte carlo ( algorithm 2 ) and ( b ) the biased monte carlo algorithm 
 ( algorithm 3 ) .
the square box represents the conformational space , 
 arrows represent perturbation moves , and filled circles represent 
 conformations . in the case of biased mc , conformations are sampled 
 independently from a biasing distribution .
replica exchange ( rex ) is an enhanced sampling method designed 
 to overcome the trapping problem in canonical simulations . in rex , multiple canonical simulations , or replicas ,
are run in 
 parallel and occasionally exchanges between one or more pairs of replicas 
 are attempted . in temperature replica exchange ( t - rex ) , all replicas 
 are run using the same potential energy function but at different 
 temperatures .
the lowest replica temperature is usually the temperature 
 of interest where the trapping problem is most severe .
exchanges with 
 the high temperature replicas help the low temperature replicas to 
 escape traps , thereby enhancing sampling in the low temperature replicas . 
 in the present work
, mc simulation is used for canonical sampling 
 for each replica so that the exchanges depend only on the potential 
 energy , and not the kinetic energy .
let 1 , 
 1 and 2 , 2 be the inverse temperatures and instantaneous conformations 
 of two replicas between which an exchange is being attempted .
the 
 exchange probability function in t - rex is given by10the 
 conformations are exchanged if ft - rex is greater than or equal to a uniform 
 random number in [ 0 , 1 ] . usually ,
exchanges are attempted between 
 replicas that are adjacent in the temperature ladder , though other 
 schemes have also been suggested.algorithm 
 4 provides the pseudo code for the present implementation of 
 t - rex . in t - rex , each replica may be considered as a reservoir 
 of boltzmann 
 distributed conformations for other replicas , albeit at different 
 temperatures .
exchanges may also be performed with a reservoir of 
 non - boltzmann distributed conformations , as long as the probability 
 distribution of the reservoir is known .
kirkwood sampling satisfies 
 this condition , as the probability of sampling a given conformation 
 can be computed ( here , using eq 2 ) .
a temperature 
 replica exchange where one ( or more ) temperature replica(s ) are coupled 
 to a kirkwood sampling scheme will be referred to as kirkwood reservoir 
 replica exchange ( kr - rex ) .
the probability for exchanges between a 
 temperature replica at inverse temperature  and a reservoir 
 with distribution bn( ) 
 is given by11where  is 
 the conformation from the temperature replica and r is a conformation drawn from the reservoir , that is , r
note that the acceptance function involves the reservoir probability 
 of the conformation from the temperature replica , bn( ) , 
 as well as the potential energy of the reservoir conformation , u(r; ) .
the likelihood of 
 acceptance from the reservoir depends on the conformational overlap 
 between the reservoir and the actual boltzmann distributions .
the acceptance ratio for biased mc simulation at a given temperature 
 gives the probability of exchange with a replica at the same temperature . 

therefore , given an acceptance ratio for exchanges with the reservoir , 
 short biased mc simulations are performed at multiple temperatures 
 to determine the lowest temperature to which the kirkwood reservoir 
 can be coupled .
schematics of the exchange protocols and exchange 
 probability function for t - rex and kr - rex are shown in figure 2 .
exchanges with the reservoir are attempted at 
 every alternate exchange cycle , as detailed in algorithm 
 5 .
schematic representation of ( a ) temperature 
 ( t - rex , algorithm 4 ) 
 and ( b ) kirkwood reservoir ( kr - rex , algorithm 5 ) replica exchange . 

the horizontal lines represent different replicas with the lower lines 
 corresponding to lower temperatures .
the green dotted box in part b represents the kirkwood 
 reservoir , filled green circles represent conformations sampled from 
 the reservoir distribution , and green arrows indicate exchange attempts 
 between the reservoir and the highest temperature replica . note that 
 conformations in the reservoir are not updated with conformations 
 from the temperature replicas .
the kirkwood replica exchange 
 strategy can efficiently sample thermodynamic states by coupling a 
 set of temperature replicas to a reservoir of kirkwood - generated structures . 
 although the trial configurations are sampled from the kirkwood distribution , 
 which is different from the underlying boltzmann distribution , the 
 acceptance criterion , eq 11 above , guarantees 
 that the move satisfies detailed balance , ensuring that each of the 
 temperature replicas samples its correct equilibrium distribution .
we employ the overlapping distribution method to verify that the kirkwood exchange moves do not disturb 
 the underlying equilibrium ensembles in the temperature replicas . 

this approach was originally introduced by bennett for free energy calculations but can also be used as a 
 consistency check on equilibrium sampling for replica exchange simulations .
we wish to 
 verify that each of our temperature replicas samples 
 its corresponding equilibrium distribution .
let us consider a pair 
 of temperature replicas , a and b , at inverse temperatures a and b , respectively . in the canonical ensemble , 
 a particular replica ,
a , for example , samples energy e with distribution12where (e ) is the energy 
 density of states and f = a1 ln z is the helmholtz free energy of the system . 

( we emphasize that the relevant quantity to analyze with this method 
 is the energy e , not the quantity u which contains the jacobian factors . )
we have a similar expression 
 for p corresponding to replica b. although p and p are generally 
 unknown ( due to the unknown density of states ) , they are related to 
 each other in a simple way as13here ,  = b  a and f = b  f  
 a  f is the reduced 
 free energy difference between the two replicas . hence , these two 
 energy distributions are connected via the unknown constant f . this relation can be seen as a manifestation of crooks 
 fluctuation theorem , and is in general 
 true only when both distributions sample their respective equilibrium 
 states , eq 12 above . hence , a verification of 
 eq 13 for the temperature replicas in a kirkwood 
 reservoir replica exchange simulation provides a consistency check 
 on the sampling . to proceed with this consistency check ,
we follow 
 bennett and define two functions:14a14balthough la and lb are functions 
 of e , from 
 eq 13 , we see that15is a constant , whose value corresponds to 
 the reduced free energy difference .
this relationship can be verified 
 from simulation data by binning the sampled energies from our various 
 replicas , constructing the functions la and lb , and plotting their difference 
 in each bin .
a plot of l as a function of e should , within statistical errors , provide a horizontal 
 line with a slope of zero . in section 3.1.4
, we apply this test to our replica exchange simulations in order 
 to verify equilibrium sampling .
the equilibrium sampling algorithms , described in section 2 , were implemented in matlab and octave .
we first applied 
 the methods on a model nine - atom chain molecule ( figure 3 ) with a simplified force field .
the model system was constructed 
 for computational efficiency and for ease of implementing internal 
 coordinate monte carlo moves .
kirkwood based biased mc and reservoir 
 replica exchange are applied to the model system , and results are 
 compared with conventional alternative approaches .
the simulations 
 were designed with the objective of validating the new sampling algorithms 
 and describe the considerations for setting up the simulations .
the 
 last section presents biased mc results for alanine dipeptide illustrating 
 the applicability to small biomolecules .
( a ) bond - angle - torsion 
 ( bat ) coordinate system for an m atom chain molecule 
 and ( b ) electrostatic charge and lennard - jones 
 well - depth of the test system studied in this work .
atoms without 
 charge and lennard - jones interaction are represented by unfilled dashed 
 circles .
this 
 section describes 
 the force field and coordinate system for the model system , characterizes 
 its energy landscape , and generates benchmark results using an independent 
 well - converged md replica exchange simulation
. we then present results 
 for perturbation and biased mc simulations applied to our system at 
 progressively lower temperatures .
we will see that both of these methods 
 fail to converge below a certain temperature .
finally , we present 
 results from t - rex and kr - rex replica exchange simulations , which 
 enhance convergence at the lower temperatures .
the number of mc steps 
 is 5  10 for all mc and replica exchange simulations 
 described in this section . the conformation of the model system is 
 specified by 21 bat coordinates , 
 which include eight bond lengths ( bi ) , seven bond angles ( ai ) , and six torsion angles ( ti ) .
figure 3 gives the atom labeling 
 and the definition of the bat coordinates . the functional form of 
 the potential energy of the molecule employs a molecular mechanics - type 
 force field , which includes bonded and nonbonded terms16where m = 9 is the number 
 of atoms .
the terms associated with the bonded , angular , and torsional 
 degrees of freedom are eb , ea , and et , respectively .
a 
 harmonic functional form was used for the bonded and angular contributions , 
 while the ryaekart 
bellman 
 functional form was used for the torsional term . in eq 16 , elj and
ee denote the nonbonded lennard - jones ( lj ) and electrostatic 
 interactions , respectively , and ri , j denotes the distance between atoms i and j. to compute distances between atoms , 
 the bat coordinates were transformed to anchored cartesian coordinates . in the present study ,
the potential energy of the test system included 
 all bonded terms but only selected nonbonded contributions ( see figure 3 and the supporting information ) .
lj interactions were considered only 
 for atom pairs ( 1,7 ) and ( 3,7 ) , while electrostatic interactions were 
 included only for atom pairs ( 1,7 ) , ( 3,7 ) , ( 1,9 ) , and ( 3,9 ) .
atoms 
 1 , 3 , and 9 carried positive charges of + 0.1 , + 0.1 , and + 0.2 , respectively . 

atom 7 carried a negative charge of 0.4 , and the remaining 
 atoms were neutral .
fewer nonbonded interactions were included to 
 reduce the cost of the energy evaluation and obtain accurately converged 
 results for the benchmarking more efficiently .
the full energy function 
 and 
 its parameters and the input files for gromacs and amber are provided as supporting information .
although the potential 
 that we have defined is relatively simple , 
 the corresponding landscape can still provide a useful benchmark to 
 compare different sampling schemes . in order to gain some idea of 
 this complexity , we sampled local minima using basin - hopping global 
 optimization , and then computed pathways connecting the global 
 minimum with all other minima using the doubly nudged elastic band method with accurate refinement of transition states 
 by hybrid eigenvector - following .
the 
 resulting database of minima and transition states was used to plot 
 the disconnectivity graph shown in figure 4 , which suggests a largely funnelled potential energy 
 surface .
the low energy minima adopt a more compact 
 conformation with similar geometry of the oppositely charged atoms . 

the energy of the global minimum was 22.6 kj / mol , and there 
 are three prominent funnels at energies below 16 kj / mol , with 
 barriers between funnels of roughly 4.14 kj / mol , corresponding to 
 a temperature of 500 k. thus , at temperatures significantly lower 
 than 500 k , the system is likely to be trapped in one of the low energy 
 funnels .
in this section , we investigate convergence of thermodynamic 
 quantities for temperatures ranging from 20 to 500 k. disconnectivity graph 
 constructed from a database of 236 minima 
 and 1454 transition states . branches leading to select minima ( a  e ) 
 are colored red .
conformations of these labeled minima are shown with 
 the positively charged atoms in blue and the negatively charged atom 
 in red ( see figure 3b for atom numbering and 
 the supporting information for potential 
 energy 
 definition ) .
distances between atom pairs ( 1,7 ) and ( 3,7 ) are also 
 shown .
note that the low energy minima ( b  e ) are more compact 
 than the high energy minimum ( a ) and have a similar arrangement of 
 atoms 1 , 3 , and 7 .
the energy of the global minimum ( d ) is 22.6 
 kj / mol , and that of minimum ( a ) is 9.9 kj / mol
. typical barriers 
 between minima belonging to the three low energy funnels are around 
 4.8 kj / mol , corresponding to a temperature of roughly 1000 k. reference canonical potential energy 
 distributions were generated 
 at different temperatures using md replica exchange .
twelve replica 
 temperatures were used : 10 , 20 , 30 , 50 , 75 , 100 , 130 , 165 , 200 , 300 , 
 400 , and 500 k. the md rex simulations were performed with gromacs 
 4.6.5 , using langevin dynamics with a 
 time step of 1 fs and a friction coefficient of 5 ps .
conformations and energies 
 were saved every 0.5 ps ( 500 time steps ) , giving 10 data 
 points per replica .
figure 5 shows the 
 potential energy distribution observed during the production run , 
 and table 1 displays the average energy and 
 heat capacity computed for the simulation temperatures using the potential 
 energy values observed in each replica . the heat capacity for the 
 replica at temperature t was computed using17the average energy and heat 
 capacity were 
 interpolated between each replica exchange simulation temperature 
 by applying the multihistogram method to the md rex potential energy distributions , using 
 an energy bin width of 0.1 kj / mol .
the peak heat capacity from the 
 multihistogram heat capacity curve is 26.6 kb/2 at 104.5 k. figure 6 shows the two 
 quantities computed at 50 equally spaced temperatures between 10 and 
 500 k. at the simulation temperatures , the interpolated values are 
 in good agreement with the values directly computed from the simulation 
 data ( see table 1 ) .
potential energy distributions 
 obtained from reference md replica 
 exchange simulation described in section 3.1.2 .
the distributions correspond to temperatures 10 , 20 , 30 , 50 , 75 , 
 100 , 130 , 165 , 200 , 300 , 400 , and 500 k , from left to right .
the inset 
 shows the energy distribution of doublet level kirkwood samples from 
 refpdfst500 ( line marked by boxes ) and refpdfst200 ( solid line by 
 circles ) overlaid on select canonical energy distributions ( unmarked 
 lines ) .
kirkwood samples generated using refpdfst200 and refpdfst500 
 have significant overlap with the canonical distribution at 200 and 
 500 k , respectively .
average energy and configurational heat capacity computed using 
 the md rex potential energy distributions shown in figure 5 .
bat coordinates were extracted 
 from all 12 trajectories , giving 
 12  10 data points for each coordinate , with 30 
 equally spaced bins between the minimum and maximum observed values 
 in each case .
the discretized coordinates were used to populate the 
 1-d ( singlet ) probability distribution for each coordinate and 2-d 
 ( doublet ) probability distribution for all pairs of coordinates . in 
 all , 21 singlet and 210 doublet distributions were populated and used 
 to sample conformations using the doublet level kirkwood algorithm 
 ( algorithm 1 ) .
note that , since data from all 
 replicas were pooled together , the conformations used to populate 
 the pdf s are not boltzmann distributed for a single temperature . 
 nevertheless , since more conformational space is covered at higher 
 temperatures , the configurations will be more representative of the 
 highest temperature replica at 500 k. we refer to this set of pdf s 
 as refpdfst500 .
another set of pdf s , referred to as refpdfst200 , 
 was generated using 2  10 conformations from a separate 
 200 ns md simulation at a lower temperature of 200 k. kirkwood 
 sampling was performed to generate 10 conformations 
 for each pdf set .
the figure 5 inset shows 
 the distribution of potential energies of the kirkwood samples overlaid 
 on the canonical potential energy distributions at different temperatures . 

the energy distribution of kirkwood samples generated from refpdfst500 
 ( curve marked by boxes ) has significant overlap with the canonical 
 distribution at 500 k , and the distribution corresponding to refpdfst200 
 ( curve marked by open circles ) has significant overlap with the canonical 
 distribution at 200 k. these results are consistent with earlier work , where the energy distribution of the kirkwood samples was found 
 to overlap strongly with the original boltzmann distribution from 
 which the conformations used to populate the reference pdf s 
 were generated .
indeed , kirkwood sampling was originally developed 
 for approximating the boltzmann distribution at a given temperature . 

note that the energy distributions of the kirkwood samples are shifted 
 to higher values relative to the original canonical distribution , 
 due to the greater coverage of the conformational space , and the fact 
 that there are more conformations at higher energies .
perturbation mc 
 ( algorithm 2 ) and biased mc 
 ( algorithm 3 ) simulations were performed for 
 successively lower temperatures starting from 500 k. both simulations 
 were performed at temperatures of 200 , 300 , 400 , and 500 k ; perturbation 
 mc was also conducted at the lower temperatures of 100 , 50 , and 20 
 k. the number of mc steps for all simulations was 5  10 . for the biased mc simulations
table 2 gives the average 
 energy and heat capacity computed from the potential energies from 
 different simulations and their difference with respect to the reference 
 values in table 1 .
the acceptance ratio observed 
 in the biased mc simulations and reference values from table 1 are also given .
considering the biased mc results first , figure 7 shows the distribution of potential energies for 
 each simulation 
 overlaid on the reference distributions from figure 5 .
the biased mc and reference distributions are in good 
 agreement for temperatures 300 k , consistent with the average 
 energy and heat capacity results in table 2 .
these results show that kirkwood sampling with a fixed set of input 
 pdf s can be used to obtain boltzmann distributed conformations 
 at different temperatures .
the biased mc algorithm may be viewed as a resampling approach wherein a set of kirkwood distributed conformations are resampled 
 to generate a boltzmann distributed set . in figure 7 , resampling effectively shifts the right - most distribution 
 ( line marked by boxes ) to match with the different boltzmann distributions . 

note that boltzmann distributions can be generated for temperatures 
 lower than that of the md simulation used to populate the input pdf s . 

biased mc simulations in figure 7 used input 
 pdf s effectively populated with 500 k md simulation data but 
 were able to generate boltzmann distributions at 300 k. potential energy 
 distributions ( in black ) from biased move mc simulations 
 ( section 3.1.3 ) at t = 200 , 
 300 , 400 , and 500 k using kirkwood moves with refpdfst500 .
the 
 potential energy distribution of the original kirkwood samples is 
 also shown marked by boxes .
acceptance ratios for the biased mc simulations were 0.29 , 
 0.25 , 
 0.07 , and 0.009 for t = 500 , 400 , 300 , and 200 k , 
 respectively .
the acceptance ratio falls as the temperature is reduced , 
 consistent with the decreasing overlap between the kirkwood potential 
 energy distribution and the canonical energy distributions ( see figure 5 ) .
note that the conformations obtained by kirkwood 
 sampling are completely determined by the input pdf s . as a 
 result , in contrast to perturbation mc , the acceptance ratio can not 
 be adjusted for a given temperature . in other words ,
a given acceptance 
 ratio would impose a lower limit on the temperatures for which biased 
 mc simulations can be run .
note that the overlap of the kirkwood distribution 
 with a target canonical distribution can be determined even in the 
 absence of the target distribution .
the acceptance ratio of a biased 
 mc simulation is a direct measure of the overlap .
moreover , since 
 no equilibration is required for biased mc simulations , a relatively 
 short run can provide an estimate of the acceptance ratio .
we 
 now characterize the perturbation mc simulations , which are 
 used in the next section for canonical sampling in the replica exchange 
 simulations . in each pmc simulation ,
a 10  000 step preliminary 
 run was performed to adjust the step size for an acceptance ratio 
 between 0.2 and 0.3 .
figure 8 shows the comparison 
 of the pmc potential energy distributions with the reference distributions 
 ( see figure 5 ) .
table 2 and figure 8 show that the pmc simulations 
 are in good agreement with the reference values for t  100 k , with higher temperatures producing better agreement . 

the deviations at low temperatures are due to the trapping of the 
 markov chain in local potential energy wells .
potential energy distributions 
 ( in black ) from the perturbation 
 mc simulations ( section 3.1.3 ) overlaid on 
 the reference distributions ( in blue ) .
the successful convergence for pmc at the low ( relative to 
 typical 
 barriers of the pes ) temperature of 100 k is likely due to the efficiency 
 of the internal coordinate bat moves . for comparison
, we performed 
 a 100 k pmc simulation using random cartesian moves , which was 10 
 times longer than the above bat move pmc simulation .
figure 9 shows that the energy distribution 
 from the cartesian move simulation is in much poorer agreement with 
 the reference distribution .
we note that , even though the coordinates 
 are perturbed independently , bat coordinate moves are particularly 
 effective here because of the unbranched chain topology of the molecule . 

absence of side chains greatly reduces the chances of steric clashes , 
 even if a torsion angle in the middle of the chain is perturbed substantially .
potential 
 energy distributions from perturbation mc simulations 
 at 100 k using bond - angle - torsion ( bat ) and cartesian coordinate moves . 

the cartesian move simulation was run for 5  10 steps , 
 while the bat move simulation was 10-fold shorter at 5  10 steps .
the distributions show that bat coordinate moves are 
 move efficient than cartesian coordinate moves .
we now compare temperature replica exchange ( algorithm 4 ) and reservoir replica exchange ( algorithm 5 ) , focusing on the low temperatures ( 100 
 k ) for which the mc simulations of the previous section failed to 
 converge .
replica exchange simulations were performed for two sets 
 of temperatures , namely , a high temperature set , 20 , 50 , 100 , and 
 200 k , and a low temperature set , 20 , 30 , 50 , and 100 k. the high 
 temperature set was chosen to inspect the convergence of the heat 
 capacity peak at 104.5 k , and the low temperature set was used to 
 inspect the impact of the reservoir on the convergence . in the case 
 of reservoir replica exchange simulations , for both sets
the kirkwood 
 reservoir for the high temperature set was based on refpdfst500 , and 
 that for the low temperature set was based on refpdfst200 .
perturbation 
 mc was used for canonical simulation in the temperature replicas , 
 and for each temperature , the step size was adjusted during equilibration 
 to achieve an acceptance ratio between 0.2 and 0.3 .
the length of 
 the production run was 5  10 steps , and exchanges 
 were attempted every 20 steps using the exchange protocols described 
 in algorithm 4 and algorithm 
 5 .
energy values were saved at every mc step generating 5  
 10 energy values per replica .
we first discuss the rex simulations 
 using the high temperature pdf set . for all rex simulations , we performed 
 the overlap test ( section 2.6 ) to validate 
 the implementation of the algorithms and , in the case of kr - rex simulations , 
 also verify that coupling to the reservoir did not disturb the boltzmann 
 distribution in the temperature replicas .
figure 10 shows the analysis for replicas 2 and 3 , corresponding to 
 50 and 100 k , respectively , for the kr - rex simulation .
the inset in 
 the lower panel of figure 10 also shows the 
 plot of l for the corresponding replicas in 
 the t - rex simulation .
we see that the function l is flat within statistical errors , indicating that the fluctuation 
 theorem , eq 13 , is satisfied , and hence both 
 replicas are sampling their respective canonical distributions correctly . 

the overlap test is 
 performed on replicas 2 and 3 , at temperatures 50 and 100 k , respectively , 
 taken from the kr - rex simulations described in section 3.1.4 using the high temperature pdf set .
lower panel : the overlap 
 functions l2 , l3 , and l ( eqs 14 and 15 ) . the function l is flat 
 within statistical error , indicating that the replicas are sampling 
 their correct equilibrium distributions ( see section 2.6 ) .
inset : the function l expanded 
 ( red crosses ) , compared to l as calculated 
 from the same temperatures in the t - rex simulation ( blue circles ) .
we next analyze and compare the 
 convergence of the t - rex and kr - rex 
 simulations using the high temperature set .
for both simulations , 
 the exchange acceptance ratios were 0.02 , 0.08 , and 0.08 for replica 
 pairs ( 20,50 ) , ( 50,100 ) , and ( 100,200 ) , respectively .
the lowest exchange 
 acceptance ratio of 0.02 corresponds to 5000 successful exchanges . 
 for kr - rex simulation ,
the number of successful exchanges with the 
 reservoir was roughly 1100 , consistent with the expected acceptance 
 ratio of 0.009 based on the biased mc simulation at 200 k. the convergence 
 of a replica exchange simulation is limited by the time needed for 
 a trajectory to diffuse between its lowest temperature , where it is 
 trapped in a metastable state , to a high temperature one , where it 
 can overcome its energetic barriers via enhanced thermal fluctuations 
 or the kirkwood reservoir . to analyze this convergence , we investigate 
 the diffusion of replica trajectories among the ladder of temperatures . 

figure 11 shows the trajectory that began in 
 the lowest temperature replica in the kr - rex simulation .
trajectory of the 20 k replica from the 
 kr - rex simulation described 
 in section 3.1.4 using the high temperature 
 pdf set .
the 
 trajectory for the first 5  10 exchange attempts 
 is depicted here .
figure 12 compares 
 the energy distributions of the t - rex and kr - rex simulations with 
 the reference distribution from md rex .
the three distributions are 
 essentially identical for all temperatures , including 20 and 50 k , 
 which failed to converge in the case of single temperature perturbation 
 mc simulations ( see figure 8) .
the t - rex and 
 kr - rex simulations were also compared in terms of convergence of the 
 average energy and heat capacity as a function of temperature , which 
 is a more stringent test of convergence than comparison of energy 
 distributions .
the average energy and heat capacity were obtained 
 by applying the multihistogram method to the potential energy distributions 
 computed using successively larger fractions of the simulation data . 

figure 13 shows the convergence of the average 
 energy and heat capacity overlaid on the reference curves , and figure 14 shows the difference between each curve and the 
 reference . for both t - rex and kr - rex
, the average energy converges 
 to within 0.5 kj / mol over the full temperature range of 20 and 200 
 k using only 10% ( 500  000 steps ) of the production data .
the 
 convergence of the heat capacity is , as expected , slower than that 
 for the average energy .
nevertheless , both the t - rex and kr - rex curves 
 converge to within 0.25 kb of the reference 
 with 50% of the production data .
these results show that the t - rex 
 and kr - rex simulations were consistent with each other .
however , from 
 these simulations , it is difficult to gauge the impact of the reservoir 
 on the convergence .
potential energy distributions for t - rex ( blue ) and kr - rex 
 ( red ) 
 simulations using the high temperature set described in section 3.1.4 .
distributions are shown for t = 20 ( left - most ) , 50 , 100 , and 200 k ( right - most ) .
the three distributions are nearly overlapping 
 for all temperatures , with t = 50 k ( inset ) showing 
 the largest deviations .
convergence of average energy ( left panels ) and heat capacity ( right 
 panels ) for t - rex ( top ) and kr - rex ( bottom ) simulations using the 
 high temperature pdf set ( section 3.1.4 ) .
the 
 legend indicates the fraction of data used from the 5  10 step production run .
the reference curves ( from figure 6 ) are in black . differences between the average energy and heat capacities in figure 13 and the reference values from figure 6 .
the horizontal dotted lines in the left panels 
 indicate energy differences of 0.2 kj / mol , and those in the 
 right panel indicate heat capacity differences of kb/4 .
the effect of the reservoir 
 in the high temperature set rex simulations was not discernible  probably 
 because the highest temperature of 200 k was sufficiently high to 
 avoid trapping , so that additional jumps in the conformation space 
 through exchanges with the reservoir did not improve convergence . 
 to investigate the effect of the reservoir
however , for temperatures 
 below 200 k , refpdfst500 could not be used to set up the reservoir 
 because of a low acceptance ratio .
we therefore created the refpdfst200 
 set of pdf s ( see section 3.1.2 ) using 
 data from a 200 k md simulation . recall that the potential energy 
 distribution of the kirkwood samples from this set has high overlap 
 with the canonical distribution at temperatures around 100 k ( see 
 figure 5 , inset ) .
the acceptance ratio for 
 a preliminary 10 step biased mc simulation at 100 k using 
 refpdfst200 was found to be 0.05 .
pdf set refpdfst200 could , therefore , 
 be used to set up the reservoir for kr - rex simulation using the lower 
 temperature set where the highest replica temperature is 100 k. figure 15 shows the energy distributions obtained 
 from the t - rex and kr - rex simulations for the low temperature set . 

the distributions from the kr - rex simulation ( lines marked by circles 
 in the bottom panel of figure 15 ) are substantially 
 closer to the reference distribution than those from the t - rex ( top 
 panel of figure 15 ) simulation .
table 3 gives the average energy and heat capacity at the 
 simulation temperatures computed from the raw simulation data , without 
 multihistogram analysis .
for the lowest temperature replica at 20 
 k , the average energy from the t - rex simulation deviates from the 
 reference by 0.62 kj / mol , while the deviation of average energy from 
 kr - rex is much smaller at 0.09 kj / mol .
the number of successful exchanges with 
 the reservoir was over 23  000 , consistent with the biased mc 
 acceptance ratio of 0.05 .
the results show that coupling with the 
 reservoir enhances sampling in the low temperature replicas .
potential 
 energy distributions ( lines marked by circles ) from the 
 5  10 steps ( a ) t - rex and ( b ) kr - rex simulations 
 described in section 3.1.4 using the low temperature 
 pdf set with replica temperatures of 20 k ( black ) , 30 k ( red ) , 50 
 k ( blue ) , and 100 k ( green ) .
the distributions from the kr - rex simulation 
 are much closer to the reference than those from the t - rex simulation , 
 indicating that coupling to the reservoir enhanced convergence .
the pdf set 
 refpdfst200 was used 
 to set up the reservoir in kr - rex simulation .
the numbers in parentheses 
 are differences with respect to the reference values . in this section
, 
 we apply biased monte carlo ( algorithm 3 ) to 
 alanine dipeptide .
similar to the model system , the bat internal coordinate 
 system was used for kirkwood sampling as described elsewhere .
the conformation of the 22-atom molecule was 
 defined using 21 bond lengths , 20 bond angles , and 19 bond torsions . 

the range for angular coordinates was [ 0,2 ] radians , 
 while that for bond coordinates was [ 0.6,2 ]  .
md data for populating 
 the pdf s was generated by a 1000 ns replica exchange simulation 
 with replica temperatures of 300 , 400 , and 500 k. in each replica , 
 coordinates were saved at regular intervals to generate 10 conformations for each temperature .
the all - atom amber99sb force field was used , and the simulations were 
 performed with gromacs 4.6.5 .
kirkwood 
 sampling was performed in matlab and 
 octave , and force field energy was computed 
 using openmm 6.2 . in terms of computational 
 cost
, most ( > 90% ) of the time was spent in generating the kirkwood 
 samples .
md data were used to populate the 60 singlet pdf s 
 and 1770 doublet pdf s corresponding to all pairwise combinations 
 of the coordinates .
two separate sets of pdf s , namely , refpdfst400 
 and refpdfst500 , were constructed using data from the 400 and 500 
 k replicas of the md simulation , respectively .
each pdf set was used 
 to generate 5 million conformations using doublet level kirkwood sampling . 

figure 16 shows the distribution of the potential 
 energy of the kirkwood samples overlaid on the distribution from the 
 md replica exchange simulation .
kirkwood samples from the two sets 
 of pdf s have good overlap with the reference distributions 
 for the corresponding temperatures , suggesting overlap between the 
 md and kirkwood sampled conformational spaces .
the energy distribution 
 for kirkwood samples follows the reference distribution such that 
 the distribution for refpdfst400 ( red dashed line ) is shifted to lower 
 energies as compared to the distribution using refpdfst500 ( blue dashed 
 line ) , consistent with the observations for the model system .
reference 
 distributions from md replica exchange simulation are in solid black 
 lines with 300 k being the left - most , followed by 400 and 500 k. the 
 dashed lines show distributions for 5  10 kirkwood 
 samples with red corresponding to refpdfs400 and blue to refpdfs500 
 pdf sets .
the marked lines are distributions from different biased 
 mc simulations . blue squares and circles correspond to biased mc using 
 refpdfs500 at 300 and 400 k , respectively .
we performed biased mc simulations to obtain a canonical 
 distribution 
 for temperatures lower than that of the original md simulation used 
 to populate the reference pdf s .
doublet level kirkwood sampling 
 is employed to set up the biasing distribution , and 5  10 mc steps are taken for all simulations in this section .
we 
 first present results for 300 and 400 k biased mc simulations with 
 refpdfst500 used for kirkwood sampling .
figure 16 shows the energy distribution obtained from the mc simulations overlaid 
 on the reference md potential energy distributions .
the 400 k distribution 
 ( blue circles ) is in good agreement with the target distribution , 
 unlike the 300 k distribution ( blue squares ) which is shifted to higher 
 energies .
the acceptance ratio at 300 k was lower ( 0.007 ) , consistent 
 with the poorer agreement of the 300 k energy distribution .
it is instructive to compare the ramachandran plots of the alanine 
 residue for ensembles generated by the different sampling methods . 
 in figure 17
, each two - dimensional plot shows 
 normalized joint distribution of  and  angles computed 
 using 10 conformations .
the left - most plot represents 
 the 500 k replica of the 1000 ns md rex simulation which was used 
 to populate the refpdfst500 pdf s .
the plot for kirkwood generated 
 conformations is in good agreement with the 500 k md plot , though 
 the kirkwood plot has a slightly smaller coverage .
recall that the 
 kirkwood sampling was performed while accounting for the pairwise 
 correlations of the  and  angles with all remaining 
 59 bat coordinates . by relaxing some correlations , for example
, ones 
 with side chain coordinates , a greater conformational space is likely 
 to become accessible by kirkwood sampling .
note that the number of 
 conformations ( 10 ) for the kirkwood ramachandran plot 
 is several orders of magnitude smaller than the number of time steps 
 ( 3  10 ) in the md rex simulation .
ramachandran plots for the alanine residue of alanine dipeptide 
 with  torsion on the horizontal axis and  on the vertical 
 axis .
the other plots employ refpdfst500 for kirkwood sampling , which was 
 populated using data from the 500 k md replica .
figure 17 also shows the ramachandran 
 plot 
 for the biased mc simulations using refpdfst500 and the reference 
 plots from the md replica exchange simulation .
biased mc effectively 
 resamples from the kirkwood distribution to reproduce a target distribution , 
 here the boltzmann distributions at 300 and 400 k. all regions of 
 the ramachandran plot that are sampled by the md simulation are represented 
 in the biased mc simulation . at both temperatures , the basins in the 
  < 0 region of the plot , which are highly populated in the 
 md simulation ,
the less populated 
 basins in the  > 0 region are poorly sampled by the biased 
 mc simulations , and may explain the deviations in the energy distribution 
 ( figure 16 ) . nevertheless , in spite of the 
 low acceptance ratio , since all regions are sampled in the biased 
 mc simulation , it would be advantageous to combine kirkwood sampling 
 with md or mc based local sampling , as demonstrated for the model 
 system above .
finally , we performed a t = 300 
 k biased mc simulation 
 using the refpdfst400 pdf set which was populated using conformations 
 from the 400 k md replica .
an acceptance ratio of 0.024 was obtained 
 in a 5  10 step simulation .
the energy distribution 
 of the mc samples using refpdfst400 ( marked by red triangles in figure 16 ) is in much better agreement with the reference 
 distribution as compared to that using refpdfst500 ( marked by blue 
 boxes ) .
note that the acceptance ratio for t = 400 
 k biased mc simulation using refpdfst500 is comparable to that of t = 300 k simulation using refpdfst400 .
one could envision 
 a multistage simulation to enable simulating arbitrarily low temperature 
 where each stage would consist of biased mc simulation followed by 
 repopulation of the pdf s using biased mc sampled conformations . 

the successive stage would simulate lower temperature with the first 
 stage employing pdf s populated by high temperature md or other 
 enhanced sampling methods .
in this contribution , we have presented two 
 methods for barrierless 
 equilibrium sampling of molecular systems .
our approach builds upon 
 kirkwood sampling , which employs low - order correlations among internal 
 coordinates of a molecule for random sampling of the conformation 
 space .
both methods make use of the property of kirkwood sampling 
 whereby the normalized probability of generating a given conformation 
 can be computed .
we have presented proof - of - concept results for the 
 new methods using a model system with nine atoms where the intramolecular 
 force field can be adjusted to highlight particular sampling issues . 

we also showed results for alanine dipeptide , a commonly used model 
 system for benchmarking sampling algorithms .
the first of the 
 two algorithms is based on biased monte carlo 
 where the kirkwood samples are used as mc moves .
in contrast to standard 
 molecular dynamics or monte carlo simulations , no equilibration is 
 required in the biased mc simulation , since successive kirkwood samples 
 are generated independently . as a result ,
biased mc using kirkwood 
 moves can be trivially parallelized in a distributed computing environment , 
 with no communication required between the compute nodes .
furthermore , 
 since kirkwood sampling is a geometrical sampling method , independent 
 of the energy landscape , the same set of samples can be used to generate 
 boltzmann distributions for different temperatures and potential energy 
 functions .
the convergence of a biased mc simulation depends 
 primarily on 
 the overlap of the kirkwood sampling distribution and the target boltzmann 
 distribution . the conformational space covered by kirkwood sampling 

is determined by the input set of probability distribution functions . 
 in the present work
the results showed that kirkwood samples 
 could be used to generate a boltzmann distribution , not only for the 
 temperature of the original md simulation but also for lower temperatures , 
 although the convergence slows as the temperature is reduced .
one 
 can imagine an iterative scheme where the initial set of pdf s 
 is constructed in a manner that provides coverage of a wide conformational 
 space .
for instance , the pdf s could be populated using high 
 temperature md , or using a database of conformations from the pdb 
 or a fragment pdf library . given this initial set of pdf s 
 and a potential energy function , one could then perform successive 
 stages of biased mc simulations and repopulation of the pdf s 
 to reach arbitrarily low temperatures .
the second algorithm 
 introduced in this work is a modification 
 of temperature replica exchange where a temperature replica is coupled 
 with a kirkwood reservoir of conformations .
exchanges with the reservoir 
 help to enhance sampling for that replica . as for biased mc ,
the acceptance 
 ratio for exchanges with the reservoir is determined by the overlap 
 of the kirkwood distribution and the boltzmann distribution for the 
 coupled replica .
note that in the absence of a reservoir the highest 
 temperature in a replica exchange simulation needs to be high enough 
 to overcome the energy barriers and avoid trapping . by coupling to 
 a kirkwood reservoir
, the highest temperature can be set independent 
 of the energy barriers . in the presence of a kirkwood reservoir ,
the 
 criterion for setting the highest temperature becomes the requirement 
 of sufficient overlap to facilitate frequent exchanges with the reservoir . 

thus , coupling with a reservoir effectively places a limit on the 
 highest temperature required , which may be lower than the highest 
 temperature dictated by the barriers of the pes . indeed ,
if the reservoir 
 has good overlap with the boltzmann distribution corresponding to 
 the temperature of interest , then just a single replica would suffice . 
 in that case ,
the temperature replica essentially performs local sampling , 
 while the reservoir facilitates global sampling of the conformational 
 space .
finally , we also note that the sampling for the lowest temperature 
 replicas can be enhanced by coupling to a reservoir constructed from 
 a set of low energy minima from the potential energy landscape . in this scenario ,
kirkwood samples could be 
 used for seeding basin - hopping simulations in order to generate 
 the low energy minima of the system , or in the mc part of a basin - sampling 
 calculation . in the present work , 
 canonical sampling for the temperature replicas 
 of a kr - rex simulation
the kirkwood 
 reservoir can also be coupled with temperature replicas evolving via 
 heat bath coupled molecular dynamics .
coupling to an md replica would 
 require a relatively minor modification of the current algorithm to 
 assign velocities to the reservoir generated conformations .
the velocities 
 would be sampled from the maxwell  boltzmann distribution corresponding to the temperature of the coupled 
 temperature replica .
it should therefore be straightforward to enhance 
 the performance of existing md based replica exchange codes by coupling 
 them to a kirkwood sampler .
therefore , it can be boltzmann 
 inverted to construct a reference energy 
 function .
we can then define intermediate potentials , which progress 
 from the potential energy of the physical force field to the kirkwood 
 reference energy .
the intermediate potentials could be used to perform 
 hamiltonian replica exchange with the 
 two boundary replicas corresponding to the reference and physical 
 energy functions .
all replicas may now be simulated at the temperature 
 of interest using standard perturbation move monte carlo .
this approach 
 would facilitate coupling of the kirkwood sampling with existing monte 
 carlo codes .
the number and spacing of the 
 intermediate replicas will depend on the overlap between the reference 
 and the physical energy surfaces .
we note that molecular dynamics 
 can not be used for canonical simulation with hamiltonian replicas , 
 since the reference energy surface is piecewise continuous and its 
 gradients are not well - defined .
this limitation arises because the 
 kirkwood distribution is constructed over a discretized conformational 
 space .
note that in this scheme kirkwood sampling is only used for 
 defining a reference energy function and not for generating mc moves . 

we are currently working on hamiltonian replica exchange using a kirkwood 
 reference energy and will present the results in a separate publication . 
 finally , nonequilibrium simulations , which switch the hamiltonian 
 between the physical and reference energy functions , can also be used to enhance replica 
 overlap and reduce the number of intermediate replicas .
application 
 of the present methods to systems of practical interest  from 
 small drug - like molecules to large proteins  would require a 
 flexible and scalable implementation of kirkwood sampling .
some of 
 the strategies that could be pursued in this direction are as follows . 
 in the present work
, the low order pdf s were populated using 
 conformations for the full molecule .
the pdf s can also be 
 generated for smaller molecular fragments and assembled for an arbitrary 
 molecule .
we are currently working on developing such a fragment pdf 
 library for proteins using peptide fragments up to five residues long . 

the fragment pdf s will be populated using pdb data and/or exhaustive md simulations . for small fragments , 

it is computationally feasible to calculate triplet or even quadruplet 
 pdf s , which would help better account for local packing .
we 
 note that the kirkwood framework can be used to sample a subset of 
 coordinates while keeping the other coordinates fixed . in the context 
 of conformational sampling of polymers ,
local packing can be accounted 
 for by incorporating correlations among coordinates of adjacent subunits . 

therefore , kirkwood sampling is likely to be particularly useful for 
 sampling small molecules , short peptides , or for sampling side - chain 
 conformations for a fixed backbone .
construction of pdf libraries 
 for proteins is relatively straightforward , 
 since proteins are built from a fixed set of amino acid residues . 
 in the case of small molecules
, it may be necessary to enumerate the 
 different possible bonded topologies for a given number of atoms and 
 atom types .
certain simplifications , such as treating rings or other 
 groups as rigid bodies , may also be beneficial . 
 for certain families of molecules , such as macrocycles
, it may be 
 necessary to explore alternative internal coordinate systems that 
 better capture the constraints imposed by the cyclic bonded topologies . 

we note that auxiliary variables could be used in the kirkwood framework 
 to impose constraints on the samples .
for example , for the system 
 studied in this work , by using end - to - end distance as a sampled variable , 
 one could generate samples with a given end - to - end distance .
in principle , kirkwood sampling can be performed using any number 
 of pdf s of different orders . in practice , however , it is desirable 
 to use as few pdf s as possible , since the computational cost 
 is proportional to the number of pdf s used and also the accessible 
 conformational space reduces with increasing pdf s
. therefore , 
 it would be useful to develop an adaptive scheme that successively 
 adds pdf s for different combinations of variables to achieve 
 a certain overlap with a target distribution .
the overlap could be 
 measured in terms of the acceptance ratio in a biased mc simulation . 

finally , we note that kirkwood sampling is a general method for sampling 
 points in high dimensional spaces and methods presented here may be 
 applicable to other areas which deal with high dimensional probability 
 distributions such as machine learning .
 we 
 present two methods for barrierless equilibrium sampling of 
 molecular systems based on the recently proposed kirkwood method ( j. chem . 
 phys.2009 , 130 , 
 134102 ) . 
 kirkwood sampling employs low - order correlations among internal 
 coordinates of a molecule for random ( or non - markovian ) sampling of 
 the high dimensional conformational space . 
 this is a geometrical sampling 
 method independent of the potential energy surface . 
 the first method 
 is a variant of biased monte carlo , where kirkwood sampling is used 
 for generating trial monte carlo moves . using this method , equilibrium 
 distributions corresponding to different temperatures and potential 
 energy functions 
 can be generated from a given set of low - order correlations . 
 since 
 kirkwood samples are generated independently , this method is 
 ideally suited for massively parallel distributed computing . 
 the second 
 approach is a variant of reservoir replica exchange , where kirkwood 
 sampling is used to construct a reservoir of conformations , which 
 exchanges conformations with the replicas performing equilibrium sampling 
 corresponding to different thermodynamic states . 
 coupling with the 
 kirkwood reservoir enhances sampling by facilitating global jumps 
 in the conformational space . 
 the efficiency of both methods depends 
 on the overlap of the kirkwood distribution with the target equilibrium 
 distribution . 
 we present proof - of - concept results for a model nine - atom 
 linear molecule and alanine dipeptide . 
Introduction
Theory
Results
Summary
and Discussion