Targeted free energy estimation via learned mappings

I Introduction

Free energy estimation is of central importance in the natural sciences. Accurate estimation of free energies, however, is challenging, as many systems are out of reach of experimental methods and analytic theory. Computer-based estimation has thus emerged as a valuable alternative. Successful application areas of in-silico free energy estimation span industry and scientific research, including drug discovery Shirts, Mobley, and Brown (2010), condensed matter physics Auer and Frenkel (2001), materials science Damasceno, Engel, and Glotzer (2012), structural biology Curk et al. (2018), and the effects of mutagenesis Hauser et al. (2018). Because of its importance and wide ranging applications, computer-based free energy estimation has been an active field of research for decades Chipot and Pohorille (2007).

At the core of many state-of-the-art estimators Shirts and Chodera (2008) lies the free energy perturbation (FEP) identity introduced by Zwanzig Zwanzig (1954) in 1954:

Here $\Delta F=F_{B}-F_{A}$ is the Helmholtz free energy difference between two thermodynamic states $A$ and $B$, each connected to a thermal reservoir at inverse temperature $\beta$. We denote $\mathbf{x}$ as a point in the system’s configuration space, and define

While Eq. (1) is exact, the convergence of this estimator for a finite number of samples strongly depends on the degree to which $A$ and $B$ overlap in configuration space Pohorille, Jarzynski, and Chipot (2010). Indeed, the dominant contributions to the above expectation will come from samples of $A$ that are typical under $B$, and such contributions become increasingly rare with decreasing overlap Jarzynski (2006).

There exist multiple strategies for mitigating the overlap requirement. Arguably the most common strategy is a multi-staged approach, also known as stratification, in which a sequence of intermediate thermodynamic states is defined between $A$ and $B$ (Fig. 1a). Here the increased convergence is facilitated by demanding that neighboring pairs of states be chosen to contain sufficient overlap. The quantity of interest, $\Delta F$, is then recovered as a sum over the pairwise differences $\Delta F_{i,i+1}$ Chipot and Pohorille (2007). The Multistage Bennett Acceptance Ratio Shirts and Chodera (2008) (MBAR) estimator is a prominent example of an estimator that follows this strategy. However, multi-staged approaches require samples from multiple states as well as a suitable order parameter to define intermediate stages. Furthermore, it is unclear a priori how best to discretize the order parameter or how many stages to use.

An alternative, elegant strategy to increasing overlap is by incorporating configuration space maps. Jarzynski developed Targeted Free Energy Perturbation Jarzynski (2002) (TFEP), a generalization of FEP whereby an invertible mapping defined on configuration space transports points sampled from $A$ to a new distribution, $A^{\prime}$ (see Fig. 1b). Jarzynski showed that a generalized FEP identity can be applied to this process, from which the free energy difference can be recovered. Importantly, if the mapping is chosen wisely, an effective overlap can be increased, leading to quicker convergence of the TFEP estimator. Hahn and Then extended TFEP to the bidirectional setting Hahn and Then (2009), whereby the mapping and its inverse are applied to samples from $A$ and $B$, respectively. Lower-error free energy estimates can then be obtained via the statistically-optimal BAR estimator Bennett (1976).

Whether in the unidirectional or bidirectional case, the main challenge for targeted approaches is crafting a mapping that is capable of increasing overlap. Unfortunately for most real-world problems the physical intuition needed to develop such a technique is simply lacking. Modern-day machine learning (ML) techniques, however, seem perfectly suited for this task.

Since the introduction of TFEP in 2002, research in ML has made remarkable progress in fields of image classification Krizhevsky, Sutskever, and Hinton (2012); He et al. (2016), playing video Mnih et al. (2015) and board games Silver et al. (2016, 2018), and generative modelling of images Brock, Donahue, and Simonyan (2019); Karras, Laine, and Aila (2019). ML has also enabled advances in the natural sciences, including state of the art protein structure prediction Senior et al. (2020), neural-network based molecular force fields Morawietz et al. (2016); Zhang et al. (2018), generative modelling of lattice field theories Albergo, Kanwar, and Shanahan (2019), new paradigms for sampling equilibrium distributions of molecules Noé et al. (2019), and variational free energy estimates Wu, Wang, and Zhang (2019); Li and Wang (2018).

In this work, we turn targeted free energy estimation into a machine learning problem. In lieu of a hand-crafted mapping, we represent our mapping by a deep neural network whose parameters are optimized so as to maximize overlap. Once trained, the free energy can then be computed by evaluating the targeted estimator with our learned mapping. Below we will consider both unidirectional and bidirectional settings, and will refer to them as Learned Free Energy Perturbation (LFEP) and Learned Bennett Acceptance Ratio (LBAR), respectively.

A key contribution of this work is the development of a mapping that respects the underlying symmetries of our system of study. In particular, our neural network is equivariant to permutation of identical particles and respects periodic boundary conditions by construction. These are particularly important considerations when modelling atomic systems, as they often obey such symmetries.

The rest of our manuscript is structured as follows. In Sec. II below we summarize the previously-developed targeted free energy estimators that we will be making use of. This is followed by development of suitable training objectives to maximize overlap (Sec. III). We then demonstrate our method by applying it to a solvation system. In Sec. IV, we describe the experimental setup and discuss inherent symmetries that are exploited to devise a model with the correct inductive biases (Sec. V). Finally, we present experimental results in Sec. VI and discuss our findings in Sec. VII.

II Theoretical background

In the following we will refer to $A$ and $B$ as the thermodynamic states, defined by equilibrium densities $\rho_{A}(\mathbf{x})=e^{-\beta U_{A}(\mathbf{x})}/Z_{A}$ and $\rho_{B}(\mathbf{x})=e^{-\beta U_{B}(\mathbf{x})}/Z_{B}$, where $Z_{A}$ and $Z_{B}$ are the normalization constants (partition functions). In the targeted scheme, configurations $\mathbf{x}$ are drawn from $A$ and mapped to new configurations $\mathbf{y}=M(\mathbf{x})$ via an invertible, user-specified mapping $M$. The set of mapped configurations can be thought of as samples from a new state $A^{\prime}$:

Similarly, we also consider the reverse case where configurations are drawn from $B$ and mapped to $B^{\prime}$ via the inverse

We refer to this pair of prescriptions as the “forward” and “reverse” process, respectively. For each process we denote generalized energy differences as

where $J_{M}$ and $J_{M^{-1}}=J_{M}^{-1}$ are the Jacobian determinants associated with the mappings. As originally shown by Jarzynski Jarzynski (2002), an identity exists which relates $\Delta F$ to an ensemble of realizations of $\Phi_{F}$:

Eq. (7) can be regarded as a generalization of FEP, as it holds for any invertible $M$, and reduces to Eq. (1) if $M$ is the identity. An analogous equation holds for the reverse process. Derivations of Eq. (7) can be found in Refs. Jarzynski, 2002; Hahn and Then, 2009 or Appendix A.

Hahn and Then extended the above result to the bidirectional case Hahn and Then (2009), showing a fluctuation theorem (FT) exists between the forward and reverse processes:

can be thought of as generalized work distributions associated with the mapping processes and $\delta$ is the Dirac delta function. With these bidirectional estimates, Bennett’s Acceptance Ratio (BAR) method Bennett (1976) can be employed as an alternative estimator of ${\Delta F}$ Hahn and Then (2009). BAR estimation of $\Delta F$ can be formulated as a self-consistent iteration of the equation

where $f(x)=1/(1+e^{x})$ is the Fermi function. For simplicity in Eq. (11) we restrict ourselves to the case where the number of samples in the forward and reverse directions are equal, but more general formulations exist. BAR has a statistical advantage over FEP as it has been shown to be the minimum variance free energy estimator for any asymptotically-unbiased method Shirts et al. (2003). Because of this property, BAR is generally the method of choice when samples from both $A$ and $B$ are available.

In summary, our method proceeds in two stages by first computing optimized work values using Eqs. (5)–(6) and then estimating $\Delta F$. For the latter we can employ the generalized FEP estimator (7) in the unidirectional setting (LFEP), or solve the BAR equations (11) in the bidirectional setting (LBAR). This highlights an important difference between LBAR and other maximum likelihood free energy estimators which assume the work values to be fixed Bennett (1976); Maragakis, Spichty, and Karplus (2006); Shirts and Chodera (2008). Instead, LBAR learns to optimize the work values and subsequently combines them optimally to predict $\Delta F$.

Crucially, the targeted estimators above hold for every invertible mapping. That is, given an infinite number of samples, any invertible choice of $M$ will produce a consistent estimate of $\Delta F$. Of course, the finite-sample convergence properties are of more practical importance and will strongly depend on the choice of $M$.

III Training objective

In a distributional sense, the forward and reverse processes act to transform $\rho_{A}$ and $\rho_{B}$ into $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$, as depicted in Fig 1. In what follows we will refer to the distribution of mapped configurations as the “images” (i.e. $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$), and the distributions we want them mapped towards as the “targets” (i.e. $\rho_{A}$ ($\rho_{B}$) for the forward (reverse) process). Due to the deterministic mapping, the bases and images are related by the change of variable formula,

The crucial consideration for convergence of our estimators (Eq. (7) and Eq. (11)) is the overlap between the image and target distributions Jarzynski (2002). Indeed, in the limit that images and targets coincide, Jarzynski showed Jarzynski (2002) that $p_{F}(\phi)\rightarrow\delta\mathopen{}\left(\phi-\Delta F\right)\mathclose{}$. This implies that the convergence of Eq. (7) is immediate (i.e. only one sample is needed). In this limit, it is also the case that $p_{R}(\phi)\rightarrow\delta\mathopen{}\left(\phi+\Delta F\right)\mathclose{}$ implying that the expectation values on either side of Eq. (11) converge immediately. This overlap argument is further reinforced in the Appendix A, where we show that the TFEP estimator can be interpreted as a FEP estimator between $A^{\prime}$ and $B$.

We now turn our attention to the construction of a loss function that accurately judges the quality of $M$. Guided by the considerations of overlap, we consider the Kullback–Leibler (KL) divergence between the image and target. For the forward process we have:

In the above derivation, we invoked a change of variable formula in going from the first to third lines, and used the identity $-\beta\Delta F=\log{Z_{B}}-\log{Z_{A}}$ to get to the last. An analogous equation can be derived for the reverse process yielding

While from Eqs. (14–15) it is clear that the KL cannot be accurately estimated unless ${\Delta}F$ is known, in terms of optimizing, ${\Delta}F$ and $\beta$ can be disregarded as they are constants.

Below we consider two separate training regimes for our model. In the unidirectional case, the model was trained only on the forward process, with a loss function

In the bidirectional case, the model was trained using both forward and reverse processes with the loss

When samples from both states are available, bidirectional training is preferable. Unlike the unidirectional loss, $\mathcal{L}_{\text{LBAR}}$ explicitly encourages both $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$ to be mass-covering (Minka, 2005), which is important for good performance of importance-sampling estimators Neal (2005).

IV Experimental setup

To test our method, we consider a system similar to the one used by Jarzynski Jarzynski (2002) consisting of a repulsive solute immersed in a bath of $N=125$ identical solvent particles. The task is to calculate the free energy change associated with growing the solute radius from $R_{A}$ to radius $R_{B}$ (see Fig. 2). In contrast to a hard-sphere solute as used in Ref. Jarzynski, 2002, we modeled our solute as a soft sphere. This is because any finite particle overlap would lead to infinite forces, which our training method cannot handle.

Intuitively, an effective mapping should push solvent particles away from the center to avoid high-energy steric clashes with the expanding solute. Jarzynski followed this intuition, defining a mapping that uniformly compresses the solvent particles amidst an expanding repulsive solute. Although this mapping gave a significant convergence gains when applied to a hard solute Jarzynski (2002), it is not directly applicable to soft solutes. This is because the phase space compression results in a transformed density whose support is not equal to that of the target density, violating the assumption of invertibility.

Below we demonstrate that we no longer need to rely on physical intuition to hand-craft a tractable mapping; this process can be fully automated using the general framework proposed in this work. In order to learn an effective mapping, however, it is crucial that the model be compatible with the inherent symmetries of the underlying physical system.

IV.2 Training data

IV.3 Symmetries

The system under consideration exhibits properties that are widely encountered in the atomistic simulation community: periodic boundary conditions (PBCs) and permutation invariance. PBCs are usually employed to reduce finite-size effects. Permutation invariance arises as a consequence of the energy being invariant to particle permutations—a condition that is satisfied by the solvent particles as they are all identical.

We design the mapping $M$ to respect PBCs and permutation invariance by construction. This means that the state $A^{\prime}$ obtained by transforming $A$ via $M$ is guaranteed to have the 3D torus geometry stipulated by PBCs, and to be symmetric with respect to any permutation of solvent particles. The next section discusses in detail how these symmetries are implemented in the model architecture.

Our model architecture does not obey the octahedral symmetries, in the sense that $A^{\prime}$ is not guaranteed to be symmetric with respect to the $48$ permutations and/or reflections of the three coordinate axes. However, since the size of this symmetry group is small, we account for the octahedral symmetries via training-data augmentation instead. That is, during training we transform every training data point (system configuration) by a random element of the octahedral group. As training progresses, all $48$ transformations of each data point are likely to be seen by the model, thus the model is trained to learn these symmetries from data. In comparison, exhausting the $N!=125!$ permutation symmetries by training-data augmentation would be infeasible in any reasonable training time.

V Model

We implement the mapping $M$ using a deep neural network, parameterized by a set of learnable parameters $\theta$. In designing the architecture of the network, we take into account the following considerations.

The mapping $M$ must be bijective, and the inverse mapping $M^{-1}$ should be efficient to compute, for any setting of the parameters $\theta$.

The Jacobian determinant $J_{M}$ should be efficient to compute for any setting of $\theta$.

The network should be flexible enough to represent complex mappings.

The transformed distributions $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$ should respect the boundary conditions and symmetries of the physical system.

The first three requirements are satisfied by a class of deep neural networks known as normalizing flows Papamakarios et al. (2019), which are invertible networks with efficient Jacobian determinants. Since bijectivity is a closed property under function composition, multiple normalizing flows (or “layers”) can be composed into a deeper flow, yielding a model with increased flexibility. We implement $M$ as a normalizing flow composed of $K$ invertible layers, that is,

For the mapping $M$ to be bijective, a sufficient condition is that $G(\cdot;\psi):[-L,L]\rightarrow[-L,L]$ be strictly increasing for any setting of $\psi$. In that case, the inverse $M_{k}^{-1}$ is obtained by simply replacing $G$ with $G^{-1}$ in Eq. (20), and the Jacobian determinants can be computed efficiently as follows:

Finally, the inverse and Jacobian determinant of the composite mapping $M$ can be computed by

To ensure that the transformed distributions $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$ obey the required boundary conditions, the implementation of $G$ must reflect the fact that $r_{i}^{\nu}=-L$ and $r_{i}^{\nu}=L$ are identified as the same point. For this to be the case, a sufficient set of conditions is the following:

for any setting of the parameters $\psi$. To satisfy the above conditions, we implement $G$ using circular splines, which were recently proposed by Rezende et al. Jimenez Rezende et al. (2020) and are based on the rational-quadratic spline flows of Durkan et al. Durkan et al. (2019). Our implementation of the circular splines is detailed in Appendix C.

Finally, to ensure that the transformed distributions $\rho_{A^{\prime}}$ and $\rho_{B^{\prime}}$ are invariant to particle permutations, it is necessary that the Jacobian determinant $J_{M}$ also be invariant to particle permutations. In our architecture, this can be achieved by taking $C$ to be equivariant to particle permutations. Specifically, let $\sigma$ be a permutation of the set $\{1,\ldots,N\}$. We say that $C$ is equivariant with respect to $\sigma$ if

that is, if permuting the particles has the effect of permuting the parameter outputs $\mathopen{}\left(\psi_{1}^{\nu},\ldots,\psi_{N}^{\nu}\right)\mathclose{}$ in exactly the same way. From Eq. (22), we can easily see that the above property implies that $\sigma$ leaves $J_{M_{k}}$, and hence $J_{M}$, invariant, because we sum over all particles and the sum is permutation invariant. Previous studies have made similar observations Bender et al. (2020); Köhler, Klein, and Noé (2019). Our implementation of $C$ is based on the architecture proposed by Vaswani et al. Vaswani et al. (2017), often referred to as the transformer, which we use in a permutation-equivariant configuration. The implementation details of our transformer architecture are in Appendix C.

VI Results

In this section, we evaluate the performance of our method for the solvation system illustrated in Fig. 2. We focus on the bidirectional BAR and LBAR estimators in the main text, due to their advantages over unidirectional approaches as discussed in Sec. II. We refer to Appendix D for a discussion of the unidirectional counterparts.

To capture statistical variation, our training and analysis procedure was performed $10$ times, each using independent training and evaluation datasets. Our loss profiles and free energy estimates below report averages as well as statistical variation across these runs.

We first report training results in Fig. 4, where the full-batch loss is plotted as a function of the number of training steps. We observe a pattern commonly encountered in ML: after an initial decrease of both training and test loss, the latter develops a minimum. At around the minimum, the model stops generalizing and starts to overfit to the training data. We therefore employ a technique called early stopping Caruana, Lawrence, and Giles (2001) and use the model parameters corresponding to the minimum test loss for all further evaluations. It is worth emphasizing, however, that the precise location of the minimum does not have an appreciable effect on the quality of the free energy estimates reported for the bidirectional estimator (results not presented). We also note the small variation among the independent runs, suggesting that there is no significant dependence of the performance on a particular dataset.

We now turn to the statistical convergence of the free energy estimates. In Fig. 5b we plot a running average of the estimate as a function of the number of evaluated samples per stage. The solid lines report averages of the estimate over the independent runs, and shaded regions represent one standard deviation of the runs. We first validate the correctness of our method against a converged MBAR estimator. Here MBAR employed $15$ stages and thus $7.5$ times more samples in total. From the figure we see that the variation of our estimate overlaps nicely with the MBAR error estimates. We next compare the efficiency of our method against the baseline BAR estimator, where we see in Fig. 5b clear variance reduction of LBAR across a wide range of evaluated sample sizes. The full-batch LBAR standard deviation we observe is approximately $19\%$ of that reported by BAR. Moreover, training and evaluation of the model occurred on the same dataset, demonstrating that an effective mapping can be learned in a data-efficient manner. This is an important practical consideration but not at all obvious a priori. We could, in principle, even combine samples from the training and test datasets for estimation of $\Delta F$ but have used the test set only to detect overfitting.

VII Discussion

In this work, we turn TFEP into a machine learning problem by combining it with state-of-the-art ML techniques. TFEP previously required hand-crafting a tractable mapping on configuration-space, a significant challenge for many realistic systems. We proposed to represent the mapping by a suitable neural network, and identified training objectives for unidirectional and bidirectional cases. We then tested their performance on a prototype solvation system – the growth of a soft sphere in a fluid of solvent particles in periodic boundary conditions. While this system is relatively simplistic from a physical standpoint, it poses a significant challenge for ML models due to the system’s underlying permutational symmetry and periodic boundary conditions. Our experimental results indicate that both LFEP and LBAR estimators can lead to a significant variance reduction compared to their respective baselines and therefore clearly highlight the potential of this approach. Interesting directions for future work include a systematic analysis of the error of the learned estimators and a detailed comparison with MBAR. We believe it is possible that optimal estimation strategies on complex systems will contain a combination of staging and mapping.

Improving TFEP via learned mappings relates to the general idea of improving importance sampling by learning the proposal distribution, which has been explored substantially in machine learning and statistics. For instance, recent works in machine learning have proposed training a flexible deep-learning model of the proposal distribution to improve importance sampling Müller et al. (2019) or more sophisticated variants such as bridge sampling Papamakarios and Murray (2015) and sequential Monte Carlo Gu, Ghahramani, and Turner (2015); Paige and Wood (2016); Le, Baydin, and Wood (2017). In turn, these approaches can be traced back to methods for adaptive importance sampling Cappé et al. (2008) and adaptive sequential Monte Carlo Cornebise, Moulines, and Olsson (2008) in statistics. One recent instance of these approaches that relates closely to our work is Neural Importance Sampling Müller et al. (2019), which uses expressive normalizing flows to learn good proposal distributions for importance sampling. Many of the above works have noted that the choice of loss function is important, with the forward KL divergence and chi-squared divergence being standard choices. These observations are in line with our observations of the differences between the unidirectional and bidirectional training losses.

We note that our learned free energy estimators and equivariant, periodic model architecture can be combined with the work on Boltzmann Generators Noé et al. (2019), which uses flow-based models in combination with statistical re-weighting to sample from desired Boltzmann distributions. In particular, it would be interesting to see how our targeted estimators compare to the ones considered in Ref. Noé et al., 2019. Furthermore, we note that neural network based free energy estimation is an active field of research. For example, two recent studies have independently proposed targeted unidirectional Nicoli et al. (2020) and bidirectional Ding and Zhang (2019) estimators similar to the ones suggested here. Both studies employ autoregressive networks to compute free energy estimates of a lattice spin model with discrete states. In addition, Ref. Ding and Zhang, 2019 also estimates free energies of a small protein in the gas phase (no periodicity) using normalizing flows. As noted in that study, however, their model lacks permutation equivariance which is likely a drawback when applying it to the type of solvation system we consider here. We also believe that permutation equivariant, periodicity-respecting networks, such as the one considered here, will be key to scaling up the approach to system sizes commonly used in atomistic simulations.

Finally, our results demonstrate that we can estimate free energy differences between two states directly and data efficiently, i.e. using fewer MD samples for training and evaluation than the base estimator would require to converge. Other studies, for example Refs. Noé et al., 2019; Ding and Zhang, 2019, follow a different strategy and estimate free energy differences by learning two separate mappings that share a common reference state. We believe that our direct approach may be preferable in cases where one state is a small perturbation of the other. It will be interesting to see how these different approaches compare with respect to data efficiency and variance reduction, and under which circumstances one is preferable to the other.

Appendix A Alternate derivation of TFEP and interpretation

In this section we derive and interpret the unidirectional TFEP estimator as a multi-staged FEP estimator. This interpretation allows us to reason about TFEP using intuition from FEP.

Given explicit densities for $A$, $A^{\prime}$ and $B$, we can formally decompose $\Delta F$ into a sum over two terms,

which can each be computed separately using the FEP estimator Eq. (1). We next define the energy of $A^{\prime}$ as

such that $\rho_{A^{\prime}}\propto e^{-\beta U_{A^{\prime}}}$, and where we have used Eq. (12) in going from the first to second line. Conveniently, one of these stages comes for free, as $\Delta F_{AA^{\prime}}=0$:

In going from the first to second line we have used Eq. (30). Combining Eqs. (29–31), as well as the FEP estimator Eq. (1), we arrive at our final result:

is the energy difference between $A^{\prime}$ and $B$. Although Eq. (32) is an explicit estimate between $A^{\prime}$ and $B$, it is just a reformulation of the TFEP estimator [Eq. (7) above] as can be seen by the equivalence of $\Delta U^{\prime}$ and $\Phi_{F}$ [compare Eqs. (5) and (33)]. This interpretation of TFEP allows us to apply the intuition on convergence we have built for FEP. Specifically, we can accelerate convergence if $A^{\prime}$ shares large overlap with $B$.

Appendix B System

To generate the training data, we performed MD simulations of the system illustrated in Fig. 2 using the simulation package LAMMPS Plimpton (1995). The system is similar to the one studied in Ref. Jarzynski, 2002 but we replaced hard solute-solvent interactions by a Weeks–Chandler–Andersen Weeks and Chandler (1971) (WCA) potential. Below, we represent our energy, length, and mass units in terms of the LJ well depth $\epsilon$, the LJ diameter $\sigma$, and the solvent particle mass $m$ Frenkel and Smit (2002). From this our unit of time is defined as $\tau=\sigma\sqrt{m/\epsilon}$. Quantities expressed in these reduced units are denoted with an asterisk. We used a cubic simulation box with edge length $2L^{*}=6.29$, employed cutoff radii of $L^{*}$ and $\sqrt{2}R_{\alpha}$ for LJ and WCA interactions, where $\alpha\in\{A,B\}$ labels the state. The solute radii were taken to be $R_{A}^{*}=2.5974$ and $R_{B}^{*}=2.8444$. Both LJ and WCA potentials shared the same value for $\epsilon$, and we set $\sigma_{\text{WCA}}=R_{\alpha}$.

To simulate a specific state, we first assigned random positions to all solvent particles. The solute was placed at the origin of the box and kept stationary throughout the simulation. After performing energy minimization, the system was equilibrated in the canonical ensemble for a period of $5\times 10^{4}\tau$. The equations of motion were integrated using the velocity Verlet algorithm Swope et al. (1982) with a timestep $0.002\tau$. We employed a Langevin thermostat with a relaxation time of $0.5\tau$ to keep the system at a temperature of $T^{*}=3.2155$. To prevent drift of the center of mass motion, we set the total random force to zero. During the $5\times 10^{5}\tau$ long production run, configurations were sampled every $10$ reduced time units, yielding a total of $5\text{\times}{10}^{4}$ samples. For each state, $20$ such simulations were generated starting from random initializations and random number seeds. The resulting $1\text{\times}{10}^{6}$ samples were then partitioned as described in the main text.

We followed the same protocol to generate samples for MBAR. In addition to the two states corresponding to $R_{A}$ and $R_{B}$, we considered $13$ intermediate states with a constant radial increment $\Delta R_{i,i+1}$. Using two different random seeds, we obtained $1\text{\times}{10}^{5}$ samples for each state and evaluated all $15$ energy functions on each sample. MBAR estimates of $\Delta F$ were then computed from this combined energy matrix using the package pymbar Shirts and Chodera (2008).

Appendix C Model implementation details

We implement $G$ using circular splines Jimenez Rezende et al. (2020) so that $G$ satisfies the boundary conditions in Eqs. (25–26). Briefly, $G$ is a piecewise function that consists of $S$ segments. The parameters $\psi$ are a set of $S+1$ triplets $(x_{s},y_{s},d_{s})$, where $[x_{s-1},x_{s}]$ defines the domain of segment $s$, $[y_{s-1},y_{s}]$ defines its image, and $d_{s-1},d_{s}$ define its slopes at the endpoints (all slopes are required to be positive). Each segment is a strictly increasing rational-quadratic function, constructed using the interpolation method of Gregory and Delbourgo Gregory and Delbourgo (1982). The conditions in Eqs. (25–26) are satisfied by setting $x_{0}=y_{0}=-L$, $x_{S}=y_{S}=L$ and $d_{0}=d_{S}>0$. We can increase the flexibility of $G$ by increasing the number of segments $S$. With a large enough $S$, circular splines can approximate arbitrarily well any strictly increasing function from $[-L,L]$ to itself that satisfies Eqs. (25–26).

All models were trained using the Adam optimizer Kingma and Ba (2015). A summary of the hyperparameters is provided in Table 1.

Appendix D Results for LFEP

In this section, we discuss the free energy estimates obtained with LFEP, using unidirectional training with the $\mathcal{L}_{\text{LFEP}}$ loss [Eq. (16)] as a proxy for $D_{\text{KL}}$ [Eq. (14)]. To estimate the KL, we replaced $\Delta F$ in Eq. (14) with the LFEP estimate of that quantity. Figure 6a shows the evolution of both quantities during training. The results feature an interesting behaviour in the initial training regime. While test and training loss are still decreasing, the KL already exhibits a pronounced minimum due to a drift of the $\Delta F$ estimate towards lower values. This is further illustrated in Fig. 6b, which compares the convergence of $\Delta F$ for three different mappings corresponding to specific training steps (colored symbols in Fig. 6a). We see that the variance is reduced significantly in all three cases. However, only the first mapping (training step $2\text{\times}{10}^{3}$) agrees with the MBAR baseline, while we can already observe a small bias for the second mapping (training step $1.1\text{\times}{10}^{4}$) that becomes even more pronounced as training progresses (training step $2\text{\times}{10}^{4}$). We note that this behaviour is consistent with other experiments that we performed in the unidirectional setting (data not shown).

This is problematic in that the minimum of the KL would be a natural point to stop training but does not yield the lowest bias. This is also in contrast to our findings for the bidirectional case, where the quality of the mapping was best in the vicinity of the minimum in the test loss. One possible explanation for these observations is the well-known zero-forcing property Minka (2005) of $\mathcal{L}_{\text{LFEP}}$. That is, $\mathcal{L}_{\text{LFEP}}$ does not encourage the transformed distribution to be mass-covering, which is known to negatively impact the performance of importance-sampling estimators Neal (2005).

References