# Kinetic reconstruction of free energies as a function of multiple order
parameters

Yagyik Goswami  Srikanth Sastry<EMAIL_ADDRESS>https://www.jncasr.ac.in/faculty/sastry/  Theoretical Sciences Unit and School
of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific
Research, Bengaluru, India.

###### Abstract

A vast array of phenomena, ranging from chemical reactions to phase
transformations, are analysed in terms of a free energy surface defined with
respect to a single or multiple order parameters. Enhanced sampling methods
are typically used, especially in the presence of large free energy barriers,
to estimate free energies using biasing protocols and sampling of transition
paths. Kinetic reconstructions of free energy barriers of intermediate height
have been performed, with respect to a single order parameter, employing the
steady state properties of unconstrained simulation trajectories when barrier
crossing is achievable with reasonable computational effort. Considering such
cases, we describe a method to estimate free energy surfaces with respect to
multiple order parameters from a steady state ensemble of trajectories. The
approach applies to cases where the transition rates between pairs of order
parameter values considered is not affected by the presence of an absorbing
boundary, whereas the macroscopic fluxes and sampling probabilities are. We
demonstrate the applicability of our prescription on different test cases of
random walkers executing Brownian motion in order parameter space with an
underlying (free) energy landscape and discuss strategies to improve numerical
estimates of the fluxes and sampling. We next use this approach to reconstruct
the free energy surface for supercooled liquid silicon with respect to the
degree of crystallinity and density, from unconstrained molecular dynamics
simulations, and obtain results quantitatively consistent with earlier results
from umbrella sampling.

##  I Introduction

Numerical free energy calculations have contributed immensely to our
understanding of phase transitions and activated processes. In particular,
determining the underlying landscape allows one to quantify the relative
stability of the various states a system can exist in and also to
probabilistically predict the time evolution of the system. A specific area of
interest to the present work in which numerical free energy calculations have
contributed immensely to our understanding is that of polyamorphism in
liquidsStanley (2013). Network-forming liquids such as waterPalmer  _et al._
(2014); Debenedetti, Sciortino, and Zerze (2020), silicaChen, Lascaris, and
Palmer (2017), siliconVasisht and Sastry (2013); Goswami and Sastry (2022) and
other model liquids liquidsSmallenburg, Filion, and Sciortino (2014); Ricci
and Debenedetti (2017) have been shown to exhibit multiple metastable liquid
states, with an associated liquid-liquid phase transition, that are present
alongside the globally stable crystalline state. The deeply supercooled
conditions at which these systems exhibit this liquid-liquid phase transition
pose a challenge not just in experimental investigations but also in
performing simulations. At these conditions, crystallisation can occur very
rapidly (as in the case of silicon, which we address in this work), as a
result of which the sampling of the metastable liquid state is very poor and
enhanced sampling techniques such as umbrella samplingTorrie and Valleau
(1977) need to be employed to reconstruct the free energy landscape.

A large number of the cases in which such free energy calculations are
performed belong to the class of “rare events", where a direct observation of
the event is often unfeasible, even numerically. In this context, many free
energy calculation methods have been devised to obtain accurate estimates of
the high free energy barrier and to sample the transition state at the top of
the barrierBennett ; Torrie and Valleau (1977); Van Erp, Moroni, and Bolhuis
(2003); Allen, Frenkel, and ten Wolde (2006); Valsson, Tiwary, and Parrinello
(2016). These enhanced sampling techniques are used to determine the free
energy landscape by improving the sampling efficiency in poorly visited
regions of order parameter space and imposing a condition of equilibrium or
zero net flux in the order parameter space. Key to such numerical
investigations is a suitable model of the system of interest as well as a low-
dimensional representation with appropriately chosen collective variables, or
order parameters, that effectively distinguishes the relevant statesPeters
_et al._ (2013). The landscape as well as the probabilistic rate of the
transformation are closely related, with the set of activated processes –
where free energy barriers separate the states of interest – receiving
tremendous scientific interest over many years. To a first approximation, the
rate and the barrier height can be related through an Arrhenius-like equation
in the following way:

$k=Ae^{-\frac{\Delta G}{k_{B}T}}.$ (1)

Here, $k$ is the rate, $A$ is a kinetic pre-factor, $k_{B}$ is Boltzmann’s
constant, $T$ the temperature and $\Delta G$ the height of the barrier
measured from the initial or reactant state. A prominent concept in this
context is the mean first passage time, which is the inverse of the rate.
Kramers framed the progress of a reaction as a diffusive barrier crossing
where the time evolution of the collective variable of interest obeyed the
Smoluchowski equation Kramers (1940). A key assumption in this approach is
that degrees of freedom other than the collective variable relaxes on
timescales much shorter than the barrier crossing timescale, and can thus be
averaged out. The connection between the mean first passage time and the free
energy barrier can be directly exploited to provide estimates of the free
energy from unconstrained simulations where the collective variable evolves
from an initial value between a reflecting and an absorbing boundary, through
a kinetic reconstruction, developed by Reguera and co-workers Wedekind, Strey,
and Reguera (2007); Wedekind and Reguera (2008). Such an approach does not
require one to produce a condition of zero net flux and has been used in the
context of metastable-to-stable phase transitions such as nucleation in deeply
supercooled conditionsWedekind  _et al._ (2009); Lundrigan and Saika-Voivod
(2009); Thapar and Escobedo (2015); Goswami  _et al._ (2021). However, the
above approach, based on the exact relationship between the free energy
function and the mean first passage times, is available when one considers
only a single order parameter.

A number of researchers have noted the connection between equilibrium free
energies and rates and their non-equilibrium steady state counterparts in
driven or dissipative systems conditionsCrooks (1998); Banik, Chaudhuri, and
Ray (2000); Hummer and Szabo (2001). Framing the problem of free energy
calculations in the case where a constant flux to an absorbing state alters
the steady state sampling probability is relevant to the context of the
liquid-liquid phase transition, which has been vigorously investigated since
it was first proposed as an explanation for the thermodynamic anomalies
exhibited by water based on numerical evidencePoole  _et al._ (1992). Biased
simulations such as umbrella sampling have been extensively used to perform
such free energy calculations. However, the choice of order parameter and bias
protocol are key to obtaining meaningful results, with inappropriate choices
leading to qualitatively misleading resultsGoswami  _et al._ (2021). Thus, a
method to compute free energies from unconstrained simulations initialised
from the disordered liquid and proceeding till crystallisation occurs would be
of value both in the context of metastable network-forming liquids, as well as
numerous other contexts where secondary or tertiary order parameters relax on
comparable timescales to the primary order parameter separating metastable
states from the globally stable stateDeFever and Sarupria (2019); Verma and
Nair (2022).

In this work, we develop a methodology to reconstruct multi-dimensional free
energy landscapes from unconstrained simulations evolving between reflecting
and absorbing boundary conditions separated by a free energy barrier along a
primary collective variable. We focus on reconstruction in cases where
multiple metastable states exist, separated by a barrier along an orthogonal
collective variable. We begin by defining the phenomenological rate of
transformation from reactant, $A$, to product, $B$, as the ratio of two
partition functionsChandler (1978); Bolhuis  _et al._ (2002) weighted on paths
connecting $A$ and $B$ and all paths exiting $A$, respectively. We consider
the effect of including an additional absorbing condition and discuss
conditions under which the rate is not altered. However, the effective
positive fluxVan Erp and Bolhuis (2005) between different regions of order
parameter space and the steady state sampling, which are respectively related
to the two aforementioned partition functions, are. Using this, we are able to
relate the steady state sampling in the presence of the additional boundary
(through which there is a finite flux) to the underlying equilibrium sampling
in the flux-balanced condition in the absence of the additional absorbing
boundary. This relationship between steady state sampling and equilibrium
sampling enables an estimate of the free energy landscape from unconstrained
trajectories. We demonstrate the efficacy of this approach on a model of
independent overdamped Brownian random walkers on a potential energy surface
as well as on the reconstruction of the free energy surface for supercooled
Stillinger-Weber siliconStillinger and Weber (1985) from unconstrained
molecular dynamics trajectories. Results are compared with recently published
estimates of the free energy landscape for silicon from umbrella sampling
simulationsGoswami  _et al._ (2021); Goswami and Sastry (2022).

The paper is organised in the following way: in Section II we discuss how one
can extend the single order parameter free energy to multiple order parameters
under the assumption of Boltzmann sampling along the other order parameters.
In Section III we describe the systems on which we employ our approach, the
test system of independent random walkers on a potential energy landscape in
Section III.1, as well as liquid silicon in Section III.2. In Section IV we
describe briefly the mean first passage time (MFPT) method employed to
reconstruct one dimensional free energy profiles. In Section V, we describe
the reconstruction of free energies as a function of multiple order parameters
using the MFPT method and the assumption of Boltzmann sampling along the
second order parameter, which reveals the inadequacies of such an approach. In
Section VI we derive the relationship between steady state sampling and
equilibrium sampling in a multi-dimensional order parameter space, which is
our main result. In Section VII, we describe results on the test system
demonstrating the effectiveness of our approach. We then apply this method to
the more complex case of supercooled liquid silicon in Section VIII and
reconstruct the barrier to crystal nucleation as well as the barrier profile
along a second order parameter, density, revealing the presence of two liquid
states. In Section IX we discuss possible ways to improve on our approach,
focusing on a few key shortcomings before a brief discussion in Section X
summarising our findings and promising future directions.

##  II Extending the one order parameter free energy to multiple order
parameters

Here, the steps to obtain the two order parameter free energy $\beta\Delta
G(x,y)$ from the single order parameter free energy $\beta\Delta G(x)$ and the
sampling along two order parameters, $P(x,y)$, are described. In equilibrium
the sampling probability can be related to free energy differences in the
following way,

$P_{eq}(x,y)=Ae^{-\beta\Delta G(x,y)}$ (2)

Using the single order parameter sampling probability,

$P(x)=\int_{-\infty}^{\infty}P(x,y)dy$ (3)

we can write,

$P_{eq}(x)=\int_{-\infty}^{\infty}P_{eq}(x,y)dy=A\int_{-\infty}^{\infty}e^{-\beta\Delta
G(x,y)}dy=Ae^{-\beta\Delta G(x)}$ (4)

Using this relation between the free energy along $x$, $\beta\Delta G(x)$ and
$P_{eq}(x)$ to give

$P_{eq}(x,y)=P_{eq}(x)e^{\beta\Delta G(x)}e^{-\beta\Delta G(x,y)}$ (5)

From this, one can rearrange and to get

$\beta\Delta G(x,y)=\beta\Delta
G(x)-ln\left(\frac{P_{eq}(x,y)}{P_{eq}(x)}\right)$ (6)

Eq. 6 is a relation between free energy and equilibrium probabilities. Note
that if the relative weights of sampling different $y$ for a given $x$,
$P_{st}(y;x)$, are in equilibrium, then we can substitute $P_{eq}(x,y)$ with
the measured $P_{st}(y;x)$ in Eq. 6 to obtain $\beta\Delta G(x,y)$. We need to
define $P_{st}(x)=\int dyP_{st}(y;x)$ for the denominator in Eq. 6 to do so.
In this case $\beta\Delta G(x)$ needs to be obtained independently, from some
other method like umbrella sampling along $x$ or the kinetic reconstruction
from the mean first passage timeWedekind, Strey, and Reguera (2007); Wedekind
and Reguera (2008); Wedekind  _et al._ (2009). One can compare the measured
free energy from (say) the single order parameter reconstruction along $x$, ,
with the quantity $G(x)$ which is given by:

$e^{-\beta\Delta G(x)}=\int_{-\infty}^{\infty}e^{-\beta\Delta G(x,y)}dy$ (7)

Eq. 6 is the result that allows the extension of single order parameter free
energies to multiple order parameters. However, as we will see, this can be
used as-is only when sampling along the other order parameters is Boltzmann.
When this is not the case, corrections need to be used to obtain the correct
free energy, which are discussed in detail in Sec. VI. As we shall see in the
discussion that follows, one can obtain the free energy as a function of one
or multiple order parameters by employing the identified corrections.

##  III Systems studied

In this section we describe the systems on which we test this method. We first
test the numerical reconstruction procedure for the case of independent
overdamped random walkers on an energy landscape which mimic the behaviour of
trajectories in order parameter space for which the dynamics are not
discontinuous. We consider first a landscape having two metastable minima and
a globally stable minimum, motivated by the problem of reconstructing the free
energy landscape of liquids displaying polyamorphism. We also consider other
such test cases (see Appendix B) as well as supercooled liquid silicon.

###  III.1 Test system

The model potential energy function we consider is of the form in Eq. 8, which
is a sum of $4$ Gaussian functions in two dimensions (see Fig 1 for
illustration and Table 1 in Appendix A for values of constants). An additional
harmonic cost potential, $V_{c}(y)$ is applied to ensure that random walkers
sample order parameter space within $y\in[-0.5,0.5]$.

$V(x,y)=\sum\limits_{i=1}^{4}V_{i}(x,y;x_{i},y_{i},\sigma x_{i},\sigma
y_{i})+V_{c}(y)$ (8)

Figure 1: The model potential energy landscape for Eq. 8. Shown here for a
case where the barrier along $x$ is approximately $2~{}k_{B}T$. Contour lines
are $0.25~{}k_{B}T$ apart. The reflecting boundary is at $x=0$ and absorbing
boundary at $x=0.75$ as shown. In order to determine the accuracy of sampling
along $y$, we compare slices along different values of $x$, such as
$x~{}\in~{}[0.10,0.12]$ as marked in the figure.

The surface is specified such that there are two saddles, with two metastable
states separated from the globally stable state at large $x$ by the barrier
along $x$ at $x\approx 0.4$. Such a surface is relevant in contexts such as
liquid polymorphism, where multiple metastable liquid states co-exist with the
globally stable crystalline state in a number of anomalous model liquidsPalmer
_et al._ (2014); Smallenburg, Filion, and Sciortino (2014); Debenedetti,
Sciortino, and Zerze (2020); Goswami and Sastry (2022). The flux to the
globally stable state at $x>0.4$ is controlled by the height of the scaled
barrier (in units of $k_{B}T$) along $x$. We study cases where the height of
the scaled barrier is low, $2~{}k_{B}T$ to $6~{}k_{B}T$, where the flux across
the barrier along $x$ is high, leading to larger deviations of steady state
sampling from equilibrium sampling. This is also consistent with our
expectation for deeply supercooled liquids where the barrier to
crystallisation is found to be very lowGoswami  _et al._ (2021). $N_{traj}$
non-interacting random walkers are initialised at $x=0$ and different $y$,
either sampling the Boltzmann distribution or at a specified $y$ value at
$x=0$. A reflecting boundary condition is placed at $x=0$ and an absorbing
boundary condition at $x=x_{c}=0.75$ for all $y$. For the reflecting boundary
condition at $x=0$, if a trial move places a particle at
$x^{{}^{\prime}}_{new}<0$, the trial move is modified to
$x_{new}=-x^{{}^{\prime}}_{new}$ and accepted or rejected depending on the
Boltzmann-weighted energy at $\beta\Delta V(x_{new},y_{new})$. No boundary
conditions are imposed along the $y-axis$. Each MC sweep consists of
$N_{traj}$ trial displacements of the random walkers. In each trial
displacement, a random walker is chosen with uniform probability and is
displaced by ${dx,dy}\in[-\delta:+\delta,-\delta:+\delta]$ Here, the value of
$\delta=4\times 10^{-3}$ is used while the order parameter space is divided
into equal-sized square bins of size $0.02$. While using larger step sizes,
i.e., comparable to the bin size, introduces sampling issues, we have
determined that the chosen step size does not affect our results. Trial
displacements are accepted or rejected using a Boltzmann weight for the change
in energy for every trial move. One can thus obtain the steady state sampling
probability, $P_{st}(x,y)$, from a number of trajectories that proceed to an
absorbing boundary condition. We demonstrate free energy reconstruction using
$N_{traj}=600$ such independent trajectories. We also compare results for the
reconstruction of the barrier along $x$ with those obtained for the same set
of trajectories by a kinetic reconstruction using the MFPT developed by
Reguera and co-workersWedekind, Strey, and Reguera (2007); Wedekind and
Reguera (2008). The procedure for this method and results for the single order
parameter problem are discussed in the next section, following which we
discuss the two order parameter reconstruction using Eq. 6.

###  III.2 Supercooled silicon

The other system we consider is liquid silicon modelled by the Stillinger-
Weber potentialStillinger and Weber (1985). The existence of two metastable
liquid states for this model has been investigated intenselySastry and Angell
(2003); Ganesh and Widom (2009); Beye  _et al._ (2010); Vasisht, Saw, and
Sastry (2011) with recent free energy calculations also finding a high density
liquid and a low density liquid separated by a free energy barrierGoswami and
Sastry (2022). This scenario is analogous to other network-forming liquids
such as waterPalmer  _et al._ (2014); Debenedetti, Sciortino, and Zerze
(2020), silicaChen, Lascaris, and Palmer (2017) and patchy colloidal model
liquidsSmallenburg, Filion, and Sciortino (2014); Neophytou, Chakrabarti, and
Sciortino (2022) where two liquid states have been identified. In the case of
water, silica and silicon, the globally stable crystalline state is separated
from two metastable liquid states by the free energy barrier to
crystallisation. In order to reconstruct the free energy from unconstrained
molecular dynamics simulations, we initialise $400$ independent molecular
dynamics simulations from configurations of randomly placed particles without
overlap at a density of $2.48~{}gcc^{-1}$. Molecular dynamics simulations are
performed in the isothermal-isobaric ensemble using the LAMMPS packagePlimpton
(1995) at target pressures and temperatures of $P=0.75~{}GPa$ and
$T=975K,~{}985K,~{}995K$, monitoring the size of the largest crystalline
clusterRomano, Sanz, and Sciortino (2011); Goswami  _et al._ (2021), denoted
$n_{max}$, and the density ($\rho$), with simulations being extended till a
largest cluster size of $n_{max}=80$ being reached. Trajectories in
($n_{max}$,$\rho$) space are treated as random walks on the underlying free
energy surface.

##  IV Kinetic reconstruction along one order parameter using the mean first
passage time

We use Eq. 9 and Eq. 10 as described inWedekind, Strey, and Reguera (2007);
Wedekind and Reguera (2008); Wedekind  _et al._ (2009, 2015) to obtain the 1D
barrier along $x$ from a set of unconstrained trajectories that proceed until
the absorbing boundary at $x=0.75$ is reached.

$\beta\Delta G(x)=\beta\Delta
G(x=1)+ln\left(\frac{B(x)}{B(1)}\right)-\int_{1}^{x}\frac{dx^{\prime}}{B(x^{\prime})}$
(9)
$B(x)=-\frac{1}{P_{st}(x)}\left[\int_{x}^{b}P_{st}(x^{\prime})dx^{\prime}-\frac{\tau(b)-\tau(x)}{\tau(b)}\right]$
(10)

Eq. 9 and Eq. 10 are the equations used to reconstruct the free energy from
the MFPT and the steady state probability.

Figure 2: Mean first passage time, $\tau_{MFPT}(x)$ and the steady state
probability $P_{st}(x)$. Absorbing boundary condition at $x=0.75$, reflecting
boundary condition at $x=0$. Number of walkers is $N_{traj}=600$, sufficient
to generate smooth data.

(a)  (b)  (c)

Figure 3: The free energy along $x$, $\Delta G(x)$ obtained from Eq. 9. For
reference, the potential dependence is shown by integrating out the $y$
variation (Eq. 7). The number of walkers is $N_{traj}=600$ and the absorbing
boundary condition is placed at $x=0.75$. Random walks are performed at a
temperature of $T=0.08$ on $3$ landscapes with barrier heights along $x$ of
$1-2~{}k_{B}T$ (panel (a)), $3~{}k_{B}T$ (panel (b)) and $5~{}k_{B}T$ (panel
(c)). Error is minimised by shifting the curves to minimise the difference
$|\beta\Delta G(x)-V_{p}(x)|$.

Fig. 2 show the measured mean first passage time and the steady state
probability along $x$ for $3$ scaled barrier heights along $x$. The
reconstructed free energy using Eq. 9 and Eq. 10 is shown in Fig. 3, compared
with the expected curve integrating out the $y$ dependence, given in Eq. 7.

##  V Sampling along orthogonal order parameters: Deviation for high flux
through the absorbing boundary

We next compare the reconstructed free energy along the $y$ direction to the
corresponding cross-section of the potential along $y$ using Eq. 6. The height
of the barrier is controlled by modulating the potential (see Appendix A). We
find that as the barrier along $x$ is lowered, the deviation of steady state
sampling from the target, $V(y)_{x=0.01}$, increases. This is shown in Fig. 4.
This can be rationalised as the enhanced flux across the lower barriers
driving the system away from equilibrium sampling to a non-equilibrium steady
state.

Figure 4: Slices along $y$ of the negative log of the steady state
probability, $-ln[P_{st}(y;x_{0})]$, compared with the target, underlying
landscape, $V(y;x_{0})$ at $x_{0}~{}\in[0.10,0.12]$. $N_{traj}=600$ random
walkers are initialised at $x=0$ and at different $y$ proportional to the
Boltzmann weight on the target distribution (panel (a)) or at a point in the
basin $A$ (panel (b)). As the scaled barrier along $x$ is lowered, the
deviation in steady state sampling of $y$ from the target distribution
increases as shown in panel (a). In panel (b) we observe that for a fixed
value of the scaled barrier along $x$, the deviation of the steady state
distribution from the target distribution is higher when the initialisation is
at a single point along $y$.

The errors in the sampling along the orthogonal order parameter, $y$, that are
summarised in Fig. 4 arise from compounding factors: Firstly, when the scaled
barrier along $x$ is low, the net flux to the absorbing boundary prevents
Boltzmann sampling along $y$ for $x<x^{*}$ (where $x^{*}$ is location of the
barrier). The final steady state distribution also has a dependence on the
injection probability, as seen in panel (b) of Fig. 4 where the deviation is
higher for the point injection compared to the Boltzmann weighted injection
along $y$. We next discuss how to correct for these errors by identifying a
relationship between steady state sampling and equilibrium sampling.

##  VI Relating steady state sampling to equilibrium sampling

We found in the previous section that sampling in the presence of a low
barrier along $x$ and a constant non-zero flux across it, that the measured
sampling along $y$ deviates from the underlying Boltzmann distribution. This
steady state is achieved by re-injecting or restarting a trajectory from an
injection point between the reflecting and the absorbing condition and
tracking it until it crosses the absorbing boundary, whereupon another
trajectory is started, thus conserving the number of “active" trajectories at
any point of time. In this section we will discuss our approach to correct for
the systematic deviation in sampling by considering the effect of the steady
state flux across the absorbing boundary that is established. In order to
understand this, we begin by considering the rate(s) of traversal between any
two regions of the order parameter space, $A$ and $B$, along the lines of the
development of transition path sampling, transition interface sampling and
forward flux sampling methodsChandler (1978); Dellago  _et al._ (1998);
Bolhuis  _et al._ (2002); Van Erp and Bolhuis (2005); Allen, Valeriani, and
Ten Wolde (2009); Vanden-Eijnden  _et al._ (2010). We represent the order
parameter values for simplicity as $x$, and those of $A$ and $B$ by $x_{A}$
and $x_{B}$. We consider $A$ and $B$ to be points within the grid resolution
we specify, but $x_{A}$ and $x_{B}$ can equivalently be considered to be a set
of $x$ and $y$ values. We first write functions $h_{A}(x)$ and $h_{B}(x)$
where $h_{A}(x)=1~{}if~{}x\in x_{A}$ and $0$ otherwise, and
$h_{B}(x)=1~{}if~{}x\in x_{B}$ and $0$ otherwise. The phenomenological rate of
transition from non-intersecting regions of order parameter space, reactant
$A$ and product $B$, is then given in terms of the time correlation of the
product of these functions asChandler (1978); Frenkel and Smit (2001); Bolhuis
_et al._ (2002); Van Erp and Bolhuis (2005); Allen, Valeriani, and Ten Wolde
(2009).

$k_{AB}=\frac{d}{dt}C(t)=\frac{d}{dt}\left[\frac{\langle
h_{A}(x_{0})h_{B}(x_{t})\rangle}{\langle h_{A}(x_{0})\rangle}\right]$ (11)

The assumption of a time-invariant rate, related to the inverse mean first
passage time, is typically invoked in the context of regions $A$ and $B$ that
are separated by a barrier with a steady state rate across it that is
established after an initial transient and decays on a global reaction
timescaleVan Kampen (1992); Frenkel and Smit (2001). Here, we consider a
steady state scenario where the rate of traversal between any two non-
overlapping regions of order parameter space is of interest. The average here
is over an ensemble of trajectories and weighted on the probability of
observing a path connecting $A$ and $B$, of length $t$ , denoted
$P[\\{x_{t}\\}]$. This is the probability of observing a trajectory, i.e., the
sequence $\\{x_{0},x_{1},\dots,x_{t}\\}$ For stochastic trajectories, with
transition matrix $\mathbf{T}$,

$P[\\{x_{t}\\}]=P(x_{0})\prod\limits_{0<t^{\prime}\leq
t}\mathbf{T}_{x_{t^{\prime}-1}x_{t^{\prime}}}$ (12)

For the deterministic case the initial conditions fully specify the
probability of observing a path, $P[\\{x_{t}\\}]=P(x_{0})$. In the equilibrium
case, $\langle h_{A}(x)\rangle$ is equal to the equilibrium probability of
being in $A$. Upon the introduction of an additional absorbing boundary $C$,
one expects both the sampling probabilities as well as the probability of
observing a given trajectory are altered. The probability of observing a path
is altered through an alteration of the microscopic transition matrix
$\mathbf{T}$, at some point $r$ in the vicinity of the absorbing boundary $C$.
We can write that outward transition probabilities for a state neighbouring
the new absorbing state, labelled $r$, are altered by the introduction of the
new absorbing state $C$. For simplicity, we consider a single such state $r$;
this choice should not affect our conclusions. The probability or weight of a
path labelled $\\{x_{t}\\}$ becomes:

$\displaystyle P^{\prime}[\\{x_{t}\\}]=$
$\displaystyle\prod\limits_{1<t^{\prime}\leq
t}\left[\mathbf{T}_{x_{t^{\prime}-1}x_{t^{\prime}}}+\delta_{r,t^{\prime}-1}(\mathbf{T^{\prime}}_{x_{t^{\prime}-1}x_{t^{\prime}}}-\mathbf{T}_{x_{t^{\prime}-1}x_{t^{\prime}}})\right]$
$\displaystyle\times f(x_{0})P(x_{0})$ (13)

Intuitively, the integral over all paths is changed when the fraction of paths
between $A$ and $B$ that pass through $r$ is significant. Otherwise the term
in the product remains unchanged. For this to be true, $A$ and $B$ should both
be far from $C$, such that typical paths connecting them are significantly
shorter than the typical length of a full trajectory that proceeds until it
encounters $C$. Moreover, if the region $A$ is defined such that $f(x_{0})$ is
a constant value $f$, within $A$, then it can be ignored while evaluating the
two integrals in Eq. 11. Under these two conditions, we then assume that the
phenomenological rates of traversal between $A$ and $B$, $k_{AB}$ and
$k_{BA}$, are unaltered upon the addition of the absorbing state $C$. However,
the flux between the two, as well as the steady state sampling probabilities
are altered with respect to the corresponding equilibrium states. The flux
from $A$ to $B$, $\langle\Phi_{AB}\rangle$, is defined here as the number of
trajectories entering $B$ in a given time window $[t,t+\Delta t]$ that had
their origin in $A$ at $t=0$. The flux per unit time is obtained by dividing
it by the length of the interval, $\Delta t$Van Erp, Moroni, and Bolhuis
(2003). In the steady state condition, we can exploit the following
replacement, $h_{A}(x_{0})h_{B}(x_{t})=h_{A}(x_{-t})h_{B}(x_{0})$, whose time
derivative can then be evaluated at $t=0$. Further, the time dependence for
$h_{A}$ can also be dropped, given the steady state condition, so long as the
trajectory did not visit $B$ prior to $t=0$ (to ensure that a trajectory that
leaves $A$ once only counts towards the flux entering $B$ once)Van Erp,
Moroni, and Bolhuis (2003); Van Erp and Bolhuis (2005). Thus, the steady state
flux, $\langle\Phi_{AB}\rangle=\langle
d/dt|_{t=0}h_{A}(x_{-t})h_{B}(x_{0})\rangle$ can be evaluated as the number of
trajectories entering $B$ during some interval $\Delta t$, who were last in
$A$ before $B$, and can be aggregated over a given time interval. In the
equilibrium case, the fluxes are balanced, and any random walker (or
trajectory) that visits either $A$ or $B$, visits the other as well. The
introduction of an additional absorbing boundary at $C$ introduces the added
condition that only random walkers leaving $B$ ($A$) that reach $A$ ($B$)
before reaching $C$ contribute to the flux, termed the splitting
probabilityMüller, Talkner, and Reimann (1997); Huang and Chen (2021) in the
non-equilibrium steady state case. Thus, by matching rate of injection at $A$
to the rate of first passage at $C$, one obtains a steady state characterised
by the macroscopic flux through $C$. One can then measure the number of
trajectories that make a transition from $A$ to $B$, or the reverse, before
reaching $C$, in the time it takes for $N$ such trajectories to traverse from
the injection point to $C$. This gives the flux subject to the steady state
specified by our injection rate. $\langle h_{A}\rangle$ is the probability
that a randomly chosen starting point for a trajectory is in $A$, which is
also the steady state sampling probability for $A$. Using this, we are now
able to write the phenomenological rate as:

$k_{AB}=\frac{\langle\Phi_{AB}\rangle}{\langle h_{A}\rangle}.$ (14)

We now discuss how Eq. 14 can be used to relate the steady state sampling
probability to the equilibrium sampling probability. In equilibrium,
$\langle\Phi_{AB}\rangle_{eq}=\langle\Phi_{BA}\rangle_{eq}$ (zero current),
$\langle h_{A}\rangle_{eq}=P_{eq}(A)$, giving the detailed balance condition

$k_{AB}P_{eq}(A)=\langle\Phi_{AB}\rangle_{eq}=\langle\Phi_{BA}\rangle_{eq}=k_{BA}P_{eq}(B)$
(15)

Upon addition of the absorbing boundary, $C$, trajectories exiting $A$ (or
$B$) can now be terminated at $C$. In the resulting steady state condition,
the probability that a randomly chosen starting point is in $A$ is now
$\langle h_{A}\rangle_{st}\equiv P_{st}(A)$, altered from $\langle
h_{A}\rangle_{eq}$. In steady state, the fluxes $\langle\Phi_{AB}\rangle_{st}$
and $\langle\Phi_{BA}\rangle_{st}$ are not equal. In order to relate the
steady state quantities to the equilibrium quantities, we first assume
trajectories cannot be initiated at $C$ (the new absorbing state). The ratio
of flux per unit time to sampling probability gives us the (assumed) unaltered
rate. Thus, in steady state, where
$\langle\Phi_{AB}\rangle_{st}\neq\langle\Phi_{BA}\rangle_{st}$

$\displaystyle\frac{P_{st}(A)}{P_{st}(B)}$ $\displaystyle=\frac{\langle
h_{A}\rangle_{st}}{\langle h_{B}\rangle_{st}}$
$\displaystyle=\frac{\langle\Phi_{AB}\rangle_{st}}{\langle\Phi_{BA}\rangle_{st}}\frac{k_{BA}}{k_{AB}}=\frac{\langle\Phi_{AB}\rangle_{st}}{\langle\Phi_{BA}\rangle_{st}}\frac{P_{eq}(A)}{P_{eq}(B)}$
(16)

What we want is to infer the equilibrium sampling probability from the
measured steady state sampling probability. It is helpful to then re-write the
equation above as

$\frac{P_{eq}(B)}{P_{eq}(A)}=\frac{\langle\Phi_{AB}\rangle_{st}}{\langle\Phi_{BA}\rangle_{st}}\frac{P_{st}(B)}{P_{st}(A)}.$
(17)

The quantities on the right hand side are evaluated from numerical
simulations, for all $B$ of interest, which them results in an estimate of the
free energies relative to that at $A$.

###  VI.1 Notes on numerical implementation

We will drop the $\langle\rangle$ hereafter, when describing
$\langle\Phi_{AB}\rangle$. We simulate $N_{traj}$ trajectories, injected at
$x=0$ in the basin marked $A$ in Fig. 1, that proceed to the absorbing bound,
marked $C$ in Fig. 1, through a random walk on the potential surface $V(x,y)$.
We then obtain $P_{st}(x,y)$ from the cumulative number of times each
trajectory visits each $(x,y)$ bin. The relevant normalisation factor is the
cumulative length of the $N_{traj}$ trajectories. We also compute the fluxes
to and from the injection point, $A$, and every other $(x,y)$ in order
parameter space, denoted $B$ in Eq. 17. It is important to note that the
established steady state is subject to the injection point, $A$, and the
definition of a trajectory, which is initialised at the injection point and is
terminated at the absorbing boundary, $x_{c}$. The rate of injection is chosen
to match the rate of termination in the cases we discuss.

We count the number of trajectories that, having visited $A$ at some time $t$,
subsequently visit $B$ (a given $(x,y)$ bin) at some later time
$t^{\prime}>t$. Likewise, we count each trajectory that, having visited a
given $B$ at some time $t$, subsequently visits $A$ at a time $t^{\prime}>t$
before being terminated at the absorbing boundary. We consider this count as
providing estimates of $\Phi_{A\rightarrow B}$ and $\Phi_{B\rightarrow A}$,
upon division by the total trajectory lenght, which we need not explicitly
consider since we are only interested in the ratio $\Phi_{A\rightarrow
B}/\Phi_{B\rightarrow A}$. Note that $\Phi_{B\rightarrow A}$ is sampled poorly
beyond the barrier and far enough beyond the barrier, this quantity goes to
$0$. The ratio $\hat{\pi}(x,y)=\frac{\Phi_{A\rightarrow
xy}}{\Phi_{xy\rightarrow A}}$ is what we apply as a correction factor, in
order to obtain the free energies, using Eq. 17.

$\frac{P_{eq}(x,y)}{P_{eq}(A)}=\hat{\pi}(x,y)\frac{P_{st}(x,y)}{P_{st}(A)}$
(18)

Inferring $P_{eq}(x,y)$ from $P_{st}(x,y)$ as shown earlier allows us to use
Eq. 6 which we initially arrived at as the extension of the single order
parameter free energy estimate to multiple order parameters. This is done in
the following way:

$\displaystyle\beta\Delta G(x,y)=$ $\displaystyle\beta\Delta
G(x)-ln\left(\frac{P_{eq}(x,y)}{P_{eq}(x)}\right)$ $\displaystyle=$
$\displaystyle\beta\Delta
G(x)-ln\left(\frac{P_{eq}(x,y)P_{eq}(A)}{P_{eq}(A)P_{eq}(x)}\right)$
$\displaystyle=$ $\displaystyle\beta\Delta
G(x)-ln\left(\frac{\hat{\pi}(x,y)P_{st}(x,y)P_{eq}(A)}{P_{st}(A)P_{eq}(x)}\right)$
$\displaystyle\beta\Delta G(x,y)=$ $\displaystyle\beta\Delta
G(x)-ln\left(\frac{\hat{\pi}(x,y)P_{st}(x,y)}{P_{eq}(x)}\right)+const.$ (19)

Above we have written terms dependent only on $A$ as an irrelevant constant.
$P_{eq}(x)$ is defined as

$\displaystyle P_{eq}(x)$
$\displaystyle=\int\limits_{-\infty}^{\infty}P_{eq}(x,y)dy$
$\displaystyle=\frac{P_{eq}(A)}{P_{st}(A)}\int\limits_{-\infty}^{\infty}\hat{\pi}(x,y)P_{st}(x,y)dy,$
(20)

using Eq. 18. $P_{eq}(x)$ is thus obtained upto a multiplicative constant
integrating out the $y$-dependence of $\hat{\pi}(x,y)P_{st}(x,y)$. We obtain
$\beta\Delta G(x)$ independently and use Eq. 18 to obtain the multi-
dimensional free energy $\Delta G(x,y)$.

Other schemes can be developed to define pairs for which we can apply Eq. 19.
Results shown in the subsequent sections consider the definition of $A$ used
above, as the point of injection of the trajectories. For the case where
trajectories are injected at $x=0$ with a Boltzmann-weighted injection
probability at different $y$, the same point is chosen as the state $A$ as in
the point-injection case (the basin $A$ in Fig. 1).

###  VI.2 Similar results in the literature

The treatment of the phenomenological rate employed here has established
itself as an immensely useful concept in numerical rate calculations and also
free energy calculations in concert with milestoning, path and interface
samplingBolhuis  _et al._ (2002); Van Erp and Bolhuis (2005); Allen,
Valeriani, and Ten Wolde (2009). An analysis of trajectory segments in
partial-path transition interface sampling (PPTIS) can be related to the
single order parameter equivalent of the result hereQin, Dellago, and
Kozeschnik (2019). In this method also, backward fluxes beyond the barrier are
not easily obtained. Typical interface sampling and methods that enhance
fluxes rely on creating a situation of equal forward and backward flux to
achieve equilibrium samplingValeriani  _et al._ (2007); Qin, Dellago, and
Kozeschnik (2019). In this work, we find that one can use a finite, but
unequal, backward flux to infer equilibrium sampling probability from the
measured steady state sampling probability. Such a situation arises naturally
where simulations proceed to an absorbing boundary. Forward flux sampling has
also been combined with the mean first passage time, eliminating the need for
backward trajectories and zero net fluxThapar and Escobedo (2015). Systematic
alterations in the phenomenological rate have been investigated using a
treatment of the rate expression that considers the effect of an additional
field on the path partition functions, rather than altered boundary
conditionsKuznets-Speck and Limmer (2021). Studies of first passage times and
first passage probabilities for Markov processes with specific boundary
conditions have discussed related concepts such as the splitting
probabilityMüller, Talkner, and Reimann (1997). Of particular note is the
study of random walks with stochastic resetting, where the effect of resetting
on the mean first passage time and the rates has been investigatedEvans,
Majumdar, and Schehr (2020); Huang and Chen (2021); Chen, Li, and Huang
(2022). Generalisations of Kramers’ formalism to open or driven systems in one
dimension have also been described, with similar ideasBanik, Chaudhuri, and
Ray (2000); Hummer and Szabo (2001).

##  VII Results for the test system

Writing

$\frac{P_{eq}(x,y)}{P_{eq}(A)}=\frac{\Phi_{A\rightarrow
xy}}{\Phi_{xy\rightarrow A}}\frac{P_{st}(x,y)}{P_{st}(A)},$ (21)

one obtains a correction factor that works for either a Boltzmann
initialisation or a point injection at some $y_{0}$, $x=0$. In Fig. 5, we show
the reconstructed (single order parameter) free energy surfaces along $x$ and
along $y$. We show the full free energy surface reconstruction and a
comparison with errors in Fig. 6.

Figure 5: The free energy reconstruction tested at $T=0.08$ for different
barrier heights for $N_{traj}=600$ random walkers with a fixed maximum step
size of $0.004$. In panel (a) the reconstruction along $x$ is shown for two
types of injection protocols as described in the legend. For the point
injection protocol, $x=0$, $y=0.35$ in the $A$ basin is chosen. The results
are compared with the kinetic reconstruction using the MFPT and against the
effective barrier along $x$, $V_{p}(x)$ in Eq. 7. In panel () and (c) a slice
is taken along $x_{0}=0.11$ to compare with the “target" slice $V(y;x_{0})$
for different barrier heights along $x$ ($\Delta G(x^{*})=1~{}k_{B}T$ in panel
(b) and $\Delta G(x^{*})=5~{}k_{B}T$ in panel (c)). We find that for both the
reconstruction along $x$ and along both $x$ and $y$, accurate reconstructions
are possible regardless of the distribution of initial $y$ values. Figure 6:
The full two-order parameter reconstruction for a case with a scaled barrier
$1~{}k_{B}T$ barrier along $x$ (see Appendix A for details on controlling
barrier height), as shown in the original potential energy surface in panel
(a. $N_{traj}=600$ random walkers are injected at $x=0$, $y=0.35$ and evolved
with a temperature of $T=0.08$, maximum step size of $0.004$ and histograms
binned with a square bin size of $0.02$. Panel (b) shows the reconstruction of
$\beta\Delta G(x,y)$ using Eq. 19 while panel (c) shows a reconstruction using
Eq. 18 $\beta\Delta G(x,y)$. The error between the reconstructed free energy
and the supplied potential energy surface is measured as $|\beta\Delta
G(x,y)-\beta V(x,y)|$ and shown in units of $k_{B}T$ in panel (d). The errors
are of the order of $<0.1~{}k_{B}T$ for $x<0.4$, which is approximately the
value of $x^{*}$.

The reconstruction of the single and two order parameter free energies using
both Eq. 18 (Fig. 5 (a) and Fig. 6 (b)) and Eq. 19 (Fig. 5 (a), (b), (c) and
Fig. 6 (c) , (d)) demonstrate the utility of this approach. Estimates of the
free energies of the metastable states and the barrier heights agree
quantitatively with the reference landscape. The errors in Fig. 6 (d) are
$<0.1~{}k_{B}T$ for all $(x,y)$ in the metastable regime and are thus of the
order of $1\%$.

The high errors beyond the barrier along $x$ arise from poor sampling of the
backward flux back to the region defined as source ($A$ at $x=0$, $y=0.35$ or
$x=0$ and all $y$ depending on injection protocol). In Section IX we describe
approaches to improve estimates beyond the barrier.

We also reconstruct different free energy surfaces, having either more
metastable states or multiple possible reaction pathways, to test the
generality of our scheme. Results for these alternate potential surfaces are
shown in Appendix B. We next describe results for the reconstruction of the
free energy landscape of supercooled silicon from unconstrained molecular
dynamics simulations using this methodology.

##  VIII Results for supercooled liquid silicon

We apply the methodology described above to the case of liquid silicon using
the size of the largest crystalline cluster, $n_{max}$ and the density,
$\rho$, as the order parameters with respect to which we reconstruct the free
energy. $n_{max}$ is analogous to $x$ in the test system and $\rho$ thus
corresponds to $y$. For each trajectory, labelled superscript $i$, we update
the steady state sampling frequency of the values of $n_{max}$ and $\rho$
sampled by it using
$N^{i}_{st}(n_{max},\rho)=\sum\limits_{t=0}^{t^{i}_{final}}\delta(n_{max}^{i}(t)-n_{max})\delta(\rho^{i}(t)-\rho)$.
The steady state sampling probability $P_{st}(n_{max},\rho)$ is obtained by
explicitly normalising with the sum of $N^{i}_{st}(n_{max},\rho)$ over all
$n_{max}$ and $\rho$.

The flux count is measured in the following way. First, we define the
injection point of $n_{max}~{}<~{}1$ and $2.45~{}\leq~{}\rho~{}<2.46$ as the
“source" or reference state $A$. For each trajectory, at time step $t$, we
consider the $n_{max}$ and $\rho$ values. If these are outside the region $A$,
we trace back along the trajectory to check if $A$ was visited before this
point. If it was visited (strictly always true for each $(n_{max},\rho)$), we
update the count of the flux from $A$ to the given $(n_{max},\rho)$ by $1$. It
must be ensured that multiple crossings from $A$ to a given $(n_{max},\rho)$
are not counted multiply. Likewise, at every time step at which the region $A$
is reached, we trace back along the trajectory and update the flux count from
every $(n_{max},\rho)$ that was visited prior to the given timestep and which
was not counted already.

We first compare estimates $\beta\Delta G(n_{max})$ obtained using a single
order parameter version of Eq. 18 with those obtained from the MFPT method in
Fig. 7. The small $n_{max}$ free energies are obtained by matching
$\beta\Delta G(n_{max})$ with $-ln(P_{st}(n))$ from unconstrained MD runs for
small $n$ (or $n_{max}$)Goswami  _et al._ (2021). $\beta\Delta G(n_{max})$
shows an artificial minimum, which is rectified from this comparison, as
detailed in Goswami  _et al._ (2021). We note that the results from the
present method compare rather well with those of the MFPT method.

We obtain the free energy in terms of the density $\rho$, $\beta\Delta
G(\rho)$, upto an irrelevant additive constant from the full probability
distribution $P(n_{max},\rho)$ by using

$P(\rho)=\sum\limits_{n_{max}=0}^{n_{max}=4}P(n_{max},\rho)$ (22)

and taking the negative logarithm. In Fig. 8 the reconstructed free energy as
a function of density, $\beta\Delta G(\rho)$ is obtained by using Eq. 19 to
reconstruct the free energy surface and Eq. 22 to get $\beta\Delta G(\rho)$.
These results are compared with corresponding results from umbrella sampling
runs constraining both $\rho$ and $n_{max}$ for which the data is obtained
from Ref. Goswami and Sastry, 2022. The density profiles show a shift in the
location of the metastable minimum in density from a high value of
$2.45~{}gcc^{-1}$ to a low density of $2.35~{}gcc^{-1}$ when the temperature
is changed from $T=995~{}K$ to $T=985~{}K$ at $P=0.75~{}GPa$ with the other
liquid state losing metastability at or around $T=985~{}K$. At $T=975~{}K$,
one observes a larger difference in the estimates for $\beta\Delta G(\rho)$
obtained using the two methods. This remains to be fully understood, possible
reasons being poor sampling of the high density liquid in the case of umbrella
sampling, given the high local variation observed, or limitations arising from
the assumption of invariant rate in writing Eq. 17.

In the next section, we discuss possible ways to improve on the methodology,
addressing the shortcomings of poor estimates beyond the barrier and the
possible sources of discrepancy in the results for silicon.

##  IX Strategies for improved free energy estimates

Given states $A$ and $B$ between which we want to measure the flux,
$\langle\Phi_{AB}\rangle$ and $\langle\Phi_{BA}\rangle$, the quality of
reconstruction is determined both by the extent of sampling the steady state
probability as well as the two fluxes. In this section we discuss approaches
to improve the sampling of order parameter space and therefore the resulting
free energy estimates by addressing these requirements.

### Using interfaces for accurate flux calulation

We first describe how to improve the sampling of the backward flux from $B$ to
$A$ for regions $B$ that are beyond the barrier along $x$, such that this
backward flux is low and therefore poorly sampled with a finite number of
trajectories. Better estimates of this flux can be obtained by placing an
interface between $A$ and $B$ and expressing the total flux as a product. This
approach can be used if and only if every trajectory from $A$ to $B$ and $B$
to $A$ passes through an intermediate, $I$, (different from the new absorbing
condition $C$). We can then write the following

$\displaystyle\Phi_{A\rightarrow B}$ $\displaystyle=\Phi_{A\rightarrow
I}\times\Phi_{I\rightarrow B}$ $\displaystyle\Phi_{B\rightarrow A}$
$\displaystyle=\Phi_{B\rightarrow I}\times\Phi_{I\rightarrow A}$

Choosing $I$ as a hyper-plane separating $A$ and $B$ ensures that this
condition is met. In the 2D case, $I$ is a line. We proceed by testing if
estimates can be improved for $x$ beyond the barrier, from where the flux back
to $A$ may be negligibly small. We place the line at $x=0.5$, beyond the
saddle. As a test, we can compare and check if the following equation is true

$\phi_{A\rightarrow xy}=\phi_{A\rightarrow I}\times\phi_{I\rightarrow xy},$
(23)

for each $(x,y)$ with $x>0.5$. This is easy to verify because there is a large
direct flux from $A$ to $(x,y)$ beyond the barrier. We find that the free
energy reconstruction is improved beyond $x=0.5$ by using this expression for
the flux. Fig. 9 (c) and (d) show the results from this procedure and can be
compared with panels (a) and (b) respectively in Fig. 9, which are obtained
without resolving the flux along the lines in Eq. 23. It should be noted that
for $(x,y)$ beyond the barrier and close to the absorbing state $C$, our
assumption that typical paths from $A$ to $xy$ and the reverse do not pass
through the neighbourhood of $C$ does not hold. As a consequence, the
phenomenological rate changes significantly when one or both of $A$ and $B$
are close to $C$, possibly also lacking a timescale over which its value is a
plateau.

### Sectioned reconstruction with different reference states

We next attempt to improve the free energy estimates by defining multiple
sections of the order parameter space, each with a unique reference $A$ state.
The equilibrium probability $P_{eq}(x,y)$ is obtained using Eq. 18
independently in each section using the fluxes with respect to the unique
reference state. We ensure that each reference state is also contained within
another neighbouring region so that the estimates for each reference region
can be matched with the estimate from the neighbouring region by a simple
shifting of $-ln(P_{eq}(A))$. This then gives a better estimate of fluxes
locally and the different free energy estimates are then shifted to
reconstruct the full surface. In Fig. 9 (e) and (f), we demonstrate the
results of this procedure for a case where the order parameter space is
divided into a $3\times 3$ grid. The bottom left corner of each region is
chosen as the reference state for it. The region sizes are larger than the
separation between the reference points so that there is overlap. The errors,
shown in panel (f), are significantly lower at large $x$ values than the
corresponding errors in panels (b) and (d), where the reconstruction was
performed using Eq. 18 and Eq. 23 respectively.

##  X Discussion

We have described an effective and efficient method to obtain estimates of
free energies as a function of multiple order parameters from unconstrained
simulations. With our motivation arising from the study of polymorphism in
supercooled liquid silicon, we address the problem of reconstructing a multi-
dimensional free energy surface that can distinguish the possible metastable
states as well as the globally stable crystalline state. In order for
simulation lengths to be tractable, an absorbing condition is placed at large
values of the crystallinity order parameter. This absorbing boundary
introduces a flux in the system altering sampling. By treating the
trajectories in the order parameter space as obeying Brownian motion in the
high-friction limit one can obtain the steady state sampling in order
parameter space, as well as the various point-to-point fluxes. We consider the
expression derived for the transition rate, expressed in terms of sampling
probability and fluxes, to determine a relation between steady state sampling
in the presence of a non-zero constant flux (to the imposed absorbing
boundary) to the underlying equilibrium sampling, which can otherwise only be
achieved under conditions of zero net flux, which describes detailed balance.
This relationship between the steady state sampling and the underlying
equilibrium sampling is the key aspect of our method, which allows us to
obtain the free energies without the added effort of ensuring the zero flux
condition. We show that the method works efficiently for multiple cases by
testing it on a test system of random walkers on a potential energy landscape.
We choose potential energy landscapes that have, in addition to the globally
stable state, (i) multiple metastable states, (ii) metastable states that are
not on the primary “transition tube" connecting the initial metastable state
to the final, globally stable state (see Appendix B Fig. 10), and (iii),
multiple paths to the globally stable state with multiple saddles of different
heights connecting the metastable states to the globally stable state (see
Appendix B Fig. 11). We find that the method reconstructs the free energy
accurately and efficiently in the metastable region provides good estimates of
basin depth and barrier height. Issues of poor sampling affect the calculation
close to the absorbing boundary and the deep minimum of the globally stable
state. We discuss and demonstrate multiple methods to improve our estimates.

This method is applied to unconstrained molecular dynamics trajectories of
supercooled liquid silicon, for which free energy calculations recently
performed using umbrella sampling Monte Carlo exhibited two metastable liquid
states. Using the approach described here, we are able to reproduce these free
energy estimates, confirming the applicability of such an approach in a more
realistic and challenging context.

Quantitative differences in the results for silicon between the first passage
time reconstruction and the method described here, as well as differences at
the lowest temperature studied here remain to be understood better. Moreover,
more exact treatments of the rate in the presence of perturbations away from
equilibrium remain an important open issue to address. The connection to
driven systems is also of importanceCrooks (1998); Kuznets-Speck and Limmer
(2021); Das, Parmar, and Sastry (2022), though the role of the free energy is
less clear for a system driven from equilibrium with an external field.

The framework of the population flux correlation function in defining
phenomenological rates is a significant milestone in the study of rare events
and in the subsequent development of free energy calculation methodsChandler
(1978). Here, by considering the issues arising from a low barrier and thus a
high flux, one can better address the relevant physics in regimes where a
number of approximations common to the high barrier regime do not apply.
Future work that addresses the limitations identified here and strengthens the
connection to driven systems are important avenues to explore.

Figure 7: The free energy as a function of largest crystalline cluster size
($n_{max}$) obtained using unconstrained MD simulations and reconstructed with
either the mean first passage time (curves labelled MFPT) or Eq. 18 used for a
single order parameter, $n_{max}$, (curves labelled 2D) for $T=975~{}K$,
$P=0.75~{}GPa$ (panel a), $T=985~{}K$, $P=0.75~{}GPa$ (panel b) and
$T=995~{}K$, $P=0.75~{}GPa$ (panel c). In order to avoid artefacts due to the
use of the largest cluster size as the order parameter, the curves are shifted
to match the full cluster size distribution for $n_{max}~{}\leq~{}2$ as
described in Ref. Goswami  _et al._ , 2021. The negative logarithm of the
steady state full cluster size distribution is shown for reference in each
case.

Figure 8: The free energy as a function of density alone, from the region of
order parameter space where the largest cluster size is less than the critical
cluster size, reconstructed using either umbrella sampling Monte Carlo
simulations (labelled USMC), reprinted with permission from Goswami and
Sastry, PNAS Nexus, $1$, $4$, (2022). Copyright 2022 Author(s), licensed under
a Creative Commons Attribution (CC BY) license. or Eq. 19 (labelled 2D) for
$T=975~{}K$, $P=0.75~{}GPa$ (panel a), $T=985~{}K$, $P=0.75~{}GPa$ (panel b)
and $T=995~{}K$, $P=0.75~{}GPa$ (panel c). The free energy as a function of
density is obtained from the full two-order parameter distribution using Eq.
22. Figure 9: A comparison of the the free energy surface, $\beta\Delta
G(x,y)$ using Eq. 19 (panel ( a) and error $|\beta\Delta G(x,y)-\beta V(x,y)|$
from the same procedure in units of $k_{B}T^{-1}$ (panel (b). The
corresponding free energy surface and errors using Eq. 23 with an interface
placed at $x=0.5$ are shown in panels (c) and (d) respectively. One observes
an improvement in estimates beyond the barrier due to the improved estimates
of the backward flux. Data is masked for $9x,y)$ values where errors exceeding
the scale shown. Panels (e) and (f) show results using the procedure of
combining estimates from multiple sections with unique reference states $A$.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Potential details

The values of the means and standard deviations for the Gaussian components of
the potential in Fig. 1 are shown in Table 1.

$V_{i}$ | $x_{i}$ | $y_{i}$ | $\sigma x_{i}$ | $\sigma y_{i}$ | $C_{i}$
---|---|---|---|---|---
$V_{1}$ | 0.05 | -0.4 | 0.2 | 0.2 | -0.2
$V_{2}$ | 0.01 | 0.4 | 0.2 | 0.2 | -0.2
$V_{3}$ | 1.0 | -0.4 | 0.1 | 0.1 | -2.0
$V_{4}$ | 0.7 | 0.4 | 0.08 | 0.08 | C
Table 1: The table of factors used to specify the potential energy surface
defined in Eq. 8. The value of $C$ is modulated to produce barriers along $x$
ranging from $0.08$ for $C=1$ to $0.16$ for $C=4$.

The harmonic potential to ensure sampling remains within $-[0.5,0.5]$ is
specified as $V_{c}(y)=k_{c}(|y|-0.45)^{2}$ with a $k_{c}$ value of $20$.

##  Appendix B Tests on an alternate potential - 3 metastable basins or two
saddles

A potential with 3 meta-stable basins, an effective barrier height along $x$
of $3~{}k_{B}T$ and one basin, C, not part of the reaction path (assuming
injection at A) to the globally stable D (see Fig. 10 (panel (a)). Panels (b)
and (c) show the reconstruction using Eq. 18 and Eq. 19 respectively. Errors
are shown in Fig. 10 (d) to show the degree of accuracy. Fig. 10 shows the
reconstruction of the free energy surface and the contracted free energy along
$x$ (panel (e)) and the comparison for a slice along $y$ (panel (f)).

Figure 10: Potential surface with three metastable basins and one globally
stable state (panel (a)). Reconstruction of the two order parameter free
energy using Eq. 18 (panel (b)). Reconstruction of the two order parameter
free energy using Eq. 19 (panel (c)). The relative error, $|\beta\Delta
G(x,y)-\beta V(x,y)|$, units of $k_{B}T$, from the reconstruction using Eq. 19
(panel (d)). Reconstructions of the one order parameter free energy along $x$
(panel (e)) and for a slice along $y$ (panel (f)).

A potential with two paths separating metastable basins from the globally
stable basin. The heights of the saddles along the two paths are unequal to
introduce an asymmetry (see Fig. 11). Fig. 11 shows the reconstruction of the
free energy surface and the contracted free energy along $x$. Errors are shown
in Fig. 11 to show the degree of accuracy.

Figure 11: A potential surface with 2 saddles separating metastable states,
$A,B,C$ from global minimum, $D$ (panel (a)). Reconstruction of the two order-
parameter free energy using Eq. 18 (panel (b)). Reconstruction of the two
order parameter free energy using Eq. 19 (panel (c)). The relative error,
$|\beta\Delta G(x,y)-\beta V(x,y)|$, units of $k_{B}T$ from a reconstruction
using Eq. 19 (panel (d)). Comparison of the free energy for a slice along $y$
for a fixed $x$ (panel (e)).

The results in this section show that the method to reconstruct free energies
is robust to free energy landscapes with a variety of features.

## References

  * Stanley (2013) H. E. Stanley,  _Liquid Polymorphism_ , Vol. 152 (Wiley Online Library, 2013).
  * Palmer  _et al._ (2014) J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos, and P. G. Debenedetti, Nature 510, 385 (2014).
  * Debenedetti, Sciortino, and Zerze (2020) P. G. Debenedetti, F. Sciortino, and G. H. Zerze, Science 369, 289 (2020).
  * Chen, Lascaris, and Palmer (2017) R. Chen, E. Lascaris, and J. C. Palmer, The Journal of chemical physics 146, 234503 (2017).
  * Vasisht and Sastry (2013) V. V. Vasisht and S. Sastry, Liquid Polymorphism 152, 463 (2013), arXiv:2108.13713.
  * Goswami and Sastry (2022) Y. Goswami and S. Sastry, PNAS Nexus 1 (2022), 10.1093/pnasnexus/pgac204.
  * Smallenburg, Filion, and Sciortino (2014) F. Smallenburg, L. Filion, and F. Sciortino, Nature physics 10, 653 (2014).
  * Ricci and Debenedetti (2017) F. Ricci and P. G. Debenedetti, Journal of Chemical Sciences 129, 801 (2017).
  * Torrie and Valleau (1977) G. M. Torrie and J. P. Valleau, Journal of Computational Physics 23, 187 (1977).
  * (10) C. H. Bennett, “Molecular dynamics and transition state theory: The simulation of infrequent events,” in  _Algorithms for Chemical Computations_, Chap. 4, pp. 63–97, https://pubs.acs.org/doi/pdf/10.1021/bk-1977-0046.ch004 .
  * Van Erp, Moroni, and Bolhuis (2003) T. S. Van Erp, D. Moroni, and P. G. Bolhuis, The Journal of chemical physics 118, 7762 (2003).
  * Allen, Frenkel, and ten Wolde (2006) R. J. Allen, D. Frenkel, and P. R. ten Wolde, The Journal of chemical physics 124, 024102 (2006).
  * Valsson, Tiwary, and Parrinello (2016) O. Valsson, P. Tiwary, and M. Parrinello, Annual review of physical chemistry 67, 159 (2016).
  * Peters  _et al._ (2013) B. Peters, P. G. Bolhuis, R. G. Mullen, and J.-E. Shea, The Journal of chemical physics 138, 054106 (2013).
  * Kramers (1940) H. A. Kramers, Physica 7, 284 (1940).
  * Wedekind, Strey, and Reguera (2007) J. Wedekind, R. Strey, and D. Reguera, The Journal of chemical physics 126, 134103 (2007).
  * Wedekind and Reguera (2008) J. Wedekind and D. Reguera, The Journal of Physical Chemistry B 112, 11060 (2008).
  * Wedekind  _et al._ (2009) J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, The Journal of chemical physics 131, 114506 (2009).
  * Lundrigan and Saika-Voivod (2009) S. E. Lundrigan and I. Saika-Voivod, The Journal of Chemical Physics 131, 104503 (2009).
  * Thapar and Escobedo (2015) V. Thapar and F. A. Escobedo, The Journal of chemical physics 143, 244113 (2015).
  * Goswami  _et al._ (2021) Y. Goswami, V. V. Vasisht, D. Frenkel, P. G. Debenedetti, and S. Sastry, The Journal of Chemical Physics 155, 194502 (2021).
  * Crooks (1998) G. E. Crooks, Journal of Statistical Physics 90, 1481 (1998).
  * Banik, Chaudhuri, and Ray (2000) S. K. Banik, J. R. Chaudhuri, and D. S. Ray, The Journal of Chemical Physics 112, 8330 (2000).
  * Hummer and Szabo (2001) G. Hummer and A. Szabo, Proceedings of the National Academy of Sciences 98, 3658 (2001).
  * Poole  _et al._ (1992) P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, Nature 360, 324 (1992).
  * DeFever and Sarupria (2019) R. S. DeFever and S. Sarupria, The Journal of chemical physics 150, 024103 (2019).
  * Verma and Nair (2022) R. Verma and N. N. Nair, The Journal of Physical Chemistry C 126, 19169 (2022).
  * Chandler (1978) D. Chandler, The Journal of Chemical Physics 68, 2959 (1978).
  * Bolhuis  _et al._ (2002) P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annual review of physical chemistry 53, 291 (2002).
  * Van Erp and Bolhuis (2005) T. S. Van Erp and P. G. Bolhuis, Journal of computational Physics 205, 157 (2005).
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber, Physical review B 31, 5262 (1985).
  * Sastry and Angell (2003) S. Sastry and C. A. Angell, Nature materials 2, 739 (2003).
  * Ganesh and Widom (2009) P. Ganesh and M. Widom, Physical review letters 102, 075701 (2009).
  * Beye  _et al._ (2010) M. Beye, F. Sorgenfrei, W. F. Schlotter, W. Wurth, and A. Föhlisch, Proceedings of the National Academy of Sciences 107, 16772 (2010).
  * Vasisht, Saw, and Sastry (2011) V. V. Vasisht, S. Saw, and S. Sastry, Nature Physics 7, 549 (2011).
  * Neophytou, Chakrabarti, and Sciortino (2022) A. Neophytou, D. Chakrabarti, and F. Sciortino, Nature Physics 18, 1248 (2022).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Romano, Sanz, and Sciortino (2011) F. Romano, E. Sanz, and F. Sciortino, The Journal of chemical physics 134, 174502 (2011).
  * Wedekind  _et al._ (2015) J. Wedekind, L. Xu, S. V. Buldyrev, H. E. Stanley, D. Reguera, and G. Franzese, Scientific reports 5, 11260 (2015).
  * Dellago  _et al._ (1998) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, The Journal of chemical physics 108, 1964 (1998).
  * Allen, Valeriani, and Ten Wolde (2009) R. J. Allen, C. Valeriani, and P. R. Ten Wolde, Journal of physics: Condensed matter 21, 463102 (2009).
  * Vanden-Eijnden  _et al._ (2010) E. Vanden-Eijnden _et al._ , Annual review of physical chemistry 61, 391 (2010).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Van Kampen (1992) N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Müller, Talkner, and Reimann (1997) R. Müller, P. Talkner, and P. Reimann, Physica A: Statistical Mechanics and its Applications 247, 338 (1997).
  * Huang and Chen (2021) F. Huang and H. Chen, Physical Review E 103, 062132 (2021).
  * Qin, Dellago, and Kozeschnik (2019) L. Qin, C. Dellago, and E. Kozeschnik, The Journal of chemical physics 150, 094114 (2019).
  * Valeriani  _et al._ (2007) C. Valeriani, R. J. Allen, M. J. Morelli, D. Frenkel, and P. Rein ten Wolde, The Journal of chemical physics 127, 114109 (2007).
  * Kuznets-Speck and Limmer (2021) B. Kuznets-Speck and D. T. Limmer, Proceedings of the National Academy of Sciences 118, e2020863118 (2021).
  * Evans, Majumdar, and Schehr (2020) M. R. Evans, S. N. Majumdar, and G. Schehr, Journal of Physics A: Mathematical and Theoretical 53, 193001 (2020).
  * Chen, Li, and Huang (2022) H. Chen, G. Li, and F. Huang, Journal of Physics A: Mathematical and Theoretical 55, 384005 (2022).
  * Das, Parmar, and Sastry (2022) P. Das, A. D. Parmar, and S. Sastry, The Journal of Chemical Physics 157, 044501 (2022).
  * Filion  _et al._ (2010) L. Filion, M. Hermes, R. Ni, and M. Dijkstra, The Journal of chemical physics 133, 244115 (2010).
  * Haji-Akbari (2018) A. Haji-Akbari, The Journal of chemical physics 149, 072303 (2018).
  * Anderson (1995) J. B. Anderson, Advances in Chemical Physics 91, 381 (1995).
  * Peters (2017) B. Peters,  _Reaction rate theory and rare events_ (Elsevier, 2017).
  * Yamada and Kawasaki (1967) T. Yamada and K. Kawasaki, Progress of Theoretical Physics 38, 1031 (1967).