# A Boltzmann generator for the isobaric-isothermal ensemble

Steyn van Leeuwen,1 Alberto Pérez de Alba Ortíz 1,2 and Marjolein Dijkstra
1111Email<EMAIL_ADDRESS>
(1Soft Condensed Matter, Debye Institute for Nanomaterials Science, Utrecht
University, Princetonplein 1, 3584 CC Utrecht, The Netherlands
2Computational Soft Matter, van ’t Hoff Institute for Molecular Sciences and
Informatics Institute, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands.
)

###### Abstract

Boltzmann generators (BGs) are now recognized as forefront generative models
for sampling equilibrium states of many-body systems in the canonical
ensemble, as well as for calculating the corresponding Helmholtz free energy.
Furthermore, BGs can potentially provide a notable improvement in efficiency
compared to conventional techniques such as molecular dynamics (MD) and Monte
Carlo (MC) methods. By sampling from a clustered latent space, BGs can
circumvent free-energy barriers and overcome the rare-event problem. However,
one major limitation of BGs is their inability to sample across phase
transitions between ordered phases. This is due to the fact that new phases
may not be commensurate with the box dimensions, which remain fixed in the
canonical ensemble. In this work, we present a novel BG model for the
isothermal-isobaric ($NPT$) ensemble, which can successfully overcome this
limitation. This unsupervised machine-learning model can sample equilibrium
states at various pressures, as well as pressure-driven phase transitions. We
demonstrate that the samples generated by this model are in good agreement
with those obtained through MD simulations of two model systems. Additionally,
we derive an estimate of the Gibbs free energy using samples generated by the
$NPT$ BG.

##  1 INTRODUCTION

Boltzmann generators (BGs), [20] first introduced in 2019, offer a promising
new approach for efficiently sampling the equilibrium states of many-body
systems. Traditionally, sampling the Boltzmann distribution[4] of atomic or
colloidal systems relies on a step-wise propagation, either via molecular
dynamics (MD)[2] or Monte Carlo (MC)[18] simulations. Although these methods
are from the 1950s, across decades they have been immensely boosted by
advances in computational power and by enhanced sampling algorithms, [26, 16,
10] which enable to explore and render free-energy landscapes of many-body
systems. These accelerated computer experiments are nowadays routinely used to
extract valuable insight about complex phenomena occurring in a variety of
systems, ranging from molecules, to polymers, to colloids. The ability to
sample the Boltzmann distribution is essential for understanding such many-
body systems, because it determines their equilibrium properties. Therefore
the simulation community is constantly seeking for more efficient ways to
perform this task in order to tackle increasingly complex systems. Although
coordinate transformations were proposed as a solution for barrier-less
sampling, [30] the case-by-case nature of such transformations impeded
progress. Fortunately, recent advances in machine learning have provided a
solution. BGs can provide one-shot samples from the Boltzmann distribution
without the need for time-stepping of MC moves, once a suitable transformation
has been machine-learned. As the simulation community starts to invest more in
improving BGs, in the same way that MD and MC were boosted decades ago, we
present here an extension of BGs, which to our knowledge is the first BG for
the isobaric-isothermal ensemble.

Boltzmann Generators (BGs) were proposed by Noé et al. [20] in 2019 for
efficiently sampling the equilibrium states of many-body condensed-matter
systems in one shot. BGs leverage normalizing flows (NFs), which are trainable
and invertible transformations, to transform a simple prior distribution
$\mu_{Z}$ (e.g. a Gaussian distribution) in latent space $Z$ into the
Boltzmann distribution $\mu_{X}$ in configuration space $X$ through a
coordinate transformation $\mathbf{x}=F_{zx}(\mathbf{z})$. Once the NF is
trained, one can obtain configurations $\mathbf{x}$ by sampling from the prior
distribution ($\mathbf{z}\sim\mu_{Z}$) and transforming those samples to
configurations $\mathbf{x}=F_{zx}(\mathbf{z})$. In this way, BGs bypass the
step-wise nature of previous sampling methods and efficiently obtain
configurations in one shot. As a result, BGs eliminate the need for an
iterative process of updating a single configuration and climbing over free-
energy barriers. This innovative approach to sampling provides a powerful and
scalable tool for exploring the equilibrium properties of complex systems in
various fields.

In recent years, BGs have been improved to handle the large symmetry groups
that are associated with many-body systems. [1, 29, 28, 15, 24, 27] However,
it is important to note that all these BGs sample many-body systems in the
canonical ensemble, i.e. the system is at a constant number of particles $N$,
temperature $T$, and volume $V$. This constant volume constraint can cause
problems when sampling across phase transitions, since a crystal may not be
commensurate with the fixed box dimensions. The isobaric-isothermal ensemble
circumvents this problem since the pressure $P$ is held constant instead of
the volume, allowing the volume to fluctuate such that different crystals may
be accommodated. Additionally, this ensemble enables the calculation of the
Gibbs free energy as opposed to the Helmholtz free energy in the canonical
ensemble.

In this work, we explore how BGs can be extended for sampling many-body
systems in the isobaric-isothermal ensemble, i.e. at constant number of
particles $N$, pressure $P$, and temperature $T$. More precisely, we
generalize an existing method by Ahmad and Cai[1] from the canonical ensemble
to the isobaric-isothermal ensemble. Using two case studies, namely a Lennard-
Jones (LJ) and a Hemmer-Stell-like (HSL) system, we compare the performance of
the new $NPT$ BG to MD simulations, and observe good agreement between the two
methods. Additionally, we also derive an expression to extract the Gibbs free
energy from the $NPT$ BG.

black

##  2 METHODS

###  2.1 Isobaric-isothermal ensemble

In the isobaric-isothermal ensemble, we consider configurations with a fixed
number of particles $N$, fixed pressure $P$, and fixed temperature $T$.[8]
Since $P$ is fixed, the volume $V$ of the system is allowed to fluctuate.
Therefore, a configuration in the $NPT$ ensemble is not only described by
particle positions $\mathbf{x}$, but also by the box dimensions
$\mathbf{L}=(L_{x},L_{y})$ in a two-dimensional, or
$\mathbf{L}=(L_{x},L_{y},L_{z})$ in a three-dimensional system, which are
allowed to fluctuate. More precisely, a configuration
($\mathbf{L},\mathbf{s})$ can be defined in 2D by the box dimensions
$\mathbf{L}=(L_{x},L_{y})$ and by the scaled positions of all its particles,
i.e. $\mathbf{s}=(\mathbf{s}_{1},...,\mathbf{s}_{N})$, where
$\mathbf{s}_{i}=(x_{i}/L_{x},y_{i}/L_{y})$ denotes the coordinates of particle
$i$ scaled by the box dimensions. The probability of a configuration
$(\mathbf{L},\mathbf{s})$ is given by

$\mu_{X}(\mathbf{L},\mathbf{s})\propto V^{N}e^{-\beta PV}e^{-\beta
U(\mathbf{s},\mathbf{L})},$ (1)

which we refer to as the isobaric-isothermal probability distribution, and
where $V$ denotes the volume of the system, i.e. $V=L_{x}L_{y}$ in 2D.
Furthermore, $U(\mathbf{x})$ is the potential energy of configuration
$\mathbf{x}$ and $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann constant.

The isobaric-isothermal partition function $\Delta$ is given by

$\begin{split}\Delta(N,P,T)&=\frac{C}{\Lambda^{DN}N!}\int\,dV\ V^{N}e^{-\beta
PV}\int\,d\mathbf{s}^{N}\ e^{-\beta U(\mathbf{s},\mathbf{L})},\\\
&=\frac{C}{\Lambda^{DN}N!}\Delta_{X},\end{split}$ (2)

where $C$ is some constant with units of inverse volume, usually set to
$C=\beta P$, and $\Delta_{X}$ is the configurational partition function. The
Gibbs free energy $G$ can be expressed in terms of $\Delta$ as

$\begin{split}\beta G=-\log(\Delta)=-\log\left(\frac{\beta
P}{\Lambda^{DN}N!}\right)+\beta G_{X},\end{split}$ (3)

where $\beta G_{X}=-\log(\Delta_{X})$ is the configurational contribution to
the Gibbs free energy. The difficulty in computing $\Delta$ and $G$ lies in
computing $\Delta_{X}$ and $G_{X}$, i.e. the configurational contributions. In
this work, we focus on estimating $\Delta_{X}$ and $G_{X}$ by generating
samples according to the isobaric-isothermal probability distribution

$\begin{split}\mu_{X}(\mathbf{L},\mathbf{s})&=\dfrac{1}{\Delta_{X}}e^{-\beta
U(\mathbf{L},\mathbf{s})-\beta PV+N\log(V)}.\end{split}$ (4)

Furthermore, we are interested in computing ensemble averages of observables
over the isobaric-isothermal probability distribution. The average of an
observable $O$ over $\mu_{X}$ is given by

$\displaystyle\left<O(\mathbf{L},\mathbf{s})\right>_{(\mathbf{L},\mathbf{s})\sim\mu_{X}}=\int\,d\mathbf{L}\,d\mathbf{s}\
\mu_{X}(\mathbf{L},\mathbf{s})O(\mathbf{L},\mathbf{s}).$ (5)

###  2.2 Boltzmann generators

Figure 1: The main idea behind a NF is the transformation of probability
distributions. Samples in latent space $Z$ are generated by sampling from a
simple distribution $\mu_{Z}(\mathbf{z})$ (e.g. a Gaussian distribution). Then
a trainable and invertible transformation $F_{zx}$ maps the samples in $Z$ to
samples in configuration space $X$. Samples in $X$ follow a distribution
$q_{X}(\mathbf{x})$ which is similar to the target distribution $\mu_{X}$
(e.g. a Boltzmann distribution). Furthermore, samples in $X$ (obtained from MD
or MC simulation) can also be mapped to $Z$ by $F_{xz}$. This results in a
distribution $q_{Z}$ which is similar to $\mu_{Z}$. Both directions can be
used to train the invertible transformation $F_{xz}=(F_{zx})^{-1}$.

The primary goal of a BG is to sample configurations of a particle system
without the step-wise propagation used in MD and MC methods. The aim is to
directly draw samples of particle configurations, denoted by $\mathbf{x}$,
from a configuration space $X$ that follows the Boltzmann distribution
$\mu_{X}(\mathbf{x})$. However, this is a challenging task, as the
distribution is not known a priori. Despite the unknown distribution, it is
possible to calculate the likelihood of a specific sample $\mathbf{x}$, and
use that to our advantage. Specifically, one can sample from a simple, known
distribution in a latent space $Z$ and then learn to transform a sample
$\bf{z}$ into a sample $\mathbf{x}$, which follows the Boltzmann distribution.
This transformation, known as an NF, denoted as
$\mathbf{x}=F_{zx}(\mathbf{z})$, can be machine-learned.

Let us now describe in more detail how NFs generate samples. First, samples
are drawn from a simple prior probability distribution $\mu_{Z}(\mathbf{z})$,
which could be a Gaussian distribution, in latent space $Z$, see Fig. 1. These
samples in $Z$ are then transformed into samples in configuration space $X$
through a machine-learned invertible and differentiable transformation
$\mathbf{x}=F_{zx}(\mathbf{z})$. If $F_{zx}$ has been learned correctly, the
resulting samples in $X$ follow a distribution $q_{X}(\mathbf{x})$ that
closely resembles the Boltzmann distribution $\mu_{X}(\mathbf{x})$. Moreover,
the transformation $F_{zx}$ can also be inverted as $F_{xz}=(F_{zx})^{-1}$,
allowing training in both directions. We note that the subscripts ${xz}$ or
${zx}$ indicate the direction of the transformation, i.e. $X$-to-$Z$ or
$Z$-to-$X$, respectively black

BGs transform samples in latent space to configuration space by means of an
NF. An NF is a series of learnable, invertible and differentiable
transformations

$F_{zx}=f_{n}\circ...\circ f_{2}\circ f_{1}.$ (6)

The transformation $F_{zx}$ generates the new distribution by translating and
locally compressing or stretching space. One of the advantages of BGs is that
the probability $q_{X}(\mathbf{x})$ of a generated sample $\mathbf{x}$ can be
directly evaluated. More precisely, $q_{X}(\mathbf{x})$ can be related to the
prior probability distribution $\mu_{Z}(\mathbf{z})$ through a change of
variables

$q_{X}(\mathbf{x})=\mu_{Z}(F_{xz}(\mathbf{x}))R_{xz}(\mathbf{x}),$ (7)

where
$R_{xz}(\mathbf{x})=|\det(J_{xz}(\mathbf{x}))|=|\det\nabla_{\mathbf{x}}F_{xz}(\mathbf{x},\boldsymbol{\theta})|$.
Note that $q_{X}(\mathbf{x})$ depends on the likelihood of the prior
distribution $\mu_{Z}(F_{xz}(\mathbf{x}))$, but also on how much the space is
locally compressed or stretched, which is quantified by $R_{xz}(\mathbf{x})$.
Furthermore, note that it is necessary for $F_{xz}$ to fulfill certain
conditions: $F_{xz}$ must be invertible, because this ensures $\mathbf{x}$ to
be mapped to a unique $\mathbf{z}=F_{xz}(\mathbf{x})$, and $F_{xz}$ must be
differentiable so that $R_{xz}$ is well-defined. Finally, we remark that the
advantage of such an exact-likelihood generative model is that $q_{X}$ can be
used to reweight samples.

Additionally, Fig. 1 shows that $F_{xz}=(F_{zx})^{-1}$ can also be used to map
configurations from $X$ to $Z$. This allows us to train the NF in two
ways.[20] The first way is by sampling a batch $B$ from $\mu_{Z}$, mapping $B$
to $X$ via $F_{zx}$ and improving $F_{zx}$ by minimizing a loss function that
is based on the difference between $q_{X}(\mathbf{x})$ and
$\mu_{X}(\mathbf{x})$ for all $\mathbf{x}\in F_{zx}(B)$. The second way is by
using an MC simulation to sample a batch $B$ from $\mu_{X}$, mapping $B$ to
$Z$ via $F_{xz}$ and updating $F_{xz}$ through a loss function which minimizes
the difference between $q_{Z}(\mathbf{z})$ and $\mu_{Z}(\mathbf{z})$ for all
$\mathbf{z}\in F_{xz}(B)$. More details on the loss functions used for
training the NFs are presented in Section 2.4.

Ideally, NFs should be easy to compute and invert, and the determinant of
their Jacobian should be easy to calculate. There are many implementations
that meet these requirements, [21, 13] such as coupling flows, [1, 20] auto-
regressive flows, residual flows, [29, 28] and infinitesimal flows. [24, 15]
In this work, we use the real-valued non-volume preserving (RealNVP) coupling
layers [6] to implement $f_{i}$ in Eq. 6. In a RealNVP coupling layer, the
input variable $\mathbf{z}$ is split into two channels:
$\mathbf{z}_{A}=\mathbf{z}_{1:d}$, $\mathbf{z}_{B}=\mathbf{z}_{d+1:D}$ (See
Fig. 2). Here, $D$ is the dimensionality of the system and $1\leq d<D$. The
coupling layer does not transform $\mathbf{z}_{A}$ and just copies it to
obtain the output $\mathbf{x}_{A}$. However, $\mathbf{z}_{B}$ is transformed
using an affine transformation that uses $\mathbf{z}_{A}$ as input

$f_{zx}(\mathbf{z}_{A},\mathbf{z}_{B})=\begin{cases}\mathbf{x}_{A}&=\mathbf{z}_{A}\\\
\mathbf{x}_{B}&=\exp(-S(\mathbf{z}_{A}))\odot\left(\mathbf{z}_{B}-T(\mathbf{z}_{A})\right),\end{cases}$
(8)

where $S:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D-d}$ and
$T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D-d}$ are machine-learnable functions
(usually deep neural networks), $\odot$ denotes element-wise multiplication
and $\exp$ is applied element-wise to $-S(\mathbf{z}_{A})$. To transform all
coordinates, the roles of $\mathbf{z}_{B}$ and $\mathbf{z}_{A}$ are switched
in the next affine coupling layer so that $\mathbf{z}_{A}$ is transformed and
$\mathbf{z}_{B}$ is copied. This combination of two affine coupling layers is
referred to as a RealNVP block.

Figure 2: The affine coupling layers in the RealNVP blocks. The input variable
$\mathbf{z}$ is split into two channels: $\mathbf{z}_{A}$ and
$\mathbf{z}_{B}$. The first channel is simply copied, while the second channel
is transformed using an affine transformation that uses the first channel as
input.

This scheme has all the properties we desire. It is easy to compute since it
is just a linear transformation. Furthermore, it is easy to invert because
$\mathbf{z}_{A}$ is not transformed, so that
$S(\mathbf{x}_{A})=S(\mathbf{z}_{A})$ and
$T(\mathbf{x}_{A})=T(\mathbf{z}_{A})$ can be readily computed and therefore
$S$ and $T$ do not need to be inverted. Therefore the inverse transformation
is given by

$f_{xz}(\mathbf{x}_{A},\mathbf{x}_{B})=\begin{cases}&\mathbf{z}_{A}=\mathbf{x}_{A},\\\
&\mathbf{z}_{B}=\exp(S(\mathbf{x}_{A}))\odot\mathbf{x}_{B}+T(\mathbf{x}_{A}),\end{cases}$
(9)

and finally, the Jacobian of $f_{zx}$ and $f_{xz}$ can be easily computed. For
$f_{xz}$, the Jacobian is given by

$\begin{split}\nabla_{\mathbf{x}}f_{xz}(\mathbf{x}_{A},\mathbf{x}_{B})&=\begin{pmatrix}\partial_{\mathbf{x}_{A}}\mathbf{z}_{A}&\partial_{\mathbf{x}_{B}}\mathbf{z}_{A}\\\
\partial_{\mathbf{x}_{A}}\mathbf{z}_{B}&\partial_{\mathbf{x}_{B}}\mathbf{z}_{B}\end{pmatrix},\\\
&=\begin{pmatrix}\mathbb{I}&0\\\
\partial_{\mathbf{x}_{A}}\mathbf{z}_{B}&\text{diag}(\exp{S(\mathbf{x}_{A})})\end{pmatrix},\end{split}$
(10)

and its Jacobian can be written as
$|\det\nabla_{\mathbf{x}}f_{xz}|=\prod_{i=1}^{D-d}\exp{S_{i}(\mathbf{x}_{A})}$.
Furthermore, the logarithm of the determinant is given by

$\log|\det\nabla_{\mathbf{x}}f_{xz}|=\sum_{i=1}^{D-d}S_{i}(\mathbf{x}_{A}),$
(11)

which can be computed easily. To summarize, BGs use NFs to generate samples
according to the Boltzmann distribution. NFs are differentiable bijections
that transform a simple distribution $\mu_{Z}$, e.g. a Gaussian, into a
complex distribution $q_{X}$ which approximates some desired distribution
$\mu_{X}$, e.g. the Boltzmann distribution. These NFs can be used to generate
samples according to $\mu_{X}$ by sampling from $\mu_{Z}$ and mapping these
samples to $X$ via $F_{zx}$. Furthermore, we can compute the exact likelihood
$q_{X}$ for all samples generated by the NF. Finally, the NF can be trained in
two ways: by improving $F_{zx}$ and by improving $F_{xz}$.

###  2.3 Isobaric-isothermal Boltzmann generator

In this section, we introduce the isothermal-isobaric Boltzmann generator
($NPT$ BG) for a system in 2D. Inspired by the work of Ahmad and Cai[1] we
generate displacements from a reference lattice, rather than absolute
coordinates. The $NPT$ BG generates configurations $(\mathbf{L},\mathbf{s})$
according to the $NPT$ distribution with isotropic box fluctuations. Here,
$\mathbf{L}=(L_{x},L_{y})$ is the box dimension in the $x$-direction and
$y$-direction, and $\mathbf{s}=(\mathbf{s}_{1},...,\mathbf{s}_{N})$ are scaled
deviations with respect to a scaled lattice
$\mathbf{r}_{0}=(\mathbf{r}_{0}^{1},...,\mathbf{r}_{0}^{N})$ (see Fig. 3). The
absolute coordinates of the system
$\mathbf{R}=(\mathbf{R}_{1},...,\mathbf{R}_{N})=(R_{1x}.R_{1y},...,R_{Nx},R_{Ny})$
are given by

$\begin{split}R_{i\alpha}&=L_{\alpha}((s_{i\alpha}+r_{0}^{i\alpha})\
\mathrm{mod}\ 1),\end{split}$ (12)

where $i$ is the particle number, $\alpha=x,y$ and $\mathrm{mod}\ 1$ applies
the periodic boundary conditions to the scaled coordinate
$s_{i\alpha}+r_{0}^{i\alpha}$.

(a)

(b)

(c)

Figure 3: The $NPT$ BG generates scaled deviations
$\mathbf{s}=(\mathbf{s}_{1},...,\mathbf{s}_{N})$ with respect to a scaled
reference lattice
$\mathbf{r}_{0}=(\mathbf{r}_{0}^{1},...,\mathbf{r}_{0}^{N})$. (a) The
reference lattice $\mathbf{r}_{0}$. Note that the lattice spacing of
$\mathbf{r}_{0}$ is smaller in the $x$-direction than in the $y$-direction.
This is because there are more particle layers in the $x$-direction. (b) The
deviations $\mathbf{s}$ with respect to $\mathbf{r}_{0}$. Note that
$\mathbf{s}_{3}$ places particle $3$ outside the box. The particle can be
mapped back into the box by applying $\mathrm{mod}\ 1$ to
$\mathbf{r}_{0}^{3}+\mathbf{s}_{3}$. (c) absolute particle coordinates and box
sizes, obtained by scaling by $\mathbf{L}=(L_{x},L_{y})$.

Any generated vector $\mathbf{v}\in X$ first lists $\mathbf{L}$ and then all
coordinates $\mathbf{s}$, such that
$\mathbf{v}=(\mathbf{L},\mathbf{s}_{1},...,\mathbf{s}_{N})$. Thus, Eq. 7
becomes

$\begin{split}q_{X}(\mathbf{L},\mathbf{s})=\mu_{Z}\left(F_{xz}(\mathbf{L},\mathbf{s})\right)R_{xz}(\mathbf{L},\mathbf{s}),\end{split}$
(13)

where
$R_{xz}(\mathbf{L},\mathbf{s})=|\det(\nabla_{(\mathbf{L},\mathbf{s})}F_{xz}(\mathbf{L},\mathbf{s}))|$
and $\mu_{Z}(\mathbf{z})=\exp(-|\mathbf{z}|^{2}/2\sigma^{2})/\Delta_{Z}$. Note
that, since the box fluctuations are isotropic, box proportions are fixed
(i.e. $L_{y}=\lambda L_{x}$ for some constant $\lambda$). Therefore, it is
sufficient to generate only $L_{x}$ and not the entire vector
$\mathbf{L}=(L_{x},L_{y})$. Furthermore, the area $V$ can be computed via
$V=\lambda L_{x}^{2}$ in 2D. In practice, within the NF, we add one more
dimension to the latent space $Z$ and configuration space $X$, so that
$\dim(Z)=\dim(X)=ND+1$.

It should be noted that the $NPT$ BG can generate unbounded displacements, as
well as configurations with $L_{x}<0$. To prevent this, we propose an NF
$F_{zx}:Z\rightarrow X$ with three transformations, given by

$F_{zx}=F_{\text{lin}}\circ F_{\text{scale}}\circ F_{\text{RealNVP}},$ (14)

where $F_{\text{RealNVP}}:\mathbb{R}^{ND+1}\rightarrow\mathbb{R}^{ND+1}$ is a
standard NF with the RealNVP architecture,
$F_{\text{scale}}:\mathbb{R}^{ND+1}\rightarrow\mathbb{R}^{ND+1}$ is a scaling
layer and $F_{\text{lin}}$ is a linear translation. More precisely,
$F_{\text{scale}}$ is given by

$F_{\text{scale}}(\mathbf{y})=C_{\text{scale}}\mathbf{y},$ (15)

with $C_{\text{scale}}$ a machine learnable parameter, and $F_{\text{lin}}$
reads

$F_{\text{lin}}(\mathbf{y}^{\prime})=\mathbf{y}^{\prime}+(C_{\text{lin}},0,...,0),$
(16)

where $C_{\text{lin}}$ is also a machine learnable parameter and $\mathbf{y}$
and $\mathbf{y}^{\prime}$ denote the intermediate coordinate transformations
from $\mathbf{z}$ to $\mathbf{x}$.

$F_{\text{scale}}$ must be initialized such that the displacements are not
larger than half the lattice spacing, i.e.
$|\mathbf{s}_{i}|<a_{\text{scaled}}/2$ for all $i$. To this end, we first set
$C_{\text{scale}}=1$ and $C_{\text{lin}}=0$ and generate a batch $B$ of
configurations $\mathbf{v}=(\mathbf{L},\mathbf{s})$. Then we compute the
maximum over all coordinates $v_{i}$ of $\mathbf{v}$ over all $\mathbf{v}\in
B$ (i.e., $v_{\text{max}}=\max_{\mathbf{v}\in B}\max_{1\leq i\leq
ND+1}|v_{i}|$). We then set $C_{\text{scale}}=1/(\gamma*v_{max})$, where
$\gamma$ is a constant that should be chosen such that $1/\gamma<1/2\sqrt{N}$.
This ensures that initial deviations will have a maximum around $1/\gamma$ and
that deviations larger than $1/2\sqrt{N}$ are very unlikely. Since
$a_{\text{scaled}}\sim 1/\sqrt{N}$, this condition also ensures that
$|\mathbf{s}_{i}|<a_{\text{scaled}}/2$ for all $i$.

On the other hand, $F_{\text{lin}}$ must be initialised such that the BG does
not generate configurations with $L_{x}<0$. We perform MD simulations and
initialise $C_{\text{lin}}$ to approximate the $L_{x}$ that was observed in
these simulations. Therefore, $F_{\text{lin}}$ ensures that the initial
distribution of $L_{x}$ is centered around $C_{\text{lin}}$. Furthermore,
$F_{\text{scale}}$ ensures that the maximum deviation from $C_{\text{lin}}$ is
$1/\gamma<1/\sqrt{N}$. Therefore, the minimum $L_{x}$ that can be generated is
around $C_{\text{lin}}-1/\gamma$, which is much larger than $0$ for the
systems we consider.

To ensure that we sample only one permutation and remove translational
invariance, we define an augmented system with a center-of-mass restraint

$U_{\text{aug}}(\mathbf{L},\mathbf{s})=U(\mathbf{L},\mathbf{s})+\dfrac{1}{2}k_{\text{CoM}}|\mathbf{s}_{\text{CoM}}|^{2},$
(17)

where $\mathbf{s}$ are deviations with respect to $\mathbf{r}_{0}$,
$k_{\text{CoM}}$ is the strength of the harmonic center-of-mass potential and
$\mathbf{s}_{\text{CoM}}$ is the displacement vector of the center of mass,
which is given by

$\begin{split}\mathbf{s}_{\text{CoM}}=\dfrac{1}{N}\sum_{i=1}^{N}\mathbf{s}_{i}.\end{split}$
(18)

Note that it is not necessary to apply periodic boundary conditions.

To select an appropriate value for $k_{\text{CoM}}$, we examine the average
energy resulting from the harmonic potential in $2D$

$U_{\text{avg}}=N\dfrac{D}{2}k_{B}T=Nk_{B}T,$ (19)

and require that the energy of the harmonic potential at the largest possible
deviation is much larger than the average energy of the system

$\begin{split}U_{\text{CoM}}(1/2\sqrt{N})=\dfrac{1}{2}k_{\text{CoM}}(1/2\sqrt{N})^{2}&\gg
U_{\text{avg}}=Nk_{b}T,\\\ k_{\text{CoM}}&\gg 8N^{2}k_{B}T.\end{split}$ (20)

To summarize, our $NPT$ BG generates scaled deviations $\mathbf{s}$ w.r.t. a
fixed reference lattice $\mathbf{r}_{0}$. The initial deviations can be kept
small by using a scaling layer in the NF. This ensures that only one
permutation is sampled. Furthermore, the $NPT$ BG generates the box size
$\mathbf{L}$. Absolute particle coordinates can be obtained by multiplying the
scaled particle coordinates ($\mathbf{r}_{0}+\mathbf{s}$) by
$\mathbf{L}=(L_{x},L_{y})$.

Our implementation of BGs is partially based on a repository by Hsu and Fobe,
[11] who implemented a double-well potential that we extended to other
systems. Additionally, we used utilities from the seminal work on BGs by Noé
et al. [20] and their repository on BGs and NFs. [14] black

###  2.4 Loss functions in the isothermal-isobaric ensemble

In this section, we derive the Maximum likelihood (ML) and the Kullback-
Leibler (KL) loss for the $NPT$ ensemble. The training of our invertible
transformations, $F_{zx}$ and $F_{xz}$, proceeds by minimizing specific loss
functions. In particular, when learning $F_{xz}$, we use the maximum
likelihood (ML) loss, which measures the difference between the approximate
probability distribution $q_{Z}$ and the ground-truth distribution $\mu_{Z}$.
Conversely, when learning $F_{zx}$, we use the Kullback-Leibler (KL) loss,
which measures the difference between the approximate probability distribution
$q_{X}$ and the ground-truth distribution $\mu_{X}$. While we discuss the ML
loss for completeness, in this work we focus solely on training based on the
KL loss to use the BG in an unsupervised way.

####  2.4.1 Maximum likelihood loss

To derive the ML loss function, we take the logarithm of Eq. 13 and write

$\begin{split}\log
q_{X}(\mathbf{L},\mathbf{s})&=\log(\mu_{Z}\left(F_{xz}(\mathbf{L},\mathbf{s})\right))+\log
R_{xz}(\mathbf{L},\mathbf{s}),\\\ &=-\log
Z_{Z}-\frac{|F_{xz}(\mathbf{L},\mathbf{s})|^{2}}{2\sigma^{2}}+\log
R_{xz}(\mathbf{L},\mathbf{s}),\end{split}$ (21)

where
$R_{xz}(\mathbf{L},\mathbf{s})=|\det(\nabla_{(\mathbf{L},\mathbf{s})}F_{xz}(\mathbf{L},\mathbf{s}))|$
and $\mu_{Z}(\mathbf{z})=\exp(-|\mathbf{z}|^{2}/2\sigma^{2})/Z_{Z}$, with
$\sigma$ the standard deviation of the Gaussian prior and $Z_{Z}$ the
normalization factor of the Gaussian prior. We can then write the KL-
divergence between $\mu_{X}$ and $q_{X}$ as

$\displaystyle\text{KL}\left(\mu_{X}\mid q_{X}\right)$ (22)
$\displaystyle=\int\,d\mathbf{L}\,d\mathbf{s}\
\mu_{X}(\mathbf{L},\mathbf{s})\left[\log\mu_{X}(\mathbf{L},\mathbf{s})-\log
q_{X}(\mathbf{L},\mathbf{s})\right],$
$\displaystyle=\mathbb{E}_{(\mathbf{L},\mathbf{s})\sim\mu_{X}}\left[\frac{\left|F_{xz}(\mathbf{L},\mathbf{s})\right|^{2}}{2\sigma^{2}}-\log
R_{xz}(\mathbf{L},\mathbf{s})\right]+\text{c},$

where, in the last line, we write $\int\,d\mathbf{L}\,d\mathbf{s}\
\mu_{X}(\mathbf{L},\mathbf{s})\left[\log(\mu_{X}(\mathbf{L},\mathbf{s}))+\log(Z_{Z})\right]$
as a constant term $c$ because they do not depend on machine learnable
parameters. The ML loss is then defined as

$L_{\text{ML}}=\mathbb{E}_{(\mathbf{L},\mathbf{s})\sim\mu_{X}}\left[\frac{\left|F_{xz}(\mathbf{L},\mathbf{s},\boldsymbol{\theta})\right|^{2}}{2\sigma^{2}}-\log
R_{xz}(\mathbf{L},\mathbf{s},\boldsymbol{\theta})\right],$ (23)

where $\boldsymbol{\theta}$ are the machine learnable parameters in $F_{xz}$.
We can use samples from MC simulations or MD simulations to approximate
$\mu_{X}$.

####  2.4.2 Kullback-Leibler loss in the $NPT$ ensemble

We now derive the KL loss. We start by noting that
$q_{Z}(\mathbf{z})=\mu_{X}(F_{zx}(\mathbf{z}))R_{zx}(\mathbf{z})$ and
substituting $\mu_{X}$ from Eq. 4, such that

$\displaystyle\log q_{Z}(\mathbf{z})=$
$\displaystyle\log(\mu_{X}(F_{zx}(\mathbf{z})))+\log(R_{zx}(\mathbf{z})),$
(24) $\displaystyle=$ $\displaystyle-\log(\Delta_{X})-\beta
U(\mathbf{L},\mathbf{s})-\beta PV$
$\displaystyle+N\log(V)+\log(R_{zx}(\mathbf{z})).$

We write the KL-divergence between $p=\mu_{Z}$ and $q=q_{Z}$ as

$\begin{split}\text{KL}(\mu_{Z}|q_{Z})\\\ =&\int\,d\mathbf{z}\
\mu_{Z}(\mathbf{z})\left[\log(\mu_{Z}(\mathbf{z}))-\log(q_{Z}(\mathbf{z}))\right],\\\
=&\int\,d\mathbf{z}\
\mu_{Z}(\mathbf{z})[\log(\mu_{Z}(\mathbf{z}))+\log(\Delta_{X})\\\ &+\beta
U(\mathbf{L},\mathbf{s})+\beta PV-N\log(V)-\log(R_{zx}(\mathbf{z}))],\\\
=&\mathbb{E}_{\mathbf{z}\sim\mu_{Z}}[\beta U(\mathbf{L},\mathbf{s})+\beta
PV-N\log(V)\\\ &-\log(R_{zx}(\mathbf{z}))]+c,\end{split}$ (25)

where we have written $\int\,d\mathbf{z}\
\mu_{Z}(\mathbf{z})\left[\log(\mu_{Z}(\mathbf{z}))+\log(\Delta_{X})\right]$ as
a constant term $c$ since it does not depend on $\boldsymbol{\theta}$. Note
that $V$ depends on $\mathbf{L}$. The KL loss can then be defined as

$\displaystyle L_{\text{KL}}=$ (26)
$\displaystyle\mathbb{E}_{\mathbf{z}\sim\mu_{Z}}\left[\beta
U(\mathbf{L},\mathbf{s})+\beta
PV-N\log(V)-\log(R_{zx}(\mathbf{z},\boldsymbol{\theta}))\right].$

####  2.4.3 Gibbs free energy

When a physical process takes place under constant temperature and pressure,
it can be described in terms of the Gibbs free energy. The Gibbs free energy
is a crucial thermodynamic quantity for studying phase transitions in many-
body systems since the difference between the Gibbs free energy of two states
and the Gibbs free-energy barrier determine whether a transition can occur
spontaneously. In this section, we derive the configurational contribution to
the Gibbs free energy $G_{X}$ using samples from the probability distribution
$q_{X}$ generated by the $NPT$ BG. This derivation can be divided into three
steps. First, based on the augmented system defined in Eq 17, we express the
free energy $G^{\text{aug}}_{X}$ of the distribution $\mu^{\text{aug}}_{X}$ by
using free-energy perturbation,[31] a technique in which the thermodynamic
properties of one system can be calculated based on a slightly different
system and on the difference between the interparticle potentials of the two
systems. Second, we calculate the Gibbs free-energy difference between the
augmented system, i.e. with the center-of-mass restraint defined in Eq. 17,
and the non-augmented system, i.e. without the center-of-mass restraint.
Finally, we use the previous two steps to compute the free energy $G_{X}$ of
$\mu_{X}$.

Let us first define the Gibbs free energy $G^{\text{aug}}_{X}$ of the
augmented system whose energy is given by Eq. 17. This Gibbs free energy can
be expressed in terms of the partition function $\Delta^{\text{aug}}_{X}$ as

$\beta G^{\text{aug}}_{X}=-\log(\Delta^{\text{aug}}_{X}),$ (27)

and $\Delta^{\text{aug}}_{X}$ can be estimated using samples from $q_{X}$ as
follows

$\displaystyle\Delta^{\text{aug}}_{X}$ (28) $\displaystyle=\int\,d\mathbf{L}\
\int\,d\mathbf{s}\ e^{-\beta U_{\text{aug}}(\mathbf{L},\mathbf{s})-\beta
PV(\mathbf{L})+N\log(V(\mathbf{L}))},$ $\displaystyle=N!\int\,d\mathbf{L}\
\int_{\mathbf{s}^{\dagger}}\,d\mathbf{s}\ e^{-\beta
U_{\text{aug}}(\mathbf{L},\mathbf{s})-\beta
PV(\mathbf{L})+N\log(V(\mathbf{L}))},$ $\displaystyle=N!\
\mathbb{E}_{(\mathbf{L},\mathbf{s})\sim q_{X}}\left[\frac{e^{-\beta
U_{\text{aug}}(\mathbf{L},\mathbf{s})-\beta
PV(\mathbf{L})+N\log(V(\mathbf{L}))}}{q_{X}(\mathbf{L},\mathbf{s})}\right],$

where $\mathbf{s}^{\dagger}$ denotes one permutation of $\mathbf{s}$.

Plugging Eq. 28 into Eq. 27, we find that the Gibbs free energy of the
augmented system can be computed from samples from $q_{X}$

$\displaystyle\beta G^{\text{aug}}_{X}=$
$\displaystyle-\log\left(\mathbb{E}_{(\mathbf{L},\mathbf{s})\sim
q_{X}}\left[e^{-\beta
U_{G}^{\text{aug}}(\mathbf{L},\mathbf{s})}\dfrac{1}{q_{X}(\mathbf{L},\mathbf{s})}\right]\right)$
(29) $\displaystyle-\log(N!).$

The previous equation allows us to compute $G^{\text{aug}}_{X}$ in terms of
samples from $q_{X}$. However, we are interested in the Gibbs free energy of
the non-augmented system $G_{X}$. Therefore, we write the free-energy
difference between the augmented and the non-augmented system as

$\begin{split}\beta\Delta G&=\beta G-\beta
G^{\text{aug}}_{X}=\log\left(\dfrac{\Delta^{\text{aug}}_{X}}{\Delta_{X}}\right)\\\
&=\log\dfrac{\int\,d\mathbf{L}\ \int\,d\mathbf{s}\ e^{-\beta
U_{\text{aug}}(\mathbf{L},\mathbf{s})-\beta PV(\mathbf{L})+N\log
V(\mathbf{L})}}{\int\,d\mathbf{L}\ \int\,d\mathbf{s}\ e^{-\beta
U(\mathbf{L},\mathbf{s})-\beta PV(\mathbf{L})+N\log
V(\mathbf{L})}}.\end{split}$ (30)

Because $U$ is invariant under global translations of the system, it can be
rewritten as

$U(\mathbf{x}_{1},...,\mathbf{x}_{N})=U(\mathbf{0},\mathbf{x}^{\prime}_{2},...,\mathbf{x}^{\prime}_{N})=U(\mathbf{x}^{\prime}_{2},...,\mathbf{x}^{\prime}_{N})$
(31)

where $\mathbf{x}^{\prime}_{i}=\mathbf{x}_{i}-\mathbf{x}_{1}$ for all $i$ and
we remove the first coordinate for the sake of simplicity. Let us now consider
a change of variables from
$(\mathbf{L},\mathbf{s}_{1},\mathbf{s}_{2},...,\mathbf{s}_{N})$ to
$(\mathbf{L},\mathbf{s}_{\text{CoM}},\mathbf{s}_{2}^{\prime},...,\mathbf{s}_{N}^{\prime})$,
where

$\begin{split}\mathbf{s}_{\text{CoM}}&=\dfrac{1}{N}\sum_{i=1}^{N}\mathbf{s}_{i}\\\
\mathbf{s}_{i}^{\prime}&=\mathbf{s}_{i}-\mathbf{s}_{1}.\end{split}$ (32)

Then, we can substitute $U_{\text{aug}}(\mathbf{L},\mathbf{s})$ by
$U(\mathbf{L},\mathbf{s}_{2}^{\prime}...,\mathbf{s}_{N}^{\prime})+U_{\text{CoM}}(\mathbf{s}_{\text{CoM}})$,
and $U(\mathbf{L},\mathbf{s})$ by
$U(\mathbf{L},\mathbf{s}_{2}^{\prime}...,\mathbf{s}_{N}^{\prime})$ in Eq. 30,
and integrate with respect to $\,d\mathbf{s}_{\text{CoM}}\
d\mathbf{s}_{2}^{\prime}...\mathbf{s}_{N}^{\prime}$ instead of $d\mathbf{s}$.
We thus simplify Eq. 30 as

$\begin{split}\beta\Delta
G&=\log\dfrac{\int\,d\mathbf{s}_{\text{CoM}}\exp\left(-\beta
U_{\text{CoM}}(\mathbf{s}_{\text{CoM}})\right)}{\int\,d\mathbf{s}_{\text{CoM}}\
},\end{split}$ (33)

The two integrals over $\mathbf{s}_{\text{CoM}}$ in Eq. 33 are easy to
compute, since one integrand is $1$ and the other integrand is $\exp(-\beta
U_{\text{CoM}}(\mathbf{s}_{\text{CoM}}))=\exp(-\beta
k_{\text{CoM}}\mathbf{s}_{\text{CoM}}^{2}/2)$, which is just a Gaussian.
However, the exact evaluation of these integrals is beyond the scope of this
work.

We now have an estimate of the Gibbs free energy of the augmented system via a
learned free-energy perturbation (LFEP) approach and an analytic expression
for the Gibbs free-energy contribution of the CoM constraint. Therefore, the
Gibbs free energy of the non-augmented system can be computed as

$\begin{split}\beta G_{X}=&\beta G^{\text{aug}}_{X}+\beta\Delta
G_{\text{CoM}}\\\ =&-\log\left(\mathbb{E}_{(\mathbf{L},\mathbf{s})\sim
q_{X}}\left[e^{-\beta(U_{G}^{\text{aug}}(\mathbf{L},\mathbf{s})}\dfrac{1}{q_{X}(\mathbf{L},\mathbf{s})}\right]\right)\\\
&-\log(N!)+\log\dfrac{\int\,d\mathbf{s}_{\text{CoM}}\exp\left(-\beta
U_{\text{CoM}}(\mathbf{s}_{\text{CoM}})\right)}{\int\,d\mathbf{s}_{\text{CoM}}\
1}.\end{split}$ (34)

Most importantly, when calculating Gibbs free-energy differences, the last two
terms in Eq. 34 cancel when the system’s number of particles, reference
lattice and center-of-mass restraint are kept unchanged.

###  2.5 Observables

In this section, we discuss observables that can be computed from
configurations sampled either by a BG or MD. More precisely, we discuss the
radial distribution function and the instantaneous pressure. These observables
will be computed for two case studies in Sections 3.1 and Section 3.2.

####  2.5.1 Radial distribution function

The radial distribution function (RDF) is a measure of the two-body
correlations in liquids and crystals. The RDF $g(r)$ is defined as the ratio
between the number density of particles at a certain distance $r$ from a fixed
particle and the expected number density at distance $r$ for an ideal gas with
the same density. [8] Numerically, we can compute $g(r)$ for one particle
position $\mathbf{r}_{i}$ in configuration
$\mathbf{x}=(\mathbf{r}_{1},...,\mathbf{r}_{N})$. To get more statistics, we
average over all particles ($i=1,...,N$) and over multiple configurations
($\\{\mathbf{x}_{1},...,\mathbf{x}_{B}\\}$). In the $NPT$ ensemble,
configurations are given by $(\mathbf{L},\mathbf{s})$, so the box size and
volume vary with each configuration. Therefore, we compute the RDF as follows

$g(r)=\dfrac{1}{BN}\sum_{j=1}^{B}\sum_{i=1}^{N}\dfrac{\text{count}(r,i,(\mathbf{L},\mathbf{s})_{j})}{\rho^{\prime}((\mathbf{L},\mathbf{s})_{j})V_{\text{shell}(r)}},$
(35)

where $\text{count}(r,i,(\mathbf{L},\mathbf{s})_{j})$ is the number of
particles at distance $[r,r+dr)$ from particle $i$ in configuration
$(\mathbf{L},\mathbf{s})_{j}$, and
$\rho^{\prime}((\mathbf{L},\mathbf{s})_{j})=(N-1)/V(\mathbf{L})$, where
$V(\mathbf{L})=L_{x}^{2}\lambda$ is the volume of the box in configuration
$(\mathbf{L},\mathbf{s})_{j}$ and $V_{\text{shell}}(r)=\pi(r+dr)^{2}-\pi
r^{2}$ is the volume of a circular shell around $\mathbf{r}_{i}$.

####  2.5.2 Pressure

Another quantity we can compute is the instantaneous pressure as defined by
the virial equation. For the $NPT$ ensemble, the macroscopic pressure is
fixed, because the system is coupled to a (fictitious) barostat at pressure
$P$. The instantaneous pressure $P_{\text{instant}}(\mathbf{L},\mathbf{s})$
fluctuates around $P$ and, in MD simulations, the volume fluctuates based on
the pressure difference $P-P_{\text{instant}}(\mathbf{L},\mathbf{s})$ between
the external and internal system. The instantaneous pressure can be computed
via

$P_{\text{instant}}(\mathbf{L},\mathbf{s})=\dfrac{Nk_{B}T}{V(\mathbf{L})}+\dfrac{\sum_{i=1}^{N}\mathbf{F}_{i}\cdot\mathbf{x}_{i}}{DV(\mathbf{L})},$
(36)

where $V(\mathbf{L})=L_{x}^{2}\lambda$ is the volume of the system, $D$ is the
dimensionality of the system and
$\mathbf{F}_{i}=\nabla_{\mathbf{x}_{i}}U(\mathbf{x}_{1},...,\mathbf{x}_{N})$
is the force on particle $i$.

##  3 RESULTS

###  3.1 Lennard-Jones system in the $NPT$ ensemble

Figure 4: The Lennard-Jones (LJ) and Hemmer-Stell-like (HSL) pair potentials.

In our first case study, we consider a two-dimensional system with periodic
boundary conditions in the $NPT$ ensemble. The system consists of $N=36$
Lennard-Jones (LJ) particles at pressure $P\sigma^{2}/\epsilon=5.0$ and
temperature $k_{B}T/\epsilon=0.5$, where the system has been shown to exhibit
a hexagonal crystal phase. [3] The LJ potential is given by

$U_{\text{LJ}}(r_{ij})=\sum^{N}_{j>i}4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],\\\
$ (37)

where $r_{ij}$ is the inter-particle distance between particle $i$ and $j$
computed using the nearest image convention. Furthermore, $\sigma$ is the
particle diameter and $\epsilon$ is the interaction strength (see Fig. 4). We
use $\sigma$ and $\epsilon$ as our units of length and energy, respectively.
The potential is truncated and shifted to zero at a chosen cut-off distance
$r_{\text{cut}}$. The settings for all these parameters are listed in Table 1.

####  3.1.1 Boltzmann generator protocol

We add a harmonic center-of-mass restraint in terms of the deviations
$\mathbf{s}$ to remove the translation invariance. This potential is given by

$U_{\text{CoM}}(\mathbf{s})=\dfrac{1}{2}k_{\text{CoM}}\mathbf{s}_{\text{CoM}}^{2}=\dfrac{1}{2}k_{\text{CoM}}\left(\dfrac{1}{N}\sum_{i=1}^{N}\mathbf{s}_{i}\right)^{2},$
(38)

where $k_{\text{CoM}}\gg 8N^{2}k_{B}T$ in order to only sample one permutation
(see Eq. 20 for a derivation). In this case we have $N^{2}k_{B}T=8\cdot
36^{2}\cdot 0.5\epsilon=5184\epsilon$. Hence, we set
$k_{\text{CoM}}=3*10^{5}\epsilon$.

Once the translation and permutation are handled, we must also handle
numerical instabilities during training. The divergence in
$\nabla_{r}U_{\text{LJ}}(r)$ for $r\downarrow 0$ can lead to instabilities in
training the BG. To overcome this problem, we regularize the potential.[29] We
linearize the LJ pair interaction below a distance $r_{\text{lin}}$, such that

$\displaystyle U_{\text{LJ}}^{\text{reg}}(r)$ (39)
$\displaystyle=\begin{cases}U_{\text{LJ}}(r_{\text{lin}})+\left(\dfrac{\partial
U_{\text{LJ}}}{\partial
r}\right)\Bigr{|}_{r=r_{\text{lin}}}(r-r_{\text{lin}})&r<r_{\text{lin}}\\\
U_{\text{LJ}}(r)&r\geq r_{\text{lin}}\end{cases}.$

On the one hand $r_{\text{lin}}$ should be chosen far enough away from $r=0$
in order to ensure that $\partial_{r}U_{\text{LJ}}^{\text{reg}}(r)$ is small
enough. On the other hand, $r_{\text{lin}}$ should be chosen small enough to
avoid influencing the generated distribution. More precisely, $r_{\text{lin}}$
should be chosen so that $U_{\text{LJ}}(r_{\text{lin}})\gg k_{B}T$, because
deviations of the order of $k_{B}T$ are exponentially unlikely. In practice,
we set $r_{\text{lin}}=0.8\sigma$ so that

$\begin{split}U_{\text{LJ}}(r)\Bigr{|}_{r=r_{\text{lin}}=0.8\sigma}&\approx
49\epsilon;\\\ \left(\dfrac{\partial U_{\text{LJ}}}{\partial
r}\right)\Bigr{|}_{r=r_{\text{lin}}=0.8\sigma}&\approx-759\dfrac{\epsilon}{\sigma}.\end{split}$
(40)

We train a $NPT$ BG on KL-loss only to use the BG in an unsupervised way. We
use the RealNVP architecture described in Section 2.2. We train over several
epochs, where each epoch can have different hyperparameters (e.g. batch size
or learning rate). For the specifics about the architecture and
hyperparameters used in each epoch, as well as the regularization of the LJ
interaction potential, see Section A.1.

black

####  3.1.2 Molecular dynamics protocol

We obtain data from MD simulations using LAMMPS.[25] Since the BG is trained
on KL-loss only, the MD data is not used as training data, but as the ground
truth to compare the BG-generated configurations. The MD data was obtained in
the following way. We initialise the particles on a hexagonal lattice. We then
simulate for $1\times 10^{6}$ MD steps with a time step $\Delta t=0.001$,
saving a configuration every $100$ steps. Furthermore, we remove the first
$100$ configurations. Therefore, the MD data consists of $9900$ samples. For
more details on the MD simulations, we refer the reader to Section A.1.

####  3.1.3 Potential energy distribution

In this section, we assess the quality of the $NPT$ BG samples by comparing
their potential energy distribution to that of the MD samples, which we take
to be the ground truth. In Fig. 5, we plot the potential energy distributions
obtained from (a) an untrained and (b) a trained $NPT$ BG along with the
potential energy distribution from MD simulations. Fig. 5 shows that the
trained BG significantly improves upon the untrained BG. Furthermore, Fig.
5(b) reveals that the BG and MD potential energy distributions show good
agreement.

(a)

(b)

Figure 5: (a) The potential energy distributions corresponding to samples from
MD simulations, an untrained $NPT$ BG, and a trained $NPT$ BG. (b) The
potential energy distributions corresponding to samples from the MD
simulations and the trained $NPT$ BG show good agreement. The distributions
are normalized histograms of 30000 $NPT$ BG and 9900 MD samples.

####  3.1.4 Volume distribution

Since the pressure is fixed in the isobaric-isothermal ensemble, the volume is
allowed to fluctuate. Fig. 6(a) shows the volume distribution of the MD
samples versus an untrained and a trained $NPT$ BG. The trained $NPT$ BG
clearly improves upon the untrained $NPT$ BG. Furthermore, Fig. 6(b) shows
just the samples from the trained $NPT$ BG and from MD simulations. The two
distributions show good agreement. The $NPT$ BG slightly undersamples the
larger volumes.

(a)

(b)

Figure 6: (a) The volume distributions corresponding to samples from MD
simulations, an untrained $NPT$ BG, and a trained $NPT$ BG. (b) The potential
energy distributions corresponding to samples from MD simulations and the
trained $NPT$ BG. The distributions are normalized histograms of 30000 $NPT$
BG or 9900 MD samples.

####  3.1.5 Pressure distribution

Subsequently, we compare the instantaneous pressure of MD and $NPT$ BG
generated configurations (see Section 2.5.2 for a definition of the
instantaneous pressure). Fig. 7 shows that the distribution of the pressure of
the $NPT$ BG generated configurations has good overlap with the pressure
distribution of configurations from MD simulations. Furthermore, both
distributions are centred around the macroscopic pressure
$P\sigma^{2}/\epsilon=5.0$, as expected.

Figure 7: The instantaneous pressure distributions corresponding to samples
from MD simulations and the trained $NPT$ BG. The distributions are normalized
histograms of 30000 $NPT$ BG or 9900 MD samples.

####  3.1.6 Radial distribution function

Finally, we examine the radial distribution functions as obtained by the $NPT$
BG. Figure 8 shows that the BG radial distribution function (RDF) has peaks at
the same distances as the MD one. This indicates that the BG samples the
hexagonal crystal structure properly. However, the peaks of the BG RDF are
slightly narrower and higher. This could indicate that the BG generates
deviations that are closer to the lattice sites. Reweighting the BG generated
samples could remove this discrepancy[20].

Figure 8: The radial distribution function corresponding to samples from MD
simulations and to samples generated by the trained BG. The distributions are
generated using 30000 BG or 9900 MD samples.

###  3.2 Hemmer-Stell-like system in the $NPT$ ensemble

We now turn our attention to the second case study, which is a Hemmer-Stell-
like (HSL) system as discussed in Section A.2. We consider a 2D system with
periodic boundary conditions in the $NPT$ ensemble. The system consists of
$N=36$ particles interacting with an HSL potential[23] at a temperature
$k_{B}T/\epsilon=0.1$. This system has shown to exhibit a metastable
isostructural phase transition.[23] We train two $NPT$ BGs, one at pressure
$P\sigma^{2}/\epsilon=5.0$ and the other at pressure
$P\sigma^{2}/\epsilon=2.0$. The HSL potential is given by

$\displaystyle U_{\text{HSL}}(r)$ $\displaystyle=$ $\displaystyle
4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]$
(41) $\displaystyle-$ $\displaystyle A\exp[-w(r_{ij}-r_{0})^{2}],$

where the first term is the same as in the LJ system, and the second term
corresponds to an exponential well with parameters $A$ and $w$, centered at
$r_{0}$[23] (see Fig. 4). The settings for these parameters are listed in
Table 5.

####  3.2.1 Boltzmann generator protocol

To sample only one permutation and remove translation invariance, we add a
harmonic center-of-mass restraint in terms of the deviations $\mathbf{s}$ with
the same parameters as for the LJ system (see Eq 38). Similarly, to avoid
numerical instabilities during training, we apply the same regularization
strategy as we used for the LJ system (see Eq. 39).

black

####  3.2.2 Molecular dynamics protocol

We again collect data from MD simulations using LAMMPS.[25] The lattice
initialization and the number of MD samples is the same as in the LJ case. For
additional details on the MD simulations, we refer to Section A.2.

The BG is again trained on KL-loss only, which means that the MD data is
solely used for comparison. We also use the same RealNVP architecture as for
the LJ system, described in Section 2.2. Specifics about the architecture and
training hyperparameters can be found in Section A.2.

####  3.2.3 Potential energy distribution

Figure 9: The potential energy distributions corresponding to samples from MD
simulations, an untrained $NPT$ BG, and a trained $NPT$ BG at
$P\sigma^{2}/\epsilon=5.0$ (MD1 and BG1) and $P\sigma^{2}/\epsilon=2.0$ (MD2
and BG2). The distributions are generated using 10000 BG or 9900 MD samples.

In Fig. 9 we present the potential energy distributions obtained from MD
samples versus those of an untrained $NPT$ BG (untrained BG) and a trained BG
for both the high-pressure (MD1 and BG1) and the low-pressure (MD2 and BG2)
HSL systems. We clearly observe that the potential energy distributions
obtained from the trained BGs improve considerably with respect to the one
from the untrained BG. For both pressures, the potential energy distributions
from the trained BGs are centered at approximately the same value as their MD
counterparts. Nonetheless, the high-pressure BG shows slightly higher
potential energies than the corresponding MD sampling. This is likely because,
by being anchored to the reference lattice, the BG sampling fluctuates less
into the exponential well of the HSL potential (see Fig. 4). This mismatch can
be solved by reweighting the BG samples according to their Boltzmann
weight.[20]

####  3.2.4 Volume distribution

Figure 10: The volume distributions corresponding to samples from MD
simulations and the $NPT$ BG at $P\sigma^{2}/\epsilon=5.0$ (MD1 and BG1) and
$P\sigma^{2}/\epsilon=2.0$ (MD2 and BG2), as well as an untrained BG. The
distributions are generated using 10000 BG or 9900 MD samples.

Fig. 10 presents the volume distribution of MD samples versus those of an
untrained $NPT$ BG (untrained BG) and a trained $NPT$ BG for both the high-
pressure (MD1 and BG1) and the low-pressure (MD2 and BG2) cases. Again the
volume distributions improve significantly upon training the BG. The volume
distributions from the trained BG are centered at the same value as their
corresponding MD distributions for both pressures. However, the high-pressure
BG distribution shows better overlap with its MD counterpart compared to the
low-pressure one, i.e. the low-pressure BG distribution is not as peaked as
its MD equivalent. This is likely due to the larger fluctuations from the
reference lattice that occur at lower pressures. After this observation, the
number of training epochs of the low-pressure BG was increased (see Section
A.2). However, the improvement was marginal. Since the distributions are
considerably close, one could still reweight the BG samples to improve the
match to MD.

####  3.2.5 Pressure distribution

Figure 11: The pressure distributions corresponding to samples from MD
simulations and the $NPT$ BG at $P\sigma^{2}/\epsilon=5.0$ (MD1 and BG1) and
$P\sigma^{2}/\epsilon=2.0$ (MD2 and BG2), as well as an untrained BG. The
distributions are generated using 10000 BG or 9900 MD samples.

Fig. 11 shows the instantaneous pressure distribution of MD samples versus BG
samples for both the high-pressure (MD1 and BG1) and the low-pressure (MD2 and
BG2) systems as well as the distribution of the untrained BG. Even though the
instantaneous pressure distributions from the trained BG show reasonable
overlap with those obtained from the MD samples, the low-pressure BG
distribution is not as peaked as its MD counterpart, similarly to the volume
distributions (see Fig. 10). On the other hand, the high-pressure BG samples
slightly higher pressures than its MD counterpart, which is also reflected in
the higher potential energy (see Fig. 9). It is also worth noting that the BG
pressure distributions are considerably similar to each other in height and
width, with only the mean changing from $P\sigma^{2}/\epsilon=2.0$ to $5.0$,
whereas the MD distributions change significantly. This might indicate that
the NF mainly performs a translation of the pressure distribution.

####  3.2.6 Radial distribution function

(a)

(b)

Figure 12: (a) The radial distribution function (RDF) at high pressure, i.e.
$P\sigma^{2}/\epsilon=5.0$, corresponding to samples from MD simulations and
to samples generated by the BG. The distributions are generated using 10000 BG
or 9900 MD samples. (b) The RDF at low pressure, i.e.
$P\sigma^{2}/\epsilon=2.0$, corresponding to samples from MD simulations and
to samples generated by the BG. The distributions are generated using 10000 BG
or 9900 MD samples.

In Fig. 12 we show the RDF for (a) the high-pressure BG and (b) the low
pressure BG. Both are compared to the corresponding MD RDF. In both cases, the
consistent distance between peaks in the BG and the MD RDF shows that the
hexagonal phase is sampled properly. Moreover, the similarity between the
high-pressure and low-pressure RDFs also confirms that we sample the high- and
low-density hexagonal phases. Similar to the LJ system, the BG RDFs are
slightly more peaked than their MD counterparts. Again, this is a consequence
of the underlying reference lattice, from which the BG generates deviations.

black

##  4 CONCLUSION

In this work, we presented and tested the first BG for the $NPT$ ensemble,
which is based on previous NFs for crystal structures. [1, 29] The key idea
behind the $NPT$ BG is to generate not only particle positions, but also
fluctuations of the box itself. This adds one other dimension to the
configuration space $X$ and latent space $Z$ of the NF, requiring some
methodological adjustments as described in Section 2.3. We used the new $NPT$
BG to generate unsupervised one-shot configurations according to the isobaric-
isothermal probability distribution. These samples show good agreement with
configurations from MD simulations, as demonstrated for the LJ system.
Moreover, we show that the $NPT$ BG can sample low-density and high-density
isostructural phases, as exemplified for the HSL system. Additionally, we
derived an expression of the Gibbs free energy in terms of samples generated
by the $NPT$ BG.

While the BG samples in the $NPT$ ensemble have good agreement—especially
regarding mean values—with samples from MD, there are a few differences to
keep in mind. The main difference is observed in the distributions of
quantities such as volume or instantaneous pressure, or RDFs, which have
narrower and higher peaks in BG sampling than in MD. This is due to the
underlying reference lattice approach.[1] While this approach eliminates the
computationally prohibitive cost of sampling multiple permutations, it also
anchors the sampling to a fixed reference, hence leading to narrower sampling.
Increasing training does not seem to significantly improve the narrow
sampling. One possible solution for this is to simply reweight BG samples by
their Boltzmann weight. This is feasible because the BG distributions already
have good overlap with the MD distributions. Another possible solution is to
have a BG with a permutation invariant NF, which is further discussed below.
Another difference between BG and MD sampling is observed in the instantaneous
pressure distributions, for which the BG shows very similar distributions in
height and width for two different pressure values, implying that it learns a
translation of the pressure distribution. black

The $NPT$ BG presented in this work opens up exciting possibilities for future
research. One possible application is in the screening of pressure-driven
transitions or Gibbs free-energy differences in soft materials or atomic
systems without relying on computationally expensive MC or MD simulations.
Furthermore, the $NPT$ BG can also be used to investigate various molecular or
colloidal processes where pressure plays a critical role.

Several possible extensions to our proposed method can be considered. An
interesting potential improvement is to introduce anisotropic box
fluctuations, where $\lambda$ in $L_{y}=\lambda L_{x}$ is not constant. This
can be achieved by adding either $\lambda$ or $L_{y}$ to the generated
configuration space $X$. Similarly, this idea can also be extended to generate
the angles between the box vectors. Furthermore, the $NPT$ BG can be applied
to 3D systems with similar flexibility. These extensions will be explored in
future work.

One limitation of our approach is that it relies on a single reference
lattice. Although this allows us to avoid the combinatorial explosion
associated with permutational symmetry, it is also important to consider
transitions between two different lattice types. We discuss two potential
solutions to this problem.

One solution to overcome the limitation of our approach based on a single
reference lattice is to extend another type of BG using equivariant NFs [29,
28, 24, 15] from the $NVT$ ensemble to the $NPT$ ensemble. These BGs use a
prior distribution $\mu_{Z}$ that is invariant under a group $G$, such as the
permutation group, [28, 29] and an NF that preserves this invariance, such
that the generated distribution $q_{X}$ is invariant under $G$. This solves
the problem between BGs and permutation symmetry. However, a configuration in
the $NPT$ ensemble is given by
$(\mathbf{L},\mathbf{s}_{1},...,\mathbf{s}_{N})$, but the permutation
invariance only applies to $(\mathbf{s}_{1},...,\mathbf{s}_{N})$ and not to
$\mathbf{L}$. Hence, the challenge when extending these BGs to the $NPT$
ensemble is to construct an NF that is equivariant with respect to the last
$N*D$ coordinates alone.

The second solution we propose is to use two separate BGs rather than using
one BG to sample both phases simultaneously. A similar method was proposed by
Noé et al.. [20] Specifically, two BGs (BG1 and BG2) can be trained with two
different reference lattices to sample each of the two crystal phases. Using
Eq. 34, the absolute Gibb free energy for each phase can then be computed
separately. The Gibbs free-energy difference between the two phases can be
obtained by subtracting these two Gibbs free energies. While this method does
not provide information about the free-energy barrier, it is a powerful tool
for quickly comparing the Gibbs free energies of various phases. Additionally,
this method can be extended to obtain free-energy differences between
arbitrary combinations of interaction potentials and lattices, opening up
interesting possibilities for quickly screening free energies of soft
materials as the BG only requires modest training. Here we have focused on
training by KL-loss to use the BG in an unsupervised manner. However ML
training could also be used, even including examples from phases different
from the lattice. Pretraining models with ML-loss can help to ensure that BGs
show a competitive advantage over MC or MD simulations in terms of the number
of energy evaluations.

As a final outlook, the current $NPT$ BG could be integrated with other
existing methods for $NVT$ BGs, such as reaction coordinate (RC) loss [20] and
MC moves in latent space. [9] The RC loss can bias the BG to produce more
samples along a given transition descriptor. This loss could be used for the
$NPT$ BG to generate more samples along a phase transition. Furthermore, NFs
have been used to generate efficient update moves for MC simulations by
performing MC moves in latent space. [9] This method could be extended from
the $NVT$ ensemble to the $NPT$ ensemble by using an $NPT$ NF. Other
advancements in BGs that could be combined with the $NPT$ version include
annealed importance sampling, [19] temperature-steerable flows, [5]
conditioning for rare events, [7] and diffusion models. [12]

## ACKNOWLEDGEMENTS

We acknowledge financial support from the European Research Council (ERC
Advanced Grant No. ERC-2019-ADV-H2020 884902, SoftML).

## DATA AVAILABILITY

The code for the $NPT$ BG will be made available in a public GitHub repository
at https://github.com/MarjoleinDijkstraGroupUU.

##  Appendix A APPENDIX

black

###  A.1 Lennard-Jones system in the $NPT$ ensemble

In this section, we present details on the generation and training of the LJ
system with periodic boundary conditions in the $NPT$ ensemble. More
specifically, we provide the parameters for the LJ model and the MD
simulations, as well as the architecture and training schedule of the BG.

####  A.1.1 Boltzmann generator protocol

The parameters for the LJ potential, its cut-off distance, regularization and
the center-of-mass restraint are specified in Table 1.black

$\sigma$ | $\epsilon$ | $N$ | $P$ | $k_{B}T$ | $r_{\text{lin}}$ | $r_{cut}$ | $k_{\text{CoM}}$
---|---|---|---|---|---|---|---
1.0 | 1.0 | 36 | 5.0 | 0.5 | 0.8 | 2.5 | $3*10^{5}$
Table 1: Model parameters for the $NPT$ LJ system. All values are in reduced
units $\sigma$ and $\epsilon$, which are the particle diameter and the
interaction strength, respectively.

We train an $NPT$ BG on KL loss only. The architecture of the BG is discussed
in Section 2.2 and 2.3. In Table 2 we specify the hyperparameters used in this
architecture. Furthermore, the BG is trained over three epochs, where the
batch size is the same in each epoch and the learning rate is gradually
decreased over the epochs (see Table 3).

$n_{\text{blocks }}$ | $n_{\text{layers }}$ | $n_{\text{nodes }}$ | $C_{\text{lin}}^{\text{init}}$ | $C_{\text{scale}}^{\text{init}}$
---|---|---|---|---
12 | 3 | 300 | 6.0 | $1/(36*v_{\text{max}})$
Table 2: The parameters of the $NPT$ BG, where $n_{\text{blocks }}$ is the number of RealNVP blocks, $n_{\text{layers }}$ is the number of hidden layers in the neural networks $S$ and $T$ (Section 2.2 explains how $S$ and $T$ are used in the BG), and $n_{\text{nodes }}$ is the number of nodes for each hidden layer in $S$ and $T$. Furthermore, $C_{\text{lin}}^{\text{init}}$ and $C_{\text{scale}}^{\text{init}}$ are the initial settings of $C_{\text{scale}}$ and $C_{\text{lin}}$ (see Section 2.3). epoch | 1 | 2 | 3
---|---|---|---
iter | 4000 | 4000 | 2000
batch | 256 | 256 | 256
lr | $10^{-3}$ | $10^{-4}$ | $10^{-5}$
$w_{\text{M L}}$ | 0 | 0 | 0
$w_{\text{K L}}$ | 1 | 1 | 1
Table 3: Training schedule for the $NPT$ BG for the LJ system, where iter is
the number of iterations used for training, batch is the batch size, lr is the
learning rate used in the Adam optimiser and $w_{\text{ML}}$ and
$w_{\text{KL}}$ are the weights for the ML and KL loss, respectively. The
schedule consists of 3 epochs. Each epoch has different hyperparameters.

####  A.1.2 Molecular dynamics protocol

We generate MD data by performing simulations using LAMMPS. [25] Note that
this data is not used for training the BG, because we only train on KL loss.
It is used only for comparison. We initialise the particles on a hexagonal
lattice. This lattice has unit vectors $\mathbf{a}_{1}=(1,0)$ and
$\mathbf{a}_{2}=(0,\sqrt{3})$. There are two atoms within this unit cell. We
repeat this unit cell 6 times in the $x$-direction and 3 times in the
$y$-direction (this gives $6\times 3\times 2=36$ particles). After
initialising the system, we simulate for $10^{6}$ MD steps with time step
$\Delta t=0.001$, saving a configuration every $100$ steps. Furthermore, we
remove the first $100$ configurations. Therefore, the MD data consists of
$9900$ samples. We use a Nosé-Hoover thermostat[17] and barostat[22] to ensure
constant temperature and pressure, respectively. We use a time constant of
$\Delta t_{T}=100\Delta t$ for the thermostat and a time constant of $\Delta
t_{P}=1000\Delta t$ for the barostat. The simulation parameters are summarised
in Table 4.

$\Delta t$ | stride | MD data size | $\Delta t_{T}$ | $\Delta t_{P}$
---|---|---|---|---
0.001 | 100 | 9.900 | $100\Delta t$ | $1000\Delta t$
Table 4: MD simulation parameters for the LJ system with periodic boundary
conditions in the $NPT$ ensemble, $\Delta t$ is the time step size, stride is
the number of MD steps between saved configurations, MD data size is the
number of generated samples excluding discarded samples, and $\Delta t_{T}$
and $\Delta t_{P}$ are the time constants of the thermostat and barostat,
respectively.

###  A.2 Hemmer-Stell-like system in the $NPT$ ensemble

In this section, we present details on the generation and training of the HSL
system with periodic boundary conditions in the $NPT$ ensemble. We provide the
parameters for the HSL model and the MD simulation, as well as the
architecture and training schedule of the BG.

####  A.2.1 Boltzmann generator protocol

The parameters for the HSL potential, cut-off distance, regularization and the
center-of-mass restraint are the same as for the LJ system and specified in
Table 1. The parameters for the exponential well term as specified in Table
5.black

$A$ | $w$ | $r_{0}$
---|---|---
1.5$\epsilon$ | 41.22$\sigma^{-2}$ | 1.44$\sigma$
Table 5: Model parameters for the $NPT$ HSL system, where $A$ and $w$ are
parameters of the exponential well centered at $r_{0}$. All values are in
reduced units $\sigma$ and $\epsilon$, which are the particle diameter and the
interaction strength respectively. Values are taken from Ref. [23]. The LJ
paramters remain the same as in Table 2.

As for the LJ system, we train an $NPT$ BG on KL loss only. The architecture
of the BG is discussed in Section 2.2 and 2.3. The hyperparameters are the
same as in Table 2. The BG for $P\sigma^{2}/\epsilon=5.0$ is trained with the
schedule described in Table 6. For the BG at $P\sigma^{2}/\epsilon=2.0$, iter
is changed to 6000 for all epochs.

epoch | 1 | 2 | 3
---|---|---|---
iter | 4000 | 4000 | 4000
batch | 256 | 256 | 256
lr | $10^{-3}$ | $10^{-4}$ | $10^{-5}$
$w_{\text{M L}}$ | 0 | 0 | 0
$w_{\text{K L}}$ | 1 | 1 | 1
Table 6: Training schedule for the $NPT$ BG for the HSL system at
$P\sigma^{2}/\epsilon=5.0$, iter is the number of iterations used for
training, batch is the batch size, lr is the learning rate used in the Adam
optimiser and $w_{\text{ML}}$, and $w_{\text{KL}}$ are the weights for the ML
and KL loss respectively. The schedule consists of 3 epochs. Each epoch has
different hyperparameters. For the BG at $P\sigma^{2}/\epsilon=2.0$, iter is
changed to 6000 for all epochs.

####  A.2.2 Molecular dynamics protocol

We generate MD data by running simulations using LAMMPS. [25] The MD data is
used only for comparison against the BG data. We initialise the coordinates in
the same lattice as the LJ system. We simulate for $10^{6}$ MD steps with time
step $\Delta t=0.001$, saving a configuration every $100$ steps. Furthermore,
we remove the first $100$ configurations. Therefore, the MD data consists of
$9900$ samples. The thermostat [17] and barostat [22] settings are the same as
for the LJ case. We run two simulations, one at $P\sigma^{2}/\epsilon=2.0$ and
another at $P\sigma^{2}/\epsilon=5.0$

## References

  * [1] Rasool Ahmad and Wei Cai.  Free energy calculation of crystalline solids using normalizing flows.  Modelling and Simulation in Materials Science and Engineering, 30(6):065007, July 2022.  Publisher: IOP Publishing.
  * [2] Berni J Alder and Thomas Everett Wainwright.  Studies in molecular dynamics. I. General method.  31(2):459–466, 1959.
  * [3] J. A. Barker, D. Henderson, and F. F. Abraham.  Phase diagram of the two-dimensional Lennard-Jones system; Evidence for first-order transitions.  Physica A: Statistical Mechanics and its Applications, 106(1):226–238, March 1981.
  * [4] Ludwig Boltzmann.  Vorlesungen über Gastheorie: 2. Teil.  J.B. Barth, Leipzig, Germany, 1898.
  * [5] Manuel Dibak, Leon Klein, Andreas Krämer, and Frank Noé.  Temperature steerable flows and boltzmann generators.  Physical Review Research, 4(4):L042005, 2022.
  * [6] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using Real NVP, February 2017.  arXiv:1605.08803 [cs, stat].
  * [7] Sebastian Falkner, Alessandro Coretti, Salvatore Romano, Phillip Geissler, and Christoph Dellago.  Conditioning normalizing flows for rare event sampling.  arXiv preprint arXiv:2207.14530, 2022.
  * [8] Daan Frenkel and Berend Smit.  Understanding molecular simulation: from algorithms to applications.  Number 1 in Computational science series. Academic Press, San Diego, 2nd ed edition, 2002.
  * [9] Marylou Gabrié, Grant M. Rotskoff, and Eric Vanden-Eijnden.  Adaptive Monte Carlo augmented with normalizing flows.  Proceedings of the National Academy of Sciences, 119(10):e2109420119, March 2022.  Publisher: Proceedings of the National Academy of Sciences.
  * [10] Jérôme Hénin, Tony Lelièvre, Michael R Shirts, Omar Valsson, and Lucie Delemotte.  Enhanced sampling methods for molecular dynamics simulations.  arXiv preprint arXiv:2202.04164, 2022.
  * [11] Wei-Tse Hsu and Theodore Fobe.  Boltzmann Generators, September 2022.  original-date: 2020-02-28T04:28:14Z.
  * [12] Bowen Jing, Gabriele Corso, Jeffrey Chang, Regina Barzilay, and Tommi Jaakkola.  Torsional diffusion for molecular conformer generation.  arXiv preprint arXiv:2206.01729, 2022.
  * [13] Ivan Kobyzev, Simon J. D. Prince, and Marcus A. Brubaker.  Normalizing Flows: An Introduction and Review of Current Methods.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11):3964–3979, November 2021.  arXiv: 1908.09257.
  * [14] Andreas Kramer, Jonas Kohler, Leon Klein, Frank Noe, Michele Invernizzi, and Dibak Manuel.  bgflow, September 2022.  original-date: 2021-04-15T15:37:40Z.
  * [15] Jonas Köhler, Leon Klein, and Frank Noé.  Equivariant Flows: Exact Likelihood Generative Learning for Symmetric Densities.  arXiv:2006.02425 [physics, stat], October 2020.  arXiv: 2006.02425 version: 2.
  * [16] Alessandro Laio and Michele Parrinello.  Escaping free-energy minima.  Proceedings of the National Academy of Sciences, 99(20):12562–12566, October 2002.  Publisher: Proceedings of the National Academy of Sciences.
  * [17] Glenn J. Martyna, Michael L. Klein, and Mark Tuckerman.  Nosé–Hoover chains: The canonical ensemble via continuous dynamics.  The Journal of Chemical Physics, 97(4):2635–2643, August 1992.  Publisher: American Institute of Physics.
  * [18] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  21(6):1087–1092, 1953.
  * [19] Laurence Illing Midgley, Vincent Stimper, Gregor NC Simm, Bernhard Schölkopf, and José Miguel Hernández-Lobato.  Flow annealed importance sampling bootstrap.  arXiv preprint arXiv:2208.01893, 2022.
  * [20] Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu.  Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning.  Science, 365(6457):eaaw1147, September 2019.
  * [21] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing Flows for Probabilistic Modeling and Inference.  arXiv:1912.02762 [cs, stat], April 2021.  arXiv: 1912.02762.
  * [22] M. Parrinello and A. Rahman.  Polymorphic transitions in single crystals: A new molecular dynamics method.  Journal of Applied Physics, 52(12):7182–7190, December 1981.  Publisher: American Institute of Physics.
  * [23] D. Quigley and M. I. J. Probert.  Progression of phase behavior for a sequence of model core-softened potentials.  Physical Review E, 72(6):061202, December 2005.  Publisher: American Physical Society.
  * [24] Victor Garcia Satorras, Emiel Hoogeboom, Fabian B. Fuchs, Ingmar Posner, and Max Welling.  E(n) Equivariant Normalizing Flows.  arXiv:2105.09016 [physics, stat], June 2021.  arXiv: 2105.09016.
  * [25] Aidan P. Thompson, H. Metin Aktulga, Richard Berger, Dan S. Bolintineanu, W. Michael Brown, Paul S. Crozier, Pieter J. in ’t Veld, Axel Kohlmeyer, Stan G. Moore, Trung Dac Nguyen, Ray Shan, Mark J. Stevens, Julien Tranchida, Christian Trott, and Steven J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Computer Physics Communications, 271:108171, 2022.
  * [26] G. M. Torrie and J. P. Valleau.  Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling.  Journal of Computational Physics, 23(2):187–199, February 1977\.
  * [27] Robin Winter, Marco Bertolini, Tuan Le, Frank Noé, and Djork-Arné Clevert.  Unsupervised Learning of Group Invariant and Equivariant Representations, September 2022.  arXiv:2202.07559 [cs].
  * [28] Peter Wirnsberger, Andrew J. Ballard, George Papamakarios, Stuart Abercrombie, Sébastien Racanière, Alexander Pritzel, Danilo Jimenez Rezende, and Charles Blundell.  Targeted free energy estimation via learned mappings.  The Journal of Chemical Physics, 153(14):144112, October 2020.  Publisher: American Institute of Physics.
  * [29] Peter Wirnsberger, George Papamakarios, Borja Ibarz, Sébastien Racanière, Andrew J. Ballard, Alexander Pritzel, and Charles Blundell.  Normalizing flows for atomic solids.  arXiv:2111.08696 [cond-mat, physics:physics, stat], November 2021\.  arXiv: 2111.08696.
  * [30] Zhongwei Zhu, Mark E Tuckerman, Shane O Samuelson, and Glenn J Martyna.  Using novel variable transformations to enhance conformational sampling in molecular dynamics.  Physical review letters, 88(10):100201, 2002.
  * [31] Robert W Zwanzig.  High-temperature equation of state by a perturbation method. i. nonpolar gases.  The Journal of Chemical Physics, 22(8):1420–1426, 1954.