A Partition Function Estimator
Ying-Chih Chiang Frank Otto Jonathan W. Essex


AIP




 Kobilka Institute of Innovative Drug Discovery, School of Medicine, The Chinese University of Hong Kong, Shenzhen, 2001 Longxiang Boulevard, 518172, Shenzhen, China
 chiangyc@cuhk.edu.cn


 Advanced Research Computing Centre, University College London, WC1H 9BT, UK
 f.otto@ucl.ac.uk





 Department of Chemistry, University of Southampton, SO17 1BJ, UK



§ ABSTRACT 
We propose a simple estimator that allows to calculate 
the absolute value of a system's partition function
from a finite sampling of its canonical ensemble.
The estimator utilizes a volume correction term to compensate
the effect that the finite sampling cannot cover the whole
configuration space.
As a proof of concept, the estimator is applied to calculate 
the partition function for several model systems, and 
the results are compared with the numerically exact solutions.
Excellent agreement is found, demonstrating that 
a solution for an efficient calculation of partition functions
is possible. 


A Partition Function Estimator
    Jonathan W. Essex
    September 9, 2024
==============================





§ INTRODUCTION

How to calculate the partition function of a system has long been 
an important question in physics because many quantities, such as
the free energy or the entropy, can then be determined afterwards.
In the past two decades, the question has been answered with 
various methods, such as the Wang-Landau algorithm <cit.> and nested sampling <cit.>. 
Both methods compute the (cumulative) density of states of a system
and allow a direct integration over the Boltzmann factor to yield 
the partition function as well as other quantities of interest. 
In particular, the nested sampling algorithm has demonstrated
its power in studying thermodynamic properties in material science <cit.>. 
The tested systems include Lennard-Jones clusters <cit.>,
metallic systems <cit.>, alloys <cit.>, 
and small water clusters <cit.>.
In these examples, the partition function is rarely the sole goal 
of the study, but rather a basic property that one can calculate 
apart from the heat capacity or a phase-diagram. 

However, calculating the partition function itself is still meaningful, 
since the ratio between two partition functions indicates 
the free energy difference between two systems (or two states). 
For instance, comparing the partition function of a system with
a protein and a free ligand to the partition function of the protein-ligand complex
yields the binding free energy between the protein and the ligand. 
This information can guide the design of drug molecules during lead optimization.
Yet, directly calculating the partition function is computationally expensive. 
Since the desired property is not a single partition function but the ratio between two, 
calculating this ratio using theories tailored to save computational resources seems more reasonable.
For instance, using importance sampling, the ratio between
two partition functions can be obtained by a single sampling on one
of the states, a method known as Zwanzig's equation <cit.>
or free energy perturbation (FEP) <cit.>. 
Together with post-processing methods such as Bennett's acceptance ratio (BAR) <cit.> 
or multistate Bennett acceptance ratio (MBAR) <cit.> to combine data 
from sampling over different states, multistep free energy 
perturbation (mFEP) <cit.> has became a standard approach for 
binding free energy calculation.
Other commonly used free energy methods include Kirkwood's thermodynamic integration <cit.>,
the Jarzynski equality <cit.>, etc. 
A comprehensive review can be found in the literature <cit.>. 
These free energy methods enable e.g. the calculation of the binding
free energy between a protein and a ligand <cit.>,
between two proteins <cit.>, or for evaluating the permeability of ligands 
through a lipid bilayer <cit.>.

Today, free energy calculations have become an indispensable tool in computer-aided drug design <cit.>.

While the aforementioned achievements makes one question the need for
calculating a single partition function directly, one may still wonder 
whether it is possible to calculate a single partition function more effectively, 
perhaps reaching a level where the performance of the calculation is comparable 
to these advanced methods, at least in some simple systems. 
To answer this question, we propose a Partition Function Estimator 
(PFE) that can compute the value of a single partition function from a 
finite sampling.  As a proof of concept, we apply this theory 
to several model examples, including the one-dimensional harmonic 
oscillator potential, the one-dimensional double-well potential, 
the two-dimensional Müller-Brown potential, and up to 30
Lennard-Jones particles in a three-dimensional box.
These examples are chosen because their references can be easily obtained:
For model potentials, exact numerical solutions are available via
brute-force integration, while the partition functions of 
Lennard-Jones particles can be obtained using nested sampling
or standard methods such as mFEP-MBAR.



§ THEORY AND METHODS




 §.§ The Partition Function Estimator

Consider a system composed of N identical particles of mass m. 
The Hamiltonian reads,

    H(𝐩,𝐪) = 𝐩^2/2m + U(𝐪)

where 𝐩 and 𝐪 denote the momenta and the spatial
coordinates of the particles, respectively. U(𝐪) denotes
the potential of the system.
The canonical partition function Z of the system is defined as,

    Z = 1/N! h^3N∫ e^-β H(𝐩,𝐪) d𝐩 d𝐪 ,

where h is Planck's constant, and β denotes the 
inverse temperature T multiplied with the Boltzmann constant
k_B (β=1/k_BT).
Z is a function of β, but this dependency is
dropped for ease of notation.
The integration over phase space in Eq. <ref> is separable,
and the integration over the momenta can be performed analytically.
The partition function of the system hence can be rewritten as <cit.>,

    Z = 1/N!( 2π m/β h^2)^3N/2∫ e^-β U(𝐪) d𝐪 .

With the integration over momentum space carried out, we can focus on 
the integration over coordinate space, and thus define

    Q = ∫ e^-β U(𝐪) d𝐪 ,

to denote the spatial contribution to the partition function Z.
In the same vein, the expectation value for any function f that
depends solely on the coordinates 𝐪, can be expressed as

    ⟨ f ⟩ = ∫ f(𝐪) e^-β U(𝐪)/Qd𝐪 .

Our discussion below will focus on how Q can be calculated from 
a finite sampling of the system's canonical ensemble. 



When sampling a system at finite temperature, the probability for
finding a microstate i with energy E_i is proportional to
exp(-β E_i), and consequently the probability for finding
any microstate with energy E is proportional to
g(E) exp(-β E) where g(E) denotes the density of states.
Due to the exponential nature of the Boltzmann factor, the energy
distribution of the system eventually decreases with increasing
energy – see Fig. <ref> for an illustration –
such that above a certain energy level E^*, the sampling becomes
insufficient, i.e. it is impractical to obtain enough samples with
energy E > E^* to reproduce the actual probability distribution
in this energy regime.
In practice, this is unproblematic if one is interested in the
sample average of quantities that don't grow exponentially with
E, such as the system energy itself and most other quantities
of physical interest. For these, the sample average is insensitive
to the high-energy tail of the probability distribution, and it is
of no concern if it was captured insufficiently during sampling.

Yet, in the extreme case one could be interested in evaluating
the sample average of the inverse Boltzmann factor b(E) = exp(+β E).  Its
theoretical expectation value is given by

    ⟨ b(E) ⟩
    = ∫ e^+β U(𝐪)e^-β U(𝐪)/Qd𝐪
    = L^3N/Q ,

where L^3 is the volume of a cubic box that the particles are confined to.
If we could obtain an estimate for ⟨ b (E) ⟩ via sampling, then
Eq. <ref> could be readily used to calculate Q. However,
in a sample of size n with energies E_i distributed according to
the distribution p(E) ∝exp(-β E),
the sample average of b is given by

    b̅ = 1/n∑_i=1^n e^+β E_i ,


which is extremely sensitive to the high-energy samples.
In practice, the insufficient sampling of the high-energy tail leads to a very
large fluctuation of b̅, making it useless as an estimate for
⟨ b ⟩ and thus for determining Q.
Note that here as well as below, we clearly distinguish between the
expectation value (⟨ b ⟩) and the sample average (b̅);
while the former is an exact (theoretical) value, the latter is empirically
obtained from a finite sample and will fluctuate if the sampling is repeated,
perhaps even failing to converge in case of insufficient sampling.

Let us now consider another function f(E;E^*) that reads,

    f(E;E^*) = e^+β E θ(E^*-E)  ,

where θ(E^*-E) denotes the Heaviside step function, 
which is 1 for E ≤ E^* and 0 for E > E^*.
In short, f(E;E^*) is the inverse Boltzmann factor,
but truncated to zero for energies E larger than some
chosen parameter E^*.
With E set to the potential energy, i.e. E=U(𝐪),
the expectation value of f is then given by,

    ⟨ f(E;E^*) ⟩ 
    = ∫ e^+β E θ(E^*-E) e^-β U(𝐪)/Q d𝐪
    = V(E^*)/Q ,

with V(E^*) defined as,

    V(E^*) = ∫θ(E^* - U(𝐪))  d𝐪 ,

which is the volume of the coordinate space where the potential 
energy is smaller than E^*.
Consequently, Q can be obtained via,

    ln Q = ln V(E^*)  - ln⟨ f(E;E^*) ⟩ .

Eq. <ref> is the proposed estimator.
We stress that this equation itself is exact, in that no approximations have
been made, and holds for any value of E^*.
Obviously, to make use of this equation, we still need to determine values
for both ⟨ f(E;E^*) ⟩ and V(E^*), along with
finding a “good” value of E^*.



 §.§ Finding E^*


Let us first concentrate on the expectation value ⟨ f(E;E^*) ⟩.
We wish to use the sample average of f, i.e.

    f̅
      = 1/n∑_i=1^n e^+β E_iθ(E^* - E_i)  ,

as an estimate for ⟨ f(E;E^*) ⟩.  In contrast to Eq. <ref>,
this sample average is not sensitive to the high-energy tail, provided that
E^* is chosen such that energies below E^* are all sufficiently sampled, 
cf. Fig. <ref>.
Indeed, if E^* were chosen to be vary large, then Eq. <ref>
would suffer from strong fluctuation due to insufficient sampling.
On the other hand, if E^* were chosen to be very small, then the
number of samples that contribute effectively to Eq. <ref>
would be very small, which again increases its error.
This suggests that there is an optimal choice for E^* that
can minimize the relative error in f̅.

According to Eq. <ref>, ln Q is given by the difference between
ln V(E^*) and ln⟨ f(E;E^*) ⟩. We are now using
lnf̅ as an estimate for the second term. Since this is the term
that comes from the sample average (we will discuss ln V(E^*) in the next section),
it makes sense to choose an E^* that can minimize the standard deviation of
lnf̅.  This is given by,

    σ_M = √(1/n ⟨ f(E;E^*)^2 ⟩ - ⟨ f(E;E^*) ⟩^2 /⟨ f(E;E^*) ⟩^2) ,

which is just the relative standard error of f(E;E^*), with n being the number of samples.
























To minimize this error, we find E^* such that

    ∂σ_M^2/∂ E^* =
    ∂/∂ E^*[ 1/n( ⟨ f(E;E^*)^2 ⟩/⟨ f(E;E^*) ⟩^2  -1 ) ]
    = 0   .

Expressing the expectation values in the above equation as integrals over energy space
by utilizing the density of states g(E), we have

    ⟨ f(E;E^*) ⟩ = ∫ e^+β Eθ(E^*-E)   g(E) e^-β E/Q d E = 1/Q∫θ(E^*-E)   g(E)  d E  ,

and

    ⟨ f(E;E^*)^2 ⟩ = ∫ e^+2β Eθ(E^*-E)^2   g(E) e^-β E/Q d E = 1/Q∫ e^+β Eθ(E^*-E)   g(E)  d E  .

Note that in the above equation, we have used the property of the Heaviside
function that θ^2 = θ. Further, the derivative of the Heaviside
function gives the Dirac delta function,

    ∂/∂ E^*θ(E^*-E) = δ(E^*-E)  ,

which is only non-zero when E = E^*.
Using Eqs. <ref>-<ref>, Eq. <ref> results
in the condition

    0    =   1/⟨ f(E;E^*) ⟩^2∂/∂ E^*⟨ f(E;E^*)^2 ⟩ -  2 ⟨ f(E;E^*)^2 ⟩/⟨ f(E;E^*) ⟩^3∂/∂ E^*⟨ f(E;E^*) ⟩
       =   1/⟨ f(E;E^*) ⟩^2[
    g(E^*)/Q e^+β E^* - 
    2 ⟨ f(E;E^*)^2 ⟩/⟨ f(E;E^*) ⟩g(E^*)/Q]  ,

or after simplification,

    e^+ β E^* = 2 ⟨ f(E;E^*)^2 ⟩/⟨ f(E;E^*) ⟩ .


According to Eq. <ref>, E^* depends on the expectation values of
f(E;E^*) and f(E;E^*)^2. Both can be estimated from
their respective sample averages. However, as both also depend on
E^*, this equation must be solved iteratively: 
First, one collects n samples e.g. via Monte Carlo sampling,
and saves the trajectory and the associated energies E_i.
Then, one calculates f̅ via Eq. <ref>
and similarly f^2, using an initial guess for E^*,
e.g. the maximum energy encountered during sampling.
Afterwards, E^* is updated via

    E^* = 1/β( ln 2 + lnf^2 - lnf̅)  ,

and this process is repeated until E^* converges.
Notably, each iteration can employ the same trajectory; one
merely has to update E^* and recalculate the sample averages,
which is computationally cheap.
Finally, one can use the optimized E^* to estimate
⟨ f(E;E^*)⟩ and σ_M from
f̅ and f^2.






 §.§ Calculating V(E^*)


With E^* determined, we can turn to the calculation of the volume term 
V(E^*). According to Eq. <ref>, V(E^*) is the 
volume of the coordinate space with E < E^*. 
For simple low-dimensional models, such as a harmonic oscillator or 
Müller-Brown potential, this volume can be calculated directly 
via “binning”, i.e. dividing the space into small bins,
assigning each sample to its corresponding bin, and simply
counting the number of occupied bins.
However, as the dimensionality increases, no trajectory
could cover the entire admissible coordinate space via binning.
Calculating V(E^*) by this or other naïve integration
approaches is hence not possible.
Second, V(E^*) soon becomes very small in comparison with the total volume
of the coordinate space. Thus a simple Monte Carlo integration would also be
insufficient to address this issue.

Fortunately, the nested sampling algorithm <cit.> was designed 
to tackle such a problem.
Here, we employ a slightly modified version of nested sampling, making it more
suitable for the purpose of Eq. <ref>.  
Briefly, the implemented algorithm is as follows:

    
  0. Randomly generate N “walkers” with energies smaller than a starting energy ceiling E_0. Each walker is an independent copy of the system that will be propagated in the course of the algorithm. The probability distribution for the initial walkers must be uniform over the coordinate space. The starting energy ceiling in principle should be the highest potential energy possible, but for practical purposes it is normally set to a very large value, such as 10^12 k_BT. The initial volume V_0 is then given by the volume of the entire coordinate space, V_0 = L^3N.
    
  1. Select a fraction p (0 < p < 1). This determines the new energy ceiling for the current (i-th) iteration, E_i = p · E_i-1.  Use a simple Monte Carlo integration to calculate the volume V(E_i), i.e. the volume of coordinate space with energy smaller than E_i. Assuming that n_i walkers have energies below E_i, we have V(E_i) = (n_i/N) · V(E_i-1).

    
  2. Relax the walkers whose energies are larger than E_i, and propagate them freely below this energy, until they represent a sample with uniform probability over the coordinate space with energy E < E_i. Rejection sampling is employed for this step.
    
  3. Repeat Steps 1 and 2 until E_i < E^* is reached. As it is unlikely for E_i to exactly hit E^*, interpolation may be required to determine V(E^*) from V(E_i) and V(E_i-1).


Compared to the original nested sampling, our implementation is more akin
to an iterative Monte Carlo integration.  Our approach differs 
from nested sampling in two aspects. First and most importantly, 
nested sampling 
throws away a fixed number of walkers in each iteration,
leading to a fixed ratio of volume truncation. 
For instance, if one throws away the walker with the highest 
energy in each iteration, then at the i-th iteration, nested sampling 
yields a volume V(E_i) = (1/N)^i · V_0, where the energy E_i is
determined by the energy of the walker that is currently thrown away.
In contrast, our algorithm would yield a volume V(E_i) = ∏ (n_i/N) · V_0.   
Thus, any sampling fluctuation reflects on E_i in nested sampling,
whereas it reflects on V(E_i) in our implementation.

Second, nested sampling does not "relax" the walkers as we did in Step 2.
Rather, it duplicates a walker with energy smaller than E_i randomly 
and propagates the cloned walker until it becomes uncorrelated from the original one. 
This operation seems to be more efficient than ours, but as our 
calculation is truncated at E^* rather than proceeding to the lowest potential
energy, the number of iterations required is much smaller.  

Third, when applying nested sampling to integrate over the Boltzmann factor,
more walkers are needed in order to better resolve the density of state in energy.
Since PFE does not require any such integration, we can use much less walkers or even
more aggressive energy truncation to calculate V(E^*).








§ RESULTS AND DISCUSSION




 §.§ Harmonic Oscillator

Our first example is the one-dimensional harmonic oscillator potential,
U(x)=kx^2/2.  Here we set the force constant as k=300 and k_BT=0.59616. 
The potential energy curve and the corresponding distributions are depicted in Fig. <ref>(a).  
In this system, we employ a simple Metropolis Monte Carlo algorithm for sampling, utilizing a total
of 10^6 steps and a step size of 0.1. 
The trajectory and sampled potential energy are recorded every 10 steps. 
Hundreds of independent samplings are conducted and processed using PFE (Eq. <ref>) 
to derive ln Q. The corresponding average and standard deviation (fluctuation) are then depicted 
in Fig. <ref>(b). 
Here, different E^* are selected to demonstrate their impact on PFE. The term a% indicates 
the criterion for selecting E^*: it is chosen such that for the top a% of samples in energy,
the corresponding value of f(E;E^*) is zero.

The outcomes obtained from the optimal E^* (calculated using Eq. <ref>) are labeled as 
“opt" (red curve). As expected, utilizing PFE (Eq. <ref>) using the best choice of E^*
(Eq. <ref>) leads to a converged result after approximately 10^4 steps, while the sampling
fluctuation diminishing rapidly with the step count. 
Conversely, the ln Q calculated using the maximum energy sampled as E^* (black curve, 0%) 
exhibits substantial fluctuation and fails to converge.  However, upon applying the cutoff 
through the Heaviside function, PFE results converge towards the exact solution (dashed line), 
albeit with slightly inferior performance compared to using the optimal E^*.
This behavior aligns with the notion that high-energy microstates are often undersampled,
leading to convergence challenges.

The standard error of ln⟨ f(E;E^*) ⟩ can also be calculated as shown in 
Eq. <ref> and compared with the fluctuation determined from the standard deviation of 
ln Q across the 100 sampling repetitions.
These results are displayed in Fig. <ref>(c).
With the exception of the 0% cutoff, the predicted standard error closely matches the fluctuation,
and increasing the sample size results in diminished error. 
This intriguing finding indicates that the fluctuation primarily arises from
inaccuracies in computing ln⟨ f(E;E^*) ⟩, as the ln V(E^*) can be fairly
accurately computed from the trajectory histogram using 100 bins.






 §.§ Double Well Potential

Next we consider a double well potential defined as follows,

    U(x) = 16h/x_0^4 x^2 (x-x_0)^2  ,

where h is the barrier height and x_0 is the position of the second minimum.
In this example, k_BT is maintained at 0.59616, while h and x_0 are configured to 10k_BT 
and 3, respectively, to introduce some sampling challenges.  For this system, a total of 10^6 
steps are sampled using both Monte Carlo and replica-exchange (RE) <cit.> techniques.  
The latter is a typical enhanced sampling method that simulates multiple replicas with
various temperatures simultaneously. The coordinates between different replicas can be
swapped to enhance the sampling efficiency. 
In the current calculation, 10 replicas are employed in total. The temperatures are linearly
scaled, with the lowest and highest k_BT values set to 0.59616 and 1.9872, respectively. 

The potential curve and the corresponding distributions are depicted in Fig. <ref>(a).
Due to the substantial barrier, Monte Carlo sampling encounters difficulties in traversing the barrier,
as evident in the blue distribution. In contrast, RE sampling facilitates frequent barrier crossings, 
resulting in a symmetric spatial distribution even for the replica with the lowest temperature, 
as observed in the orange distribution. The energy histograms exhibit similarities between both
sampling methods, except in the high-energy region.
Shown in Fig. <ref>(b) are ln Q calculated using PFE. 
With the exception of the curve labeled “RE", all results are obtained from an extensive
Monte Carlo simulation. Once again, the optimal E^* (opt) provides the most accurate estimation
of ln Q, although the result is not yet converged after 10^6 sampling steps. 
Owing to the sampling difficulty posed by the high barrier, the ln Q plot displays a plateau 
around 10^5 steps (ln Q=-0.9), as depicted by the red curve. 
This plateau arises because the simulation struggles to surmount the barrier (Boltzmann factor 
e^-10= 4.54*10^-5) in shorter sampling durations. 
However, with longer sampling durations, simulations eventually overcome the barrier,
exploring the second well. Consequently, the fluctuation increases, and ln Q approaches 
the exact solution.


In contrast, ln Q calculated with RE sampling and the optimal E^* shows an early 
convergence well before 10^5 steps (orange curve).
This observation underscores how sampling difficulties can influence the convergence rate
of calculations. However, as enhanced sampling methods are extensively documented 
in literature <cit.> and out of scope of this study,  our discussion centers
on the behavior of PFE.
Interestingly, when the sampling remains confined to one well 
(reflected in the plateau-like ln Q values),the deviation from the exact solution is around 0.7.
This magnitude is about the same as ln(2), implying that the primary source of error stems
from halving the size of V(E^*) when the simulation is trapped in one well.
Introducing enhanced sampling techniques here primarily aids in the accurate computation
of V(E^*), rather than to ⟨ f(E;E^*) ⟩.

In Fig. <ref>(c), a comparison is made between the standard error of 
ln⟨ f(E;E^*) ⟩ and the fluctuation (standard deviation) of ln Q.  
Unlike the previous example of harmonic oscillator, the fluctuation surpasses the standard error, 
suggesting that calculating ln V(E^*) can introduce a significant error, even when directly assessed
from the trajectory.
Nonetheless, increasing the number sampling steps can reduce the fluctuation, as 
demonstrated by the orange curve in panel(b), where the fluctuation eventually dwindles to insignificance.





 §.§ Müller-Brown Potential

Another often used model potential is the Müller-Brown potential, defined as follows:

    U(x,y) = ∑^4_k=1 A_k exp[ a_k(x-x^0_k)^2+b_k(x-x^0_k)(y-y^0_k)+c_k(y-y^0_k)^2 ] + U_0  ,

with A=(-200,-100,-170,15), a=(-1,-1,-6.5,0.7), b=(0,0,11,0.6), c=(-10,-10,-6.5,0.7).
In the original Müller-Brown potential, U_0=0. Here we set U_0=147.70 to ensure  
the potential minimum is at zero. The potential energy surface is shown in Fig. <ref>(a).
The potential features three local minima, with the highest barrier approximately at 107.
To evaluate the performance of PFE, three different temperatures are employed:
k_BT=100, k_BT=10, and k_BT=2, which correspond to a well-sampled trajectory, 
a poorly sampled trajectory, and a trapped trajectory, respectively. 
The trajectory histograms (proportional to the population) are depicted in Fig. <ref>(b-d).
It is important to note that in panels (b)-(d), the number of bins employed in each direction is 100, 
but the bin width is automatically adjusted, resulting in varying count scales across the panels.



The Monte Carlo simulation is once more employed to compute ln Q (conducted over 10^7 steps with
data saved every 10 steps). The volume V(E^*) is determined based on the 2-dimensional 
trajectory histogram, where the area of the populated bins is summed to give V(E^*). 
Results are illustrated in Fig. <ref>.  Like before, 100 independent trajectories are
processed to derive the mean and standard deviation (fluctuation). Various E^* values are
employed to showcase the effect of incorporating the Heaviside function to truncate the inverse 
Boltzmann factor: PFE with the optimal E^* (opt) consistently yields the best results, while PFE 
with the highest sampled energy (0%) performs the poorest.
That said, PFE with the optimal E^* steadily converges to the exact solution, exhibiting
smaller fluctuation compared to choices of E^*.
This underscores the significance of selecting E^* to minimize error (as per Eq. <ref>).

Noteworthy observations arise when contrasting the standard error of ln⟨ f(E;E^*) ⟩
with the fluctuation of ln Q: At low temperature (2k_BT, panel f), the sampling is notably 
localized in space, enabling an accurate calculation of ln V(E^*).  
Consequently, calculating ln V(E^*) scarcely 
contributes to the error in ln Q, with data predominantly aligning along the plot's diagonal. 
As the temperature rises (10k_BT, panel e), the sampling ventures beyond barriers to explore other 
local minima. However, the sampling proves insufficient, as indicated by the histogram in Fig. <ref>(c).
This results in a larger fluctuation, since calculating V(E^*) based on the 
trajectory can be less precise. Nonetheless, once the convergence is achieved, the fluctuation 
diminishes significantly. 
Finally, at sufficiently elevated temperatures (k_BT=100, panel d), the fluctuation and standard 
error once again align, although the alignment is no longer diagonal.
This finding demonstrates that calculating V(E^*) entails inherent errors, although they are
diminishable with an increasing number of steps, as exemplified by the red curve in panel (a).






 §.§ Lennard-Jones Particles

The last example to be discussed in a more realistic system: Lennard-Jones particles in a 3-dimensional
box with periodic boundary conditions.  In this example, the box length is fixed at 25 Å and
Lennard-Jones particles are sequentially inserted into the box. 
The parameters are taken from the OpenMM example (Ar atom) <cit.>:
ϵ = 0.238 kcal/mol, σ = 3.4 Å, with the potential shifted to zero at 3σ 
for cutoff. The temperature is set to 120 K. The dispersion correction is disabled for simplicity.
Once again we calculate ln Q using PFE with the optimal E^*. A basic Monte Carlo simulation 
is employed for data collection. The simulation consists of a 50000-step equilibration 
with a step size 1.0 and a 10^6-step production run, with the data being saved every 1000 steps.
The volume V(E^*) is determined using the modified nested sampling algorithm (an iterative
Monte Carlo integration method). This approach involves 200 walkers, 2000 equilibration steps, 
and an energy fraction 0.99 at each iteration. The relaxation steps and the number of walkers 
to be relaxed are dynamically adjusted during execution.   
Moreover, the modified nested sampling algorithm is utilized for the direct integration of ln Q,
using identical parameters as those for V(E^*) calculations. 
Additionally, the standard mFEP-MBAR method <cit.> is applied for ln Q calculation
as a reference.
Briefly, mFEP-MBAR computes the free energy difference between two states, e.g. ln (Q_N/Q_N-1)
with N stands for the number of particles. Since Q_1 is simply the box volume, subsequent values
ln Q_2, ln Q_3, ..., ln Q_N are derived by accumulating the mFEP-MBAR outcomes.
The mFEP-MBAR calculation is conducted using molecular dynamics and alchemical functions implemented 
in OpenMMTools <cit.>.  
Simulation parameters are meticulously selected to ensure a comparable level of precision to PFE.
This leads to a maximum number of steps 21210000 (21 windows, each comprising 10000 equilibration steps 
and 1000000 steps production run). 

Results of ln Q are depicted in Fig. <ref>(a). This time, 10 independent calculations are
carried out for each method to determine the mean and the corresponding standard deviation (fluctuation). 
All three methods yield identical ln Q, which is a straight line that increases with the number of
particles. This outcome is anticipated since the slope provides the chemical potential of the system.  
It is noteworthy that the standard deviation is significantly smaller compared to ln Q, implying that
a single calculation is adequate to obtain a precise ln Q. 
The multiple independent calculations are primarily conducted to assess the magnitude of ln Q 
fluctuations, which are slightly larger than the anticipated standard error, as depicted in 
Fig. <ref>(b).
Considering performance metrics, a comparison is made based on the number of steps required to 
complete the calculations.  Given the use of distinct sampling techniques, viz. Monte Carlo for
PFE and molecular dynamics for mFEP-MBAR, the average number of sampling steps over 10 computations 
is reported.
Fig. <ref>(c) presents the step count necessary for ln Q calculation.
While a direct numerical integration over the Boltzmann factor (labeled NS, blue) is expected
to be computationally expensive, it seems odd that mFEP-MBAR demands even more steps.  
This discrepancy arises because mFEP-MBAR is designed to efficiently compute the free energy difference
between two states, such as ln (Q_N/Q_N-1). To get ln Q, 
the free energy differences are accumulated
over N, resulting in a higher total step count.
What should really surprise the reader is when one compares the step counts needed to compute
ln (Q_N/Q_N-1), i.e. the free energy difference between two states. 
In this instance, PFE's performance slightly surpasses that of mFEP-MBAR, as evident from
the red and black curves in Fig. <ref>(d).
This intriguing outcome highlights that calculating ln Q at a given temperature can really be achieved
with the same computational effort as finite sampling.



Last but not least, we discuss the stability of PFE as an estimator.
As shown in Fig. <ref>(a), ln Q is calculated using PFE with varying 
fractions of energy for evaluating V(E^*) in each iteration. All other parameters are kept 
as previously described. Despite the averages of ln Q being nearly identical across different 
parameters in 10 calculations, the fluctuation of ln Q actually increases, see panel (b). 
Additionally, Fig. <ref>(c) indicates that adjusting the parameter can marginally
enhance the performance of PFE, although the reduction in effort is not much.






§ CONCLUSION

In this article we proposed a partition function estimator that is composed of an
inverse Boltzmann factor and a Heaviside step function, with a parameter E^*
that can be determined from minimizing the square of the sampling error.
This results in a working equation that evaluates the expectation value of the proposed
function, along with a volume term V(E^*) to account for the energy cutoff
imposed by the Heaviside function. We demonstrate that the volume term can be directly 
evaluated from the trajectory for 1- and 2- dimensional examples, while 
a modified nested sampling integration method is proposed for examples with higher dimensions.

As a proof of concept, the performance of the estimator is tested using model examples, 
including the harmonic oscillator potential, the double-well potential, the Müller-Brown potential,
and Lennard-Jones particles.  Good agreement with the numerically exact solution is found across
all the test cases. While this result is anticipated, in the case of Lennard-Jones particles, 
a performance comparable to, or even better than, the standard FEP method is noted.

Currently, PFE is not yet ready to handle large biological systems like the standard mFEP-MBAR
can, nor can it calculate the partition function at multiple temperatures in one single run
as nested sampling can. Yet, PFE remains intriguing and warrants further development as it
offers a fresh perspective on addressing a long-standing challenging task.



This work is supported by the Kobilka Institute of Innovative Drug Discovery (KIIDD), School of Medicine, Chinese University of Hong Kong (Shenzhen), and the Shenzhen Science, Technology, and Innovation Commission (SZSTI), grant number JCYJ20230807114206014. 
YCC acknowledges the Royal Society for their support of the precursor to this work (NF171278). Special thanks to Prof. Guanglian Li and Prof. Ye Mei for their valuable discussions, to Prof. Livia B. Pártay for a fruitful discussion on nested sampling, as well as to Dr. Lantian Yao, Prof. Tzong-Yi Lee, and staff at KIIDD for providing a supportive working environment that nurtured this work.




*


