# Reweighted Jarzynski sampling: Acceleration of rare events and free energy
calculation with a bias potential learned from nonequilibrium work

Kristof M. Bal<EMAIL_ADDRESS>Department of Chemistry and NANOlab
Center of Excellence, University of Antwerp, Universiteitsplein 1, 2610
Antwerp, Belgium

###### Abstract

We introduce a simple enhanced sampling approach for the calculation of free
energy differences and barriers along a one-dimensional reaction coordinate.
First, a small number of short nonequilibrium simulations are carried out
along the reaction coordinate, and the Jarzynski equality is used to learn an
approximate free energy surface from the nonequilibrium work distribution.
This free energy estimate is represented in a compact form as an artificial
neural network and used as an external bias potential to accelerate rare
events in a subsequent molecular dynamics simulation. The final free energy
estimate is then obtained by reweighting the equilibrium probability
distribution of the reaction coordinate sampled under the influence of the
external bias. We apply our reweighted Jarzynski sampling recipe to four
processes of varying scales and complexities—spanning chemical reaction in the
gas phase, pair association in solution, and droplet nucleation in
supersaturated vapor. In all cases, we find reweighted Jarzynski sampling to
be a very efficient strategy, resulting in rapid convergence of the free
energy to high precision.

##  1 Introduction

The free energy $F$ is a central thermodynamic quantity that governs the
equilibrium properties and dynamics of a system. In order to calculate free
energy differences between states, one must be able to express the free energy
as a low-dimensional function $F(\mathbf{s})$, in which $\mathbf{s}$ are
collective variables that can properly distinguish configurations belonging to
any of the states of interest. To evaluate the free energy surface (FES)
$F(\mathbf{s})$ using atomistic simulations, enough configurations for all
values of $\mathbf{s}$ must be sampled. This is, however, challenging if high
free energy barriers are present on $F(\mathbf{s})$ and state-to-state
transitions are rare events.

Several techniques have been developed in the last few decades to improve
sampling and allow for an efficient reconstruction of $F(\mathbf{s})$. One
particular class is that of the so-called _adaptive_ sampling methods 1, 2, 3,
4, 5. Typically, such a method uses an external bias potential or force. This
bias is refined iteratively in such a way that the system fully explores the
whole collective variable space $\mathbf{s}$ of interest. The best possible
sampling can be achieved if the bias perfectly compensates $F(\mathbf{s})$ or
its gradient, i.e., $V(\mathbf{s})=-F(\mathbf{s})$. Therefore, a converged
bias in adaptive sampling methods is also a good estimator of the free energy.

Although the simultaneous optimization of sampling and reconstruction of the
free energy is an attractive proposition, practical application of adaptive
sampling methods can be difficult. In essence, one must find a balance between
fast exploration of configuration space (which requires rapid updates of the
bias) and a smooth convergence of the free energy (slow updates of the bias),
which may require extensive optimization of the simulation parameters 6. In
addition, even if a method’s bias is formally an estimator of the free energy,
such as in metadynamics 3 or the adaptive biasing force algorithm 2, it is not
necessarily the most efficient. Indeed, a free energy estimate from on-the-fly
sampled forces or histogram tends to converge more quickly than the bias
itself in many cases 7, 8, 9.

Convergence of the free energy requires frequent visits of all regions in
configuration space 7. In this sense, it may be better to have a roughly
reasonable bias potential as soon as possible: $V(\mathbf{s})$ need not be a
perfect estimator of $F(\mathbf{s})$, as long as it is able to induce state-
to-state transitions on an MD time scale. If sampling is good, $F(\mathbf{s})$
can then be obtained from the average force or reweighted histogram. Properly
reweighting the action of a fluctuating adaptive bias potential is however
challenging 7, 8, 9. In contrast, reweighting a trajectory obtained with a
fixed, static bias potential—a concept dating back to umbrella sampling 10—is
a much more robust approach, provided that the bias allows for sufficient
sampling. It was therefore soon realized that an adaptive method can be used
to generate $V(\mathbf{s})$, which then serves as a suitable static bias in a
subsequent reweighting 11.

Several recent sampling strategies have further explored the decoupling of the
bias and the underlying FES 12, 13, 14. A recurring concept in such methods is
that an approximate FES $\widetilde{F}$ is constructed as an interpolation of
local models of $F(\mathbf{s})$ for (well sampled) individual metastable
states. Generating an efficient sampling bias
$V(\mathbf{s})=-\widetilde{F}(\mathbf{s})$ then only entails optimizing a
limited number of interpolation parameters. Reweighting, finally, takes care
of any remaining discrepancies between $F(\mathbf{s})$ and
$\widetilde{F}(\mathbf{s})$.

In this work, we present a simple Ansatz for an efficient sampling bias that
requires no adaptive optimization before or during the sampling trajectory. In
a two-state system, characterized by a reaction coordinate $\chi$, a bias
$V(\chi)$ is learned from the work distribution of a small number of
nonequilibrium molecular dynamics trajectories along $\chi$, using the
Jarzynski equality 15. An artificial neural network (ANN) is used to represent
a smooth bias potential in a flexible manner. Although Jarzynski’s free energy
estimator is hard to converge reliably 16, 17, 18, 19, 20, 21, 22, we find
that its second order cumulant expansion can efficiently generate a bias
potential capable of driving transitions along $\chi$. A simple reweighting of
the biased trajectory suffices to reconstruct the unbiased free energy
$F(\chi)$. We demonstrate our reweighted Jarzynski sampling approach on four
systems representing gas phase chemical reactions, processes in solution, and
phase transitions. For all systems, accurate free energy differences and
barriers are obtained at a cost competitive to, or lower than, previous
enhanced sampling studies.

##  2 Theory

###  2.1 Jarzynski equality and free energy reconstruction

Suppose a transformation of interest can be described by a single collective
variable $\chi(\mathbf{R})$, in which $\mathbf{R}$ are the microscopic
coordinates of the the system. Let the system be equilibrated at a temperature
$T$, in contact with a heat bath at $T$, and defined by an initial state
$\chi=\chi_{a}$. An external force applied to $\chi$ can then be used to drive
the system to a new state $\chi_{b}$, delivering an amount of work $W$.
Jarzynski discovered the existence of a relation between the nonequilibrium
work $W$ of such a process, and the free energy difference $\Delta
F=F(\chi_{b})-F(\chi_{a})$ between the microstates $\chi_{a}$ and $\chi_{b}$
15:

$e^{-\beta\Delta F}=\langle e^{-\beta W}\rangle,$ (1)

in which $\beta=(k_{B}T)^{-1}$ and $k_{B}$ the Boltzmann constant. The average
$\langle\cdots\rangle$ indicates that from a large number of _nonequilibrium_
experiments $\chi_{a}\rightarrow\chi_{b}$, the _equilibrium_ free energy
difference $\Delta F$ can be recovered.

A practical difficulty of applying Eq. (1) for the calculation of free energy
differences, is the poor convergence of its exponential average 22. A typical
realization of the nonequilibrium experiment will require $W>\Delta F$ and
only make a small contribution to the final estimate. The average work
$\overline{W}=\Delta F+\overline{W}_{d}$, in which $\overline{W}_{d}$ is the
average dissipated work; the probability of carrying out a dominant
realization, for which $W\approx\Delta F$, is much lower. Indeed, the free
energy estimator can be very difficult to converge in the context of molecular
dynamics simulations, with unpredictable bias and variance 16, 17, 18, 19, 20,
21, 22. The large number nonequilibrium trajectories, and associated
computational cost, required for a reliable application of Eq. (1) makes its
usefulness for free energy calculation impractical. Therefore, alternative
strategies have been proposed that increase the probability of sampling
dominant, low-work trajectories or use more robust free energy estimators 23,
24, 25, 26, 27.

Within MD simulations, the transformation $\chi_{a}\rightarrow\chi_{b}$ can be
realized using steered MD (SMD), by coupling a moving Hooke spring potential
to $\chi$. Park and Schulten observed that in the limit of stiff springs and
overdamped Langevin dynamics, $W$ follows a Gaussian distribution, centered at
$\overline{W}$ 18. In this regime, the second order cumulant expansion
therefore becomes a good estimator of $\Delta F$:

$F(\chi)\approx\langle W(\chi)\rangle-\frac{\beta}{2}\left(\langle
W(\chi)^{2}\rangle-\langle W(\chi)\rangle^{2}\right),$ (2)

where we have chosen $F(\chi_{a})=0$. In principle, this expression exhibits
better convergence when only a small sample of SMD trajectories is available,
although some reports have doubted its practical accuracy 28.

###  2.2 A reweighting approach

Here, we propose to sidestep the convergence concerns regarding the
application of the Jarzynski equality to SMD simulations. We do not attempt to
directly calculate the free energy $F$ from a sample $\\{W_{i}\\}$. Rather, we
use the Jarzynski equality to learn a bias potential $V(\chi)$, and
reconstruct a free energy profile $F(\chi)$ from a reweighted molecular
dynamics trajectory on the biased potential energy surface. To wit, our
procedure is as follows:

  1. 1.

Perform a small number of SMD simulations, pulling the system from $\chi_{a}$
to $\chi_{b}$. For each trajectory $i$, record the instantaneous values of
$(\chi_{i}(t_{j}),W_{i}(t_{j}))$ at regular time intervals $t_{j}$.

  2. 2.

Learn smooth functions $W_{i}(\chi)$ from the respective set of
$(\chi_{i},W_{i})$ data. Use the cumulant expression Eq. (2) to calculate an
estimated free energy profile $\widetilde{F}(\chi)$ from $\\{W_{i}(\chi)\\}$.

  3. 3.

Perform a MD simulation while applying an external bias potential
$V(\chi)=-\widetilde{F}(\chi)$. Reconstruct the true free energy surface
$F(\chi)$ from the reweighted histogram $\langle
w(t)\cdot\delta[\chi-\chi(t)]\rangle_{b}$, in which $w(t)=e^{\beta
V(\chi(t))}$ and $\langle\cdots\rangle_{b}$ denotes a time average on the
biased potential energy surface.

In order to maximize the validity of Eq. (2) all SMD simulations should be
carried out with a strongly coupled Langevin thermostat and use spring forces
that can tightly confine $\chi$ throughout the trajectory. As a rule of thumb,
the aim was to limit fluctuations of $\chi$ to around 1% of
$|\chi_{b}-\chi_{a}|$ in all reported simulations. Still, the resulting
$(\chi_{i},W_{i})$ will be noisy. Dealing with such data is a key aspect of
machine learning techniques, and we use kernel ridge regression to learn
smooth functions $W_{i}(\chi)$ from each SMD run $i$. We can then obtain an
approximate free energy surface $\widetilde{F}(\chi)$ by inserting all
$W_{i}(\chi)$ in the cumulant expansion of Eq. (2).

We do not need $\widetilde{F}(\chi)$ to be a high-quality estimator of
$F(\chi)$. The bias potential $V(\chi)=-\widetilde{F}(\chi)$ will sufficiently
enhance sampling in the interval $[\chi_{a},\chi_{b}]$ within MD time scales
as long as it approximates the true free energy up to a few $k_{B}T$. The
curvature of $\widetilde{F}(\chi)$ is strongly dependent on the process, which
can pose a challenge when using it to calculate a biasing force. Here, too,
machine learning techniques can provide an outcome: Artificial neural networks
(ANNs) can efficiently represent arbitrary free energy surfaces in a compact
and differentiable manner 29, 30, 31, 32. We find that the one-dimensional
$\widetilde{F}(\chi)$ (or $V(\chi)$) is always well-described by an ANN with a
single hidden layer. Note that the use of an ANN is not a crucial aspect of
our approach, but is merely the most convenient option available to us that
avoids the discretization and boundary errors associated with storing
$\widetilde{F}$ on a grid.

$F(\chi)$ is defined in terms of the marginal probability density of $\chi$,
$p(\chi)=\langle\delta[\chi-\chi(t)]\rangle$:

$F(\chi)=-\frac{1}{\beta}\ln\langle\delta[\chi-\chi(t)]\rangle.$ (3)

Accordingly, if the external bias potential $V(\chi)$ is used to enhance
sampling along $\chi$, reweighting of the biased trajectory recovers $p(\chi)$
10:

$F(\chi)=-\frac{1}{\beta}\ln\langle w(t)\cdot\delta[\chi-\chi(t)]\rangle_{b},$
(4)

in which $w(t)=e^{\beta V(\chi(t))}$ and $V(\chi)=-\widetilde{F}(\chi)$, as
noted before. Because the bias potential $V(\chi)$ is static in our case, the
simple reweighting formula above suffices. Alternatively, $F$ can be recovered
directly from the biased histogram
$H(\chi)=\langle\delta[\chi-\chi(t)]\rangle_{b}$ through

$F(\chi)=\widetilde{F}(\chi)-\frac{1}{\beta}\ln\langle\delta[\chi-\chi(t)]\rangle_{b},$
(5)

which might exhibit better numerical stability if $\widetilde{F}$ is large.
The form of eq. (5) also highlights that the biased sampling run serves to
correct the imperfect nature of $\widetilde{F}(\chi)$: If
$\widetilde{F}(\chi)=F(\chi)$, $V(\chi)=-\widetilde{F}(\chi)$ is a perfect
bias potential and $H(\chi)$ will be flat.

We have restricted ourselves to a one-dimensional FES $F(\chi)$. For many
systems and transformations, good approximate reaction coordinates $\chi$ are
known and our method can be readily applied. However, even if no intuitive
reaction coordinate is available, recent data-driven approaches have made it
possible to discover candidate reaction coordinates $\chi(\mathbf{R})$ in an
automated manner 33. In one following example (sec. 3.5), we show how one such
method—harmonic linear discriminant analysis (HLDA) 34—can be integrated in
our workflow.

##  3 Examples

###  3.1 General computational details

All SMD simulations, sampling, and reweighting were carried out using PLUMED
35, 36 interfaced with different simulation codes: CP2K 37 in case of the SN2
and hydrobromination reactions and LAMMPS 38 for the nucleation process and
solvated dimer system. The temperature was controlled via Langevin
thermostats: a white noise implementation 39 in the LAMMPS-based simulations,
and a colored noise generalized Langevin thermostat 40 in CP2K, optimized for
efficient sampling 41.

All machine learning algorithms were used as implemented in the scikit-learn
library 42. Each $W_{i}(\chi)$ was interpolated with kernel ridge regression,
using a regularization strength of $0.1$ and a radial basis function (RBF)
kernel; the learned models were then used to generate a large number
$W_{i}(\chi)$ values in the interval $[\chi_{a},\chi_{b}]$ (in the order of
$10^{3}$ to $10^{4}$, depending on the system). From these data points, the
cumulant estimator of $\widetilde{F}(\chi)$ was calculated at all considered
values of $\chi$. A single hidden layer ANN with $\tanh$ activation functions
was then fitted to these $(\chi,\widetilde{F}(\chi))$ value pairs. This
training data was shuffled prior to learning the ANN, and $10$% of the data
points were held back for validation purposes. The number of neurons in the
hidden layer was chosen to minimize the error on the test set.

For each of the considered systems, only five independent SMD trajectories
were computed. Three independent reweighting simulations were carried out for
each system to assess their reproducibility. Error bars on free energy
differences or barriers are then reported at the 75 % confidence level.

Free energy differences $\Delta F$ between metastable states $A$ and $B$ are
defined as:

$\Delta F=F_{B}-F_{A}\quad\text{with}\quad
F_{X}=-\frac{1}{\beta}\int_{X}\mathrm{d}\chi\,e^{-\beta F(\chi)},$ (6)

whereas the free energy barrier of a transition $A\rightarrow B$, passing
through a dividing surface $\chi=\chi^{*}$, is

$\Delta^{\ddagger}F=F(\chi^{*})+\frac{1}{\beta}\ln\frac{\langle|\nabla\chi|\rangle^{-1}_{\chi=\chi^{*}}}{h}\sqrt{\frac{2\pi
m}{\beta}}-F_{A},$ (7)

with $h$ being the Planck constant and $m$ the mass of the reaction coordinate
43. The rate constant of the process can then be calculated within the
framework of transition state theory (TST), from the Eyring expression
$k^{\mathrm{TST}}=(h\beta)^{-1}e^{-\beta\Delta^{\ddagger}F}$. The definition
of the barrier in Eq. (7) guarantees that the choice of the dividing surface
$\chi=\chi^{*}$ that maximizes $\Delta^{\ddagger}F$ also minimizes the TST
rate $k^{\mathrm{TST}}$, which is in turn an upper bound on the true rate $k$.
A gauge correction, function of the norm of the gradient of $\chi$ with
respect to all atomic coordinates $\mathbf{R}$, is needed to obtain barriers
invariant to the functional form of $\chi(\mathbf{R})$.

Sample inputs and scripts to reproduce the reported simulations are deposited
on PLUMED-NEST, the public repository of the PLUMED consortium 36, as
plumID:21.020 44.

###  3.2 Symmetric nucleophilic substitution reaction

Figure 1:  Symmetric S${}_{\mathrm{N}}2$ reaction of methyl chloride. (a)
Definition of the distances $d_{1}$ and $d_{2}$. (b) $(\chi,W)$ of two
independent SMD runs, and the final ANN representation of
$\widetilde{F}(\chi)$. (c) Example trajectory of $\chi$ during a reweighting
run using $-\widetilde{F}(\chi)$ as external bias potential. (d) Final
estimates $F(\chi)$ from three independent reweighting runs, compared to
$\widetilde{F}(\chi)$. The curves are vertically aligned at $\chi=0$ Å.

As a first illustration of reweighted Jarzynski sampling, we use the symmetric
S${}_{\mathrm{N}}2$ reaction of methyl chloride,
$\mathrm{Cl^{-}}+\mathrm{CH_{3}Cl}\rightarrow\mathrm{ClCH_{3}}+\mathrm{Cl^{-}}$,
a common testbed for enhanced sampling techniques 45, 46, 47, 5, 48.
Intuitively, we choose $\chi=(d_{2}-d_{1})/\sqrt{2}$ (Figure 1a). We apply
harmonic restraints to keep $d_{1}$ and $d_{2}$ below 5 Å, describe the
interatomic interactions with the semi-empirical PM6 Hamiltonian 49, and
integrate the equations of motion with a time step of 0.5 fs. SMD simulations
were then performed in the interval $\chi\in[-1.5,1.5]$ Å at $T=300$ K, over a
time scale of 10 ps.

The $(\chi(t),W(t))$ trajectory of two independent SMD simulations can be
quite different (up to about $2k_{B}T$), as shown in Figure 1b. Even though
the sampled data is rather noisy (using a spring constant of $10^{3}$ kJ mol-1
Å-2), the final ANN approximation of $\widetilde{F}(\chi)$ is however a smooth
curve, for which 12 hidden neurons suffice. Because of the symmetry of the
problem, we can also see that $\widetilde{F}(\chi)$ is not a perfect estimator
of the true FES, as it is biased towards the $\chi<0$ state by about 3 kJ/mol.
Such a bias is usually hard to distinguish from the variance or other sources
of error 17, 19. The asymmetry of the putative FES also affects a reweighting
run using $V(\chi)=-\widetilde{F}(\chi)$ as a bias potential, which
preferentially samples the $\chi>0$ state (Figure 1c). The bias potential
does, however, also efficiently induce many state-to-state crossings, and an
extensive sampling of the transition state region at $\chi\approx 0$. As a
result, the final estimates of $F(\chi)$ after 100 ps reweighting runs are, on
average, free from the initial bias in $\widetilde{F}(\chi)$ (Figure 1d): From
three trajectories, we find $\Delta F=0.2\pm 1.1$ kJ/mol.

A very modest sampling time of $5\times 10$ ps (SMD runs) + 100 ps
(reweighting) can thus yield reaction free energies with rather good
precision, and repeated reweighting runs show that the uncertainty is well
below $k_{B}T$. This does not only hold true for thermodynamics, but also
kinetics. The free energy barrier of the forward reaction (from low to high
$\chi$) is predicted to be $\Delta^{\ddagger}F=36.7\pm 0.8$ kJ/mol,
corresponding to a rate $k^{\mathrm{TST}}=(2.6\pm 0.8)\times 10^{6}$ s-1. This
implies an average reaction time of $390\pm 120$ ns, which compares well to
the range 180–300 ns obtained by Fu et al. using metadynamics in an almost
identical model system 47. TST thus appears to describe this particular
process rather well.

The ability of the reweighting procedure to quickly and transparently remove
the bias error that may be present in $\widetilde{F}$ is greatly facilitated
by the static nature of the bias. Whenever many state-to-state transitions are
observed during the reweighting run, in principle the only remaining error in
$F$ will be due to variance. The size of this variance can then be estimated
with the usual techniques, such as block averaging or repeated simulations.

###  3.3 Droplet nucleation in supersaturated vapor

As an example of a process involving diffusive collective motion of a large
number of atoms, we revisit the condensation of Ar vapor at $T=80.7$ K, for
which rates have been calculated previously 50, 51. We model 512 atoms in a
cubic simulation cell of length 11.5 nm, corresponding to an initial
supersaturation of 8.68. Interatomic interaction are described by a Lennard-
Jones potential, using the same parameters as earlier work 50. We use the
number of liquidlike atoms $n$ as a reaction coordinate, that is, the number
of atoms that have more than five nearest neighbors. Equations of motion are
integrated with a time step of 5 fs. Because we are only interested in the
formation of the critical nucleus, we perform SMD from $n=0$ to $n=64$ over 10
ns, and also place a harmonic wall at $n=64$ (Figure 2a).

Figure 2:  Droplet nucleation from supersaturated Ar vapor. (a) Snapshot of a
configuration corresponding to $n\approx 64$. (b) Example $(n,W)$ curves of
three independent SMD runs, and the final ANN representation of
$\widetilde{F}(n)$. (c) Final estimates $F(n)$ from three independent
reweighting runs, compared to $\widetilde{F}(n)$. The curves are vertically
aligned at $n=0$.

The individual $(n(t),W(t))$ trajectories can differ quite strongly for the
process, and are also rather rough (Figure 2b). $\widetilde{F}(n)$ is however
again a very smooth function owing to the regularization employed in the
kernel ridge regression interpolation, and is well represented by only 12
neurons in the hidden layer. From three 100 ns reweighting runs we can see
that: (1) $\widetilde{F}(n)$ mostly recovers the shape of the FES, and while
it underestimates the nucleation barrier by more than 2 $k_{B}T$ (2) the
reweighted $F(n)$ reliably converges within this time scale in a very
consistent manner (Figure 2c).

Nucleation barriers are recovered with very high precision, yielding the
estimate $\Delta^{\ddagger}F=7.06\pm 0.15$ kJ/mol. Note that in our previous
study of this process 51, a metadynamics simulation required 1 $\mu$s to
produce a barrier estimate with a larger uncertainty, namely, $7.46\pm 0.43$
kJ/mol (error calculated from four 200 ns trajectory chunks, adjusted to the
same confidence level of 75 %). A correct rate estimate also requires
knowledge of the transmission coefficient, which is $\kappa=(2.1\pm 0.5)\times
10^{-3}$ at this particular supersaturation level 51. Any rate estimate
$k=\kappa k^{\mathrm{TST}}$ is only valid within the simulation cell volume
$V$. A global nucleation rate should therefore be calculated as $J=k/V$ 50,
51. From this data, we can estimate the nucleation rate to be $J=(6.3\pm
2.1)\times 10^{22}$ cm${}^{-3}~{}s^{-1}$, close to our previously reported
value of $(3.5\pm 2.4)\times 10^{22}$ cm${}^{-3}~{}s^{-1}$. The overall
computational cost required in this work is however three times lower.

###  3.4 Dissociation and association in solution

Processes in solution depend on an intricate interplay between solvent and
solute, which makes their analysis non-trivial 52. We introduced a simple
example of such a system in earlier work, consisting of a Morse dimer solvated
in a dense Lennard–Jones fluid 43.

Reduced units were used throughout. We performed SMD simulations on this
system at a reduced temperature of $T=1.5$, manipulating the interatomic
distance in the dimer from $r=0.8$ to $r=8.0$, over a time scale of
$10^{3}\tau$ ($5\times 10^{5}$ time steps). We again fit $\widetilde{F}(r)$ to
a {12} ANN and converge $F(r)$ in reweighting runs of $10^{4}\tau$, depicted
in Figure 3.

Figure 3:  Free energy of a strongly interacting dimer, solvated by a Lennard-
Jones fluid, as a function of the interatomic separation $r$. The resultant
$F(r)$ of three independent reweighting runs are depicted in different colors.

$F(r)$ is highly asymmetric, and it seems that narrow deep free energy wells
(i.e., the bound state) are better-sampled than wide, shallow states (the
dissociated state). This is also reflected by the unbinding free energy
barrier $\Delta^{\ddagger}F=4.69\pm 0.07$, which is determined with a higher
precision than the unbinding free energy, $\Delta F=1.17\pm 0.33$. These
estimates are in agreement with previously calculated values, which are
$\Delta^{\ddagger}F=4.83\pm 0.14$ and $\Delta F=1.39\pm 0.16$, respectively
43.

###  3.5 Four-center addition reaction

It is not always trivial to find an appropriate one-dimensional reaction
coordinate $\chi$. For example, the hydrobromination reaction of propene, HBr
+ C3H6, involves forming a C–Br bond, forming a C–H bond, breaking the H–Br
bond, and changing the bonding character of the vinyl moiety. In addition, two
possible products can be formed: the dominant 2-bromopropane (the Markovnikov
product) and the less likely 1-bromopropane (the anti-Markovnikov product).
Due to this manybody nature, the reaction has recently been recognized as an
interesting test case for targeted sampling of reaction pathways. Methods that
aim to explore high-dimensional configuration spaces have been tested on this
system 14, as well as dimensionality reduction techniques 53.

Figure 4:  Hydrobromination of propene. (a) Definition of the six degrees of
freedom included in the reaction coordinate $\chi(\mathbf{R})$. (b) Free
energy surface $F(\chi)$ of the formation of the Markovnikov product
2-bromopropane. The resultant $F(\chi)$ of three independent reweighting runs
are depicted in different colors and line styles.

Here, we wish to specifically target the Markovnikov pathway. For this
purpose, we must identify a suitable reaction coordinate $\chi(\mathbf{R})$. A
method that integrates rather well in our established workflow, is harmonic
linear discriminant analysis (HLDA) 34. In HLDA, one first identifies the
metastable states of interest, which are here the HBr + propene system and
2-bromopropane, respectively. A set of simple microscopic descriptors that may
be able to discriminate between these states must then be selected. In order
to exploit chemical intuition as little as possible, we here opt to use all
internal degrees of freedom of the four reaction centers, namely, the six
interatomic distances between H, Br and the vinylic C atoms (Figure 4a).
Finally, short MD runs are carried out in the two metastable state of interest
so that for each state $i$, the mean $\bm{\mu}_{i}$ and multivariate variance
$\bm{\Sigma}_{i}$ of the descriptors $\mathbf{d}=(d_{1},d_{2},\ldots,d_{6})$
can be calculated. Finally, the coordinate $\chi$ that best separates the two
states $A$ and $B$ is calculated as a linear combination of the candidate
descriptors:

$\chi=(\bm{\mu}_{A}-\bm{\mu}_{B})^{T}\left(\frac{1}{\bm{\Sigma}_{A}}+\frac{1}{\bm{\Sigma}_{B}}\right)\mathbf{d},$
(8)

in a modification of Fisher’s linear discriminant analysis (LDA). An advantage
of HLDA is that it only requires us to know a priori the metastable states of
interest, without having to sample any $A\rightarrow B$ transition pathway
first.

As in previous studies 53, 14, harmonic restraints were used to keep all
components of $\mathbf{d}$ below 3 Å, interatomic forces were calculated at
the semi-empirical PM6 level of theory 49, and sampling was carried out at
$T=300$ K. The integration time step was 0.5 fs. We first parameterized
$\chi=0.661d_{1}-0.656d_{2}-0.328d_{3}+0.011d_{4}-0.021d_{5}+0.157d_{6}$ from
20 ps sampling runs in the metastable states. Then, SMD runs were carried out
over 10 ps from $\chi=-2$ Å and $\chi=1$ Å. A larger {48} ANN was found to be
needed for an accurate representation of the sharply varying
$\widetilde{F}(\chi)$. Finally, reweighting over just 100 ps yields very
consistent estimates of the FES $F(\chi)$ (Figure 4b).

The predicted reaction free energy $\Delta F=-33.8\pm 1.9$ kJ/mol, and the
barrier $\Delta^{\ddagger}F=131.7\pm 0.8$ kJ/mol. This very high free energy
barrier implies that the reaction is very slow: $k^{\mathrm{TST}}=(7.46\pm
2.34)\times 10^{-11}$ s-1. However, the bias potential $-\widetilde{F}(\chi)$
is of sufficient quality to induce several transitions in both directions
within short MD runs.

Due to its high barriers and asymmetric FES this particular system constitutes
the most challenging application of reweighted Jarzynski sampling up till now.
It is therefore an instructive test case for comparing the efficiency of our
sampling recipe to state-of-the-art adaptive techniques. For this purpose, we
choose metadynamics, 3 variationally enhanced sampling (VES), 4 and the on-
the-fly probability-enhanced sampling (OPES) method 54. A fair comparison is
easiest to make with OPES due to its limited number of parameters: Only the
bias update frequency and an estimate of the highest barrier is required for
it to be applied. We can stack the deck strongly in the favor of OPES because
we already have an estimate of the barriers by now: We set the OPES barrier
parameter to 200 kJ/mol, well above either barrier of interest, and update the
bias in 100 fs intervals. Due to their larger number of user-set parameters, a
straightforward comparison to metadynamics and VES is somewhat more tricky. We
opt to use a rather aggressive biasing strategy (along the lines of previous
metadynamics simulations on the system 53): A hill height/bias update stepsize
of 5 kJ/mol, a well-tempered bias factor 55, 56 of $\gamma=50$, and bias
update stride of 100 fs.

Figure 5:  Comparing sampling efficiency of different methods in the
hydrobromination reaction. Trajectory of $\chi$ in (a) a reweighting run using
$-\widetilde{F}(\chi)$ as static bias and adaptive runs with (b) OPES, (c)
metadynamics, and (d) VES. Note that the depicted OPES trajectory is five
times longer.

We plot the trajectory of $\chi$ in a reweighting run under influence of the
bias $-\widetilde{F}(\chi)$ as well as during the various dynamic biasing runs
(Figure 5). As noted earlier, the Jarzynski-derived $\widetilde{F}(\chi)$ is a
good model for both barriers and leads to nearly uniform sampling from the
get-go (Figure 5a). This is not true for OPES (Figure 5b). The initial
$A\rightarrow B$ transition occurs rather quickly in the OPES run, because its
barrier height is approximately known and an efficient bias can be easily
constructed. Then, however, the simulation remains stuck in the $B$ state: The
OPES algorithm has no way of knowing the free energy difference $\Delta F$
between the metastable states or, equivalently, the barrier height of the
$B\rightarrow A$ transition. The bias for escaping state $B$ is therefore only
accumulated slowly. Similar behaviour is also observed with slower bias update
frequencies (1 ps) and a lower barrier parameter (150 kJ/mol).

Table 1: Comparison of Hydrobromination Free Energy Estimates from Different
Sampling Methods _a_

Method |  $t$ (ps) |  $\Delta F$ (kcal/mol) |  $\Delta^{\ddagger}F$ (kcal/mol)
---|---|---|---
This work | 100 | $-33.8\pm 1.9$ | $131.7\pm 0.8$
Metadynamics | 100 | $4.0\pm 0.7$ | $153.9\pm 3.8$
| 500 | $-22\pm 14$ | $134.1\pm 1.8$
VES | 100 | $-56.8\pm 4.9$ | $124\pm 16$
| 500 | $-37\pm 13$ | $125.7\pm 6.6$

_a_ Reaction free energy $\Delta F$ and barrier $\Delta^{\ddagger}F$ from
reweighting runs of varying length $t$.

Metadynamics and VES fare better (Figs. 5c–d). Both algorithms can reliably
drive transitions within 100 ps, which is about on par with the Jarzynski
bias. We can therefore also compare the quality of their respective free
energy estimates. For fair comparison we use reweighting techniques 57, 4 to
reconstruct the free energy estimates also for these methods. Three
independent simulations are performed with each algorithm. Reweighting the
action of a fluctuating bias potential is more difficult 7, 8, 9, as can be
seen in Table 1. A 100 ps sampling time is not sufficient to produce
quantitatively correct free energy estimates. Longer 500 ps sampling runs
improve the quality of $\Delta F$ and $\Delta^{\ddagger}F$, but still with a
substantially higher variance than reweighted Jarzynski sampling; better
reweighting approaches or biasing strategies might however improve the
performance of metadynamics and VES.

###  3.6 Limitations and future prospects

Although the good performance of reweighted Jarzynski sampling in the above
examples is very encouraging, we must address two obvious limitations of the
approach.

First, the approach outlined here is intended for the targeted calculation of
free energy differences and barriers of transformations along a single
reaction coordinate. The reaction path of interest must be known a priori. The
study of high-dimensional free energy surfaces is therefore not the main
application domain. Nor is the unconstrained exploration of configuration
space: To characterize previously unknown states and pathways, one can use
specialized “blind” methods 58, 59, 14, 60.

Figure 6:  The ANN $\widetilde{F}$ approximates $F$ up to a few $k_{B}T$ in
the considered systems. Using the the bias $V=-\widetilde{F}$, the FES is
almost fully flattened, as demonstrated by the biased histogram
$H=\langle\delta[\chi-\chi(t)]\rangle_{b}$. The contributions to the FES
estimate $F(\chi)=-V(\chi)-k_{B}T\ln H(\chi)$ are shown in units of $k_{B}T$
for (a) the symmetric SN2 reaction, (b) droplet nucleation, (c) the solvated
dimer, and (d) the hydromination reaction.

Second, although $-\widetilde{F}$ has turned out to be a perfectly cromulent
bias potential in the above examples, this might not be generally true. As
noted earlier, it is difficult to consistently eliminate bias and variance
errors in applications of the Jarzynski equality 17, 19. If $\widetilde{F}$
deviates from $F$ by more than a few $k_{B}T$, the resultant bias may not lead
to sufficient sampling. The above examples therefore represent more or less
ideal scenarios: the ANN $\widetilde{F}(\chi)$ is a good model of $F(\chi)$,
and the biased histogram $H(\chi)=\langle\delta[\chi-\chi(t)]\rangle_{b}$ is
nearly flat (Figure 6). The correction $-k_{B}T\ln H(\chi)$ to
$\widetilde{F}(\chi)$ (following Eq. (5)) is therefore always only a few
$k_{B}T$. Because it is however not a priori obvious what an optimal SMD set
up is for a specific system, such a high-quality $\widetilde{F}$ may not
possible to obtain reliably for all processes. One possible way to make
$-\widetilde{F}$ an even more robust bias would be to use methods with
improved convergence 23, 24, 25, 26, 27. In addition, the toolbox of
nonequilibrium nanoscale thermodynamics 61 could be further explored: For
example, if also reverse trajectories ($B\rightarrow A$) are sampled, the
Crooks fluctuation theorem 62 can be invoked as an alternative means of
estimating free energies.

Even if $\widetilde{F}$ is a poor approximation of $F$ it can still serve as
an initial step in an enhanced sampling procedure: An adaptive sampling method
like metadynamics, variationally enhanced sampling, or OPES could operate on
top of a preconditioned, mollified FES $F-\widetilde{F}$. Such preconditioning
may also help prevent the adaptive method from getting stuck in very deep
minima. Alternatively, following other recent strategies, our approach could
be the basis of an iterative procedure, in which the initial estimate
$\widetilde{F}$ is further refined in a sequence of sampling runs until
convergence is reached 30, 31, 54—a strategy that dates back to at least
adaptive umbrella sampling 1.

Currently, we have only employed a _parallel_ strategy to improve the estimate
of $F$, in the sense that multiple independent reweighting runs were carried
out. If we observed that several transitions occurred within each reweighting
run, we assumed that the bias error in $F$ was small, so that we could
directly quantify the inherent variance of the FES estimates by comparing the
reweighting simulations. If either check proved unsatisfactory for our trial
runs, we simply opted to carry out the reweighting runs for more steps. We
also did not explicitly asses the quality of the employed reaction
coordinates: A poorly chosen $\chi$ negatively affects the performance of
adaptive sampling methods, and this is equally true for our approach.

Finally, we would like to note one very practical convenience of the
reweighted Jarzynski sampling approach as described here: All required
simulation steps can be readily carried out with unmodified versions of
standard sampling codes. The only requirement is that moving harmonic
restraints and the collective variables of choice are implemented, and that
arbitrary bias potentials are allowed. Recent versions of PLUMED 35, 36
possess these features.

##  4 Conclusions

We have introduced a recipe for the calculation of a one-dimensional free
energy surface (FES) from atomistic simulations. The final FES estimate
$F(\chi)$ is obtained from the reweighted _equilibrium_ probability density of
the coordinate $\chi$, sampled under influence of an external bias potential
$V(\chi)$. The bias $V(\chi)$, in turn, is learned first from a small
collection of _nonequilibrium_ trajectories by employing the Jarzynski
equality. $V(\chi)$ is smoothly and compactly represented in the form of an
artificial neural network (ANN).

We find reweighted Jarzynski sampling to exhibit excellent efficiency and
accuracy in the considered systems. Free energy differences and barriers with
error bars well below $k_{B}T$ are routinely obtained within very economical
simulation time scales, and accurate kinetics are simultaneously accessible
through a recent numerical framework 43, 51.

The approach is also compatible with adaptively biased sampling methods so
that a joint application may exhibit faster convergence. The method can
readily use any type of collective variable as reaction coordinate $\chi$, and
only requires an off-the-shelf version of the PLUMED plugin or any other free
energy code with similar capabilities.

## Acknowledgments

K.M.B. was funded as a junior postdoctoral fellow of the FWO (Research
Foundation – Flanders), Grant 12ZI420N. The computational resources and
services used in this work were provided by the HPC core facility CalcUA of
the Universiteit Antwerpen, and VSC (Flemish Supercomputer Center), funded by
the FWO and the Flemish Government. HLDA calculations were performed with a
script provided by GiovanniMaria Piccini. K.M.B. thanks Erik Neyts for
proofreading this manuscript, the many discussions over the years, and for his
continuous support.

## References

  * Mezei 1987 Mezei, M. Adaptive umbrella sampling: Self-consistent determination of the non-Boltzmann bias. _J. Comput. Phys._ 1987, _68_ , 237–248
  * Darve and Pohorille 2001 Darve, E.; Pohorille, A. Calculating free energies using average force. _J. Chem. Phys._ 2001, _115_ , 9169–9183
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proc. Natl. Acad. Sci._ 2002, _99_ , 12562–12566
  * Valsson and Parrinello 2014 Valsson, O.; Parrinello, M. Variational Approach to Enhanced Sampling and Free Energy Calculations. _Phys. Rev. Lett._ 2014, _113_ , 090601
  * Zhang et al. 2019 Zhang, J.; Yang, Y. I.; Noé, F. Targeted Adversarial Learning Optimized Sampling. _J. Phys. Chem. Lett._ 2019, _10_ , 5791–5797
  * Dama et al. 2014 Dama, J. F.; Rotskoff, G.; Parrinello, M.; Voth, G. A. Transition-Tempered Metadynamics: Robust, Convergent Metadynamics via On-the-Fly Transition Barrier Estimation. _J. Chem. Theory Comput._ 2014, _10_ , 3626–3633
  * Cuendet and Tuckerman 2014 Cuendet, M. A.; Tuckerman, M. E. Free Energy Reconstruction from Metadynamics or Adiabatic Free Energy Dynamics Simulations. _J. Chem. Theory Comput._ 2014, _10_ , 2975–2986
  * Mones et al. 2016 Mones, L.; Bernstein, N.; Csányi, G. Exploration, Sampling, And Reconstruction of Free Energy Surfaces with Gaussian Process Regression. _J. Chem. Theory Comput._ 2016, _12_ , 5100–5110
  * Marinova and Salvalaglio 2019 Marinova, V.; Salvalaglio, M. Time-independent free energies from metadynamics via mean force integration. _J. Chem. Phys._ 2019, _151_ , 164115
  * Torrie and Valleau 1977 Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. _J. Comput. Phys._ 1977, _23_ , 187–199
  * Babin et al. 2006 Babin, V.; Roland, C.; Darden, T. A.; Sagui, C. The free energy landscape of small peptides as obtained from metadynamics with umbrella sampling corrections. _J. Chem. Phys._ 2006, _125_ , 204909
  * McCarty et al. 2016 McCarty, J.; Valsson, O.; Parrinello, M. Bespoke Bias for Obtaining Free Energy Differences within Variationally Enhanced Sampling. _J. Chem. Theory Comput._ 2016, _12_ , 2162–2169
  * Invernizzi and Parrinello 2019 Invernizzi, M.; Parrinello, M. Making the Best of a Bad Situation: A Multiscale Approach to Free Energy Calculation. _J. Chem. Theory Comput._ 2019, _15_ , 2187–2194
  * Debnath and Parrinello 2020 Debnath, J.; Parrinello, M. Gaussian Mixture-Based Enhanced Sampling for Statics and Dynamics. _J. Phys. Chem. Lett._ 2020, _11_ , 5076–5080
  * Jarzynski 1997 Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. _Phys. Rev. Lett._ 1997, _78_ , 2690–2693
  * Hummer 2001 Hummer, G. Fast-growth thermodynamic integration: Error and efficiency analysis. _J. Chem. Phys._ 2001, _114_ , 7330–7337
  * Gore et al. 2003 Gore, J.; Ritort, F.; Bustamante, C. Bias and error in estimates of equilibrium free-energy differences from nonequilibrium measurements. _Proc. Natl. Acad. Sci. U.S.A._ 2003, _100_ , 12564–12569
  * Park and Schulten 2004 Park, S.; Schulten, K. Calculating potentials of mean force from steered molecular dynamics simulations. _J. Chem. Phys._ 2004, _120_ , 5946–5961
  * Shirts and Pande 2005 Shirts, M. R.; Pande, V. S. Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration. _J. Chem. Phys._ 2005, _122_ , 144107
  * Lua and Grosberg 2005 Lua, R. C.; Grosberg, A. Y. Practical Applicability of the Jarzynski Relation in Statistical Mechanics: A Pedagogical Example. _J. Phys. Chem. B_ 2005, _109_ , 6805–6811
  * Oberhofer et al. 2005 Oberhofer, H.; Dellago, C.; Geissler, P. L. Biased Sampling of Nonequilibrium Trajectories: Can Fast Switching Simulations Outperform Conventional Free Energy Calculation Methods? _J. Phys. Chem. B_ 2005, _109_ , 6902–6915
  * Jarzynski 2006 Jarzynski, C. Rare events and the convergence of exponentially averaged work values. _Phys. Rev. E_ 2006, _73_ , 046105
  * Ytreberg and Zuckerman 2004 Ytreberg, F. M.; Zuckerman, D. M. Efficient use of nonequilibrium measurement to estimate free energy differences for molecular systems. _J. Comput. Chem._ 2004, _25_ , 1749–1759
  * Echeverria and Amzel 2012 Echeverria, I.; Amzel, L. M. Estimation of Free-Energy Differences from Computed Work Distributions: An Application of Jarzynski’s Equality. _J. Phys. Chem. B_ 2012, _116_ , 10986–10995
  * Wolf and Stock 2018 Wolf, S.; Stock, G. Targeted Molecular Dynamics Calculations of Free Energy Profiles Using a Nonequilibrium Friction Correction. _J. Chem. Theory Comput._ 2018, _14_ , 6175–6182
  * Arrar et al. 2019 Arrar, M.; Boubeta, F. M.; Szretter, M. E.; Sued, M.; Boechi, L.; Rodriguez, D. On the accurate estimation of free energies using the Jarzynski equality. _J. Comput. Chem._ 2019, _40_ , 688–696
  * Roussey and Dickson 2020 Roussey, N. M.; Dickson, A. Enhanced Jarzynski free energy calculations using weighted ensemble. _J. Chem. Phys._ 2020, _153_ , 134116
  * Noh and Notman 2020 Noh, S. Y.; Notman, R. Comparison of umbrella sampling and steered molecular dynamics methods for computing free energy profiles of aromatic substrates through phospholipid bilayers. _J. Chem. Phys_ 2020, _153_ , 034115
  * Schneider et al. 2017 Schneider, E.; Dai, L.; Topper, R. Q.; Drechsel-Grau, C.; Tuckerman, M. E. Stochastic Neural Network Approach for Learning High-Dimensional Free Energy Surfaces. _Phys. Rev. Lett._ 2017, _119_ , 150601
  * Galvelis and Sugita 2017 Galvelis, R.; Sugita, Y. Neural Network and Nearest Neighbor Algorithms for Enhancing Sampling of Molecular Dynamics. _J. Chem. Theory Comput._ 2017, _13_ , 2489–2500
  * Sidky and Whitmer 2018 Sidky, H.; Whitmer, J. K. Learning free energy landscapes using artificial neural networks. _J. Chem. Phys._ 2018, _148_ , 104111
  * Bonati et al. 2019 Bonati, L.; Zhang, Y.-Y.; Parrinello, M. Neural networks-based variationally enhanced sampling. _Proc. Natl. Acad. Sci. U.S.A._ 2019, _116_ , 17641–17647
  * Sidky et al. 2020 Sidky, H.; Chen, W.; Ferguson, A. L. Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation. _Mol. Phys._ 2020, _118_ , e1737742
  * Mendels et al. 2018 Mendels, D.; Piccini, G.; Parrinello, M. Collective Variables from Local Fluctuations. _J. Phys. Chem. Lett._ 2018, _9_ , 2776–2781
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Comput. Phys. Commun._ 2014, _185_ , 604–613
  * The PLUMED consortium 2019 The PLUMED consortium, Promoting Transparency and Reproducibility in Enhanced Molecular Simulations. _Nat. Methods_ 2019, _16_ , 670–673
  * Kühne et al. 2020 Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F.; Golze, D.; Wilhelm, J.; Chulkov, S.; Bani-Hashemian, M. H.; Weber, V.; Borštnik, U.; Taillefumier, M.; Jakobovits, A. S.; Lazzaro, A.; Pabst, H.; Müller, T.; Schade, R.; Guidon, M.; Andermatt, S.; Holmberg, N.; Schenter, G. K.; Hehn, A.; Bussy, A.; Belleflamme, F.; Tabacchi, G.; Glöß, A.; Lass, M.; Bethune, I.; Mundy, C. J.; Plessl, C.; Watkins, M.; VandeVondele, J.; Krack, M.; Hutter, J. CP2K: An electronic structure and molecular dynamics software package – Quickstep: Efficient and accurate electronic structure calculations. _J. Chem. Phys._ 2020, _152_ , 194103
  * Plimpton 1995 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. _J. Comput. Phys._ 1995, _117_ , 1–19
  * Bussi and Parrinello 2007 Bussi, G.; Parrinello, M. Accurate sampling using Langevin dynamics. _Phys. Rev. E_ 2007, _75_ , 056707
  * Ceriotti et al. 2009 Ceriotti, M.; Bussi, G.; Parrinello, M. Langevin Equation with Colored Noise for Constant-Temperature Molecular Dynamics Simulations. _Phys. Rev. Lett._ 2009, _102_ , 020601
  * Ceriotti et al. 2010 Ceriotti, M.; Bussi, G.; Parrinello, M. Colored-Noise Thermostats à la Carte. _J. Chem. Theory Comput._ 2010, _6_ , 1170–1180
  * Pedregosa et al. 2011 Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E. Scikit-learn: Machine Learning in Python. _J. Mach. Learn. Res._ 2011, _12_ , 2825–2830
  * Bal et al. 2020 Bal, K. M.; Fukuhara, S.; Shibuta, Y.; Neyts, E. C. Free energy barriers from biased molecular dynamics simulations. _J. Chem. Phys._ 2020, _153_ , 114118
  * Bal 2021 Bal, K. M. Reweighted Jarzynski sampling. 2021; https://www.plumed-nest.org/eggs/21/020, PLUMED-NEST, plumID:21.020
  * Ensing et al. 2005 Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions. _J. Phys. Chem. B_ 2005, _109_ , 6676–6687
  * Fleming et al. 2016 Fleming, K. L.; Tiwary, P.; Pfaendtner, J. New Approach for Investigating Reaction Dynamics and Rates with Ab Initio Calculations. _J. Phys. Chem. A_ 2016, _120_ , 299–305
  * Fu et al. 2017 Fu, C. D.; Oliveira, L. F. L.; Pfaendtner, J. Assessing Generic Collective Variables for Determining Reaction Rates in Metadynamics Simulations. _J. Chem. Theory Comput._ 2017, _13_ , 968–973
  * Bal et al. 2020 Bal, K. M.; Bogaerts, A.; Neyts, E. C. Ensemble-Based Molecular Simulation of Chemical Reactions under Vibrational Nonequilibrium. _J. Phys. Chem. Lett._ 2020, _11_ , 401–406
  * Stewart 2007 Stewart, J. J. P. Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. _J. Mol. Model._ 2007, _13_ , 1173–1213
  * Salvalaglio et al. 2016 Salvalaglio, M.; Tiwary, P.; Maggioni, G. M.; Mazzotti, M.; Parrinello, M. Overcoming time scale and finite size limitations to compute nucleation rates from small scale well tempered metadynamics simulations. _J. Chem. Phys._ 2016, _145_ , 211925
  * Bal 2021 Bal, K. M. Nucleation rates from small scale atomistic simulations and transition state theory. 2021
  * Mullen et al. 2014 Mullen, R. G.; Shea, J.-E.; Peters, B. Transmission Coefficients, Committors, and Solvent Coordinates in Ion-Pair Dissociation. _J. Chem. Theory Comput._ 2014, _10_ , 659–667
  * Piccini et al. 2018 Piccini, G.; Mendels, D.; Parrinello, M. Metadynamics with Discriminants: A Tool for Understanding Chemistry. _J. Chem. Theory Comput._ 2018, _14_ , 5040–5044
  * Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: From Bias Potentials to Probability Distributions. _J. Phys. Chem. Lett._ 2020, _11_ , 2731–2736
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. _Phys. Rev. Lett._ 2008, _100_ , 020603
  * Valsson and Parrinello 2015 Valsson, O.; Parrinello, M. Well-Tempered Variational Approach to Enhanced Sampling. _J. Chem. Theory Comput._ 2015, _11_ , 1996–2002
  * Tiwary and Parrinello 2015 Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. _J. Phys. Chem. B_ 2015, _119_ , 736–742
  * Bal and Neyts 2016 Bal, K. M.; Neyts, E. C. Direct observation of realistic-temperature fuel combustion mechanisms in atomistic simulations. _Chem. Sci._ 2016, _7_ , 5280–5286
  * Fu and Pfaendtner 2018 Fu, C. D.; Pfaendtner, J. Lifting the Curse of Dimensionality on Enhanced Sampling of Reaction Networks with Parallel Bias Metadynamics. _J. Chem. Theory Comput._ 2018, _14_ , 2516–2525
  * Giberti et al. 2021 Giberti, F.; Tribello, G. A.; Ceriotti, M. Global Free-Energy Landscapes as a Smoothly Joined Collection of Local Maps. _J. Chem. Theory. Comput._ 2021, _17_ , 3292–3308
  * Jarzynski 2011 Jarzynski, C. Equalities and Inequalities: Irreversibility and the Second Law of Thermodynamics at the Nanoscale. _Annu. Rev. Condens. Matter Phys._ 2011, _2_ , 329–351
  * Crooks 1999 Crooks, G. E. Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences. _Phys. Rev. E_ 1999, _60_ , 2721–2726


# A Reduced-Complexity Maximum-Likelihood Detection with a sub-optimal BER
Requirement

Sharan Mourya
Indian Institute of Technology, Kharagpur
Hyderabad, India
<EMAIL_ADDRESS>Amit Kumar Dutta  Indian Institute of Technology,
Kharagpur
Kharagpur, India
<EMAIL_ADDRESS>
###### Abstract

Maximum likelihood (ML) detection is an optimal signal detection scheme, which
is often difficult to implement due to its high computational complexity,
especially in a multiple-input multiple-output (MIMO) scenario. In a system
with $N_{t}$ transmit antennas employing $M$-ary modulation, the ML-MIMO
detector requires $M^{N_{t}}$ cost function (CF) evaluations followed by a
search operation for detecting the symbol with the minimum CF value. However,
a practical system needs the bit-error ratio (BER) to be application-dependent
which could be sub-optimal. This implies that it may not be necessary to have
the minimal CF solution all the time. Rather it is desirable to search for a
solution that meets the required sub-optimal BER. In this work, we propose a
new detector design for a SISO/MIMO system by obtaining the relation between
BER and CF which also improves the computational complexity of the ML detector
for a sub-optimal BER.

Index Terms: Maximum Likelihood (ML) Detection, Multiple-input Multiple-output
(MIMO)

##  I Introduction

For any modern communication system, a maximum likelihood (ML) detector is
preferred due to its best data recovery performance[1] but it is
computationally very expensive to implement. This issue of computational
complexity in conventional digital machines is problematic as it increases the
latency with the increase of signal constellation size and the number of
transmitter antenna (TA). Magnificent tele-traffic growth in the last decade
has pushed the computational complexity and latency of a base station (BS) to
an alarming level. Diversified computational requirements and massive growth
in the connected user equipment (UE) in modern cellular standards make the
situation worse. From the physical layer point of view, higher complexity
operations include signal detection, parameter estimation, and various other
error corrections at the receiver. To address this, several low complexity
detection schemes were introduced in the recent years.

In [2], complexity of ML detection is reduced using Sensitive Bits (SB). [3]
proposes a pre-decoder guided local exhaustive search mechanism for V-BLAST
[4]. In [5] a similar reduction in ML complexity is achieved in M-PSK
modulated systems. Complexity reduction for a spatially modulated system is
proposed in [6]. A more generic approach with wide range of applicability is
considered in [7] where performance is traded off with computational
complexity. But this approach doesn’t address the search algorithm, it rather
deals with the complexity reduction of calculation of the cost function (CF).
In this work, we try to address the search complexity aspect of ML detection.

In the context of ML search, the extremum values from the database of
likelihood values dictate the positions of decision boundaries in the
constellation space. For example, for an additive white Gaussian noise (AWGN)
system, ML detector boils down to a minimum distance detector and decision
boundaries are the perpendicular bisectors of the lines joining any two
adjacent constellation points. Also, ML minimizes the BER for an equiprobable
source. If we choose any target BER other than the optimum value, then the
decision boundaries change accordingly and BER increases. We exploit this
trade-off in the subsequent sections in order to reduce the search complexity
in ML detection.

Contributions:

  1. 1.

In this work, we propose an ML detector which looks for a sub-optimal BER
while reducing the search complexity for ML detection.

  2. 2.

The proposed detection is first applied on a single input single output (SISO)
scenario and then extended to a multi-antenna case.

##  II System Model

###  II-A SISO

For a SISO system, a typical data model is given as

$\displaystyle y=hs+w,$ (1)

where $h,s,w$ are channel, transmitted data and AWGN, respectively. We assume
that $s$ is from the constellation space $S$ with cardinality $M$. For ML
detection, the likelihood function (LF) corresponding to any $s=s_{i}$ is
given as

$\displaystyle
LF(s_{i})=\displaystyle\frac{1}{\sqrt{2\pi\sigma^{2}}}\displaystyle\exp\left[-\frac{1}{2}\left(\frac{y-s_{i}}{\sigma}\right)^{2}\right].$
(2)

ML algorithm maximizes (2), which is the same as minimizing $(y-s_{i})^{2}$.
So we make the following search space $\Theta$ for the ML algorithm

$\displaystyle\Theta=\\{d_{0},d_{1},....,d_{M-1}\\},$ (3)

where $d_{i}=|y-s_{i}|$ for $s_{i}\in S$, i.e. $d_{i}$ is the distance between
the received symbol and $i^{th}$ constellation point. For ML detection,
finding the optimum value is actually finding the constellation point nearest
to the received symbol. Therefore, the search complexity depends on the size
of $\Theta$. For a sub-optimal BER we search for a value from a desired data
set other than $\Theta$. This data set, say $\Theta_{d}\subseteq\Theta$
contains all the values that ensure that BER is less than the maximum
allowable BER. With this motivation, we first start with a binary phase-shift
keying (BPSK) system.

####  II-A1 BPSK

Figure 1: Decision boundaries for BPSK

Choosing the minimum value from $\Theta$ will result in decision boundaries
being the perpendicular bisectors. So, when we choose a sub-optimal value, it
is logical to move the decision boundary from the midpoint. This is similar to
a Neyman Pearson (NP) criterion, where the boundary is shifted. However, the
NP criterion is best suited for binary hypothesis, not for multiple ones. As
this will lead to the determination of the decision boundary value for various
constellation points and may lead to multi-objective optimization problem. To
alleviate the issue, we propose to have symmetric decision regions. This leads
to the concept of a ”NULL” region. Hence, if the decision boundary moves in
one direction by some amount, we will add a mirror image of it on the opposite
side in order to preserve the symmetry. With this motivation, we propose the
following detector design for BPSK.

Let us choose the decision boundaries to be at a distance of $\beta$ from
$s_{0}$ and $s_{1}$ as shown in Fig. 1. Therefore, we get two decision regions
($R_{1}$ and $R_{3}$) each at a distance of $\beta$ from each constellation
point. In this case, we have an overlapping region $R_{2}$ which creates
ambiguity for detection. So, we propose the following detection criteria: If
$y$ is the received symbol and $\\{s_{0},s_{1}\\}$ form the constellation,
then

  1. 1.

If $y$ lies in the region $R_{1}$, then $\hat{s}=s_{0}$.

  2. 2.

If $y$ lies in the region $R_{3}$, then $\hat{s}=s_{1}$.

  3. 3.

If $y$ lies in the region $R_{2}$, then decide on $\hat{s}=s_{0}$ or
$\hat{s}=s_{1}$ with probability $\frac{1}{2}$.

Let us now calculate the average probability of error for this new detector.

Figure 2: Probability of error calculation for BPSK

For $s_{1}$, error occurs when $y$ falls in the region $R_{1}$ and error might
occur with probability $\frac{1}{2}$, when it falls in $R_{2}$. From Fig. 2,
let us assume that $A_{1}$ is the area of the curve in $R_{1}$ and $A_{2}$ is
the area of the curve in $R_{2}$. Then, the probability of error is given as

$\displaystyle P(e)$ $\displaystyle=A_{1}+\frac{1}{2}A_{2}$
$\displaystyle=\frac{1}{2}A_{1}+\frac{1}{2}(A_{1}+A_{2})$
$\displaystyle=\frac{1}{4}erfc\bigg{(}\frac{|h|(d_{min}-\beta)}{\sqrt{N_{0}}}\bigg{)}+\frac{1}{4}erfc\bigg{(}\frac{|h|\beta}{\sqrt{N_{0}}}\bigg{)}.$
(4)

For an ML detector, average probability of error is given by substituting
$\beta=d_{min}/2$ as

$\displaystyle
P_{min}=\frac{1}{2}erfc\bigg{(}\frac{|h|d_{min}}{2\sqrt{N_{0}}}\bigg{)}.$ (5)

We call it $P_{min}$ because it is the least value that can be achieved for
any value of $\beta$ i.e., $P_{e}\geq P_{min}$. So by choosing a different
value from $\Theta$ instead of the minimum value, we increase the probability
of error. Hence, it is possible to fix a BER (application dependent) that is
greater than the minimum BER and obtain the value of $\beta$ from equation
(4). We obtain a relationship between $\beta$ and $\Theta_{d}$ in Section-
III.B. Also, as we are addressing the issue of complexity rather than
performance, we will drop the channel coefficients from hereon.

####  II-A2 PAM

Let us follow the same procedure for obtaining the average probability of
error of pulse amplitude modulation (PAM). We choose 4-PAM for simplicity. The
received signal distributions with the proposed decision regions are shown in
Fig. 3.

Figure 3: Decision boundaries for 4-PAM

$P(e/S_{0})=A_{1}+\frac{1}{2}A_{2}=\frac{1}{2}A_{1}+\frac{1}{2}(A_{1}+A_{2})$
$P(e/S_{1})=2A_{1}+A_{2}=A_{1}+(A_{1}+A_{2})$
$P(e/S_{2})=2A_{1}+A_{2}=A_{1}+(A_{1}+A_{2})$
$P(e/S_{3})=A_{1}+\frac{1}{2}A_{2}=\frac{1}{2}A_{1}+\frac{1}{2}(A_{1}+A_{2})$

$where,$

$A_{1}=\frac{1}{2}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}$
$A_{1}+A_{2}=\frac{1}{2}erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}$
$\because P(e)=\frac{1}{4}(P(e/S_{0})+P(e/S_{1})+P(e/S_{2})+P(e/S_{3}))$
$\Rightarrow P(e)=\frac{3}{4}(A_{1}+(A_{1}+A_{2}))$
$P(e)=\frac{3}{8}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}+\frac{3}{8}erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}.$
(6)

####  II-A3 16-QAM case

Similarly, let’s derive the expression of the average probability of error for
quadrature amplitude modulation (QAM). We derive it first for 16-QAM for
simplicity as follows,

$P(e)=\frac{3}{4}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}+\frac{3}{4}erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}.$
(7)

The detailed derivation is given in Appendix-B.

####  II-A4 M-QAM case

We will generalize the expression for $M$-QAM as,

$\displaystyle
P(e)=\bigg{(}1-\frac{1}{\sqrt{M}}\bigg{)}\bigg{[}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}+erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}.$
(8)

The detail derivation is given in Appendix-C.

###  II-B UNION BOUND

Before extending to MIMO, we derive the union bound of the proposed detector
for a SISO link.

####  II-B1 SISO with M-QAM case

To derive the union bound for the proposed symmetric based decision regions,
we adopt the nearest neighbour approximation. We need to categorize the
constellation depending on the number of neighbours the symbol has as per the
following table.

Symbols | Cardinality |  $\\#$ of nearest neighbours
---|---|---
Edges | 4 | 2
Interior | $(\sqrt{M}-2)^{2}$ | 4
Other | $4(\sqrt{M}-2)$ | 3

where symbols are categorized into edges, interior points and others.
Therefore, the probability of error for $M$-QAM can be bounded as

$\displaystyle
P(e)\leq\bigg{\\{}\frac{1}{2}Q\bigg{(}\frac{\beta}{\sqrt{2N_{0}}}\bigg{)}+\frac{1}{2}Q\bigg{(}\frac{d_{min}-\beta}{\sqrt{2N_{0}}}\bigg{)}\bigg{\\}}\frac{1}{M}\sum_{i=0}^{M-1}n_{i},$
(9)

where,

$\frac{1}{M}\sum_{i=0}^{M-1}n_{i}=\frac{4\times
2+4\times(\sqrt{M}-2)^{2}+3\times 4(\sqrt{M}-2)}{M}$
$\Rightarrow\frac{1}{M}\sum_{i=0}^{M-1}n_{i}=4\bigg{(}1-\frac{1}{\sqrt{M}}\bigg{)}$
$\Rightarrow P(e)\approx
4\bigg{(}1-\frac{1}{\sqrt{M}}\bigg{)}\bigg{\\{}\frac{1}{2}Q\bigg{(}\frac{\beta}{\sqrt{2N_{0}}}\bigg{)}+\frac{1}{2}Q\bigg{(}\frac{d_{min}-\beta}{\sqrt{2N_{0}}}\bigg{)}\bigg{\\}}$

Notice that this is exactly similar to Eq (8) and it can be further simplified
to

$\displaystyle P(e)\leq
4\bigg{(}1-\frac{1}{\sqrt{M}}\bigg{)}\bigg{\\{}\frac{1}{2}e^{-\frac{\beta^{2}}{4N_{0}}}+\frac{1}{2}e^{-\frac{(d_{min}-\beta)^{2}}{4N_{0}}}\bigg{\\}}.$
(10)

####  II-B2 MIMO

Consider a $N_{r}\times N_{t}$ MIMO system (i.e., with $N_{t}$ transmitter
antennas and $N_{r}$ receiver antennas) with an AWGN channel. Let $M_{i}-QAM$
be the modulation of $i^{th}$ transmitting antenna then the union bound
expression is given as (from [8]),

$\displaystyle
P(e)\leq\displaystyle\frac{1}{\prod_{i=1}^{N_{t}}M_{i}}\sum_{x}\sum_{\hat{x}\neq
x}P_{2}(x,\hat{x}).$ (11)

Here, the first summation is over all possible input vectors (i.e.,
$\prod_{i=1}^{N_{t}}M_{i}$) and $P_{2}(x,\hat{x})$ is the pair-wise error
probability given as follows

$\displaystyle P_{2}(x,\hat{x})$
$\displaystyle=E\bigg{\\{}Q\bigg{(}\frac{|x-\hat{x}|^{2}}{2N_{0}}\bigg{)}\bigg{\\}}$
$\displaystyle=\bigg{[}\frac{1}{2}(1-\mu)\bigg{]}^{N_{r}}\sum_{k=0}^{N_{r}-1}{N_{r}-1+k\choose
k}\bigg{[}\frac{1}{2}(1+\mu)\bigg{]}^{k},$ (12)

where $\mu=\displaystyle\sqrt{\frac{\bar{\gamma_{c}}}{1+\bar{\gamma_{c}}}}$
and $\bar{\gamma_{c}}=\displaystyle\frac{|x-\hat{x}|^{2}}{2N_{0}}$. Also, from
[8], if $\bar{\gamma_{c}}\gg 1$ then $P_{2}(x,\hat{x})$ would become

$P_{2}(x,\hat{x})=\bigg{(}\frac{1}{4\bar{\gamma_{c}}}\bigg{)}^{N_{r}}{2N_{r}-1\choose
N_{r}}.$ (13)

Then, the simplified union bound expression is

$P(e)\leq\frac{(2N_{0})^{N_{r}}{2N_{r}-1\choose
N_{r}}}{\prod_{i=1}^{N_{t}}M_{i}}\sum_{x}\sum_{\hat{x}\neq
x}\bigg{(}\frac{1}{4|x-\hat{x}|^{2}}\bigg{)}^{N_{r}}.$ (14)

##  III IMPLEMENTATION

There are a few problems in this detector to deal with.
1) Obtaining $\beta$ from $P(e)$.
2) Translation of the proposed detector to a search algorithm.

###  III-A Proposed Solution for $\beta$

From Eq (7) it is clear that it is not possible to write a closed-form
expression for $\beta$ so we resort to numerical methods. We use the Newton-
Raphson method to find a solution to Eq (7). Let $P$ be the required BER,
$P_{min}$ be the BER achieved by the ML detector and $k=log_{2}(M)$.

$g(\beta)=\frac{3}{4k}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}+\frac{3}{4k}erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}-P.$
(15)

Then, the derivative of this function is

$g^{{}^{\prime}}(\beta)=\frac{3}{2k\sqrt{\pi}}e^{-\frac{(d_{min}-x)^{2}}{N_{0}}}-\frac{3}{2k\sqrt{\pi}}e^{-\frac{x^{2}}{N_{0}}}.$
(16)

Using the above two equations we start the iteration with $\beta_{0}=0$ and
the following update equation.

$\beta_{n+1}=\beta_{n}-\frac{g(\beta_{n})}{g^{{}^{\prime}}(\beta_{n})}.$ (17)

Comment: For a MIMO system, calculation remains the same except that
$g(\beta)$ becomes the Union Bound of Probability of Error we calculated in
equation (14).

$g(\beta)=\frac{(2N_{0})^{N_{r}}{2N_{r}-1\choose
N_{r}}}{\prod_{i=1}^{N_{t}}M_{i}}\sum_{i}\sum_{j\neq
i}\bigg{[}\bigg{(}\frac{1}{4|d_{ij}-\beta|^{2}}\bigg{)}^{N_{r}}$
$+\bigg{(}\frac{1}{4\beta^{2}}\bigg{)}^{N_{r}}\bigg{]}-P.$
$g^{{}^{\prime}}(\beta)=\frac{(2N_{0})^{N_{r}}{2N_{r}-1\choose
N_{r}}}{\prod_{i=1}^{N_{t}}M_{i}}8N_{r}\sum_{i}\sum_{j\neq
i}\bigg{[}-\beta\bigg{(}\frac{1}{4\beta^{2}}\bigg{)}^{N_{r}+1}$
$+|d_{min}-\beta|\bigg{(}\frac{1}{4|d_{ij}-\beta|^{2}}\bigg{)}^{N_{r}+1}\bigg{]}.$

###  III-B Proposed simplification of Detector

Implementing this detector requires a uniform distribution which adds extra
circuitry to the detector. Fortunately, we can implement this detector without
ever using a uniform distribution. Consider the BPSK again as in Fig. 4a,

(a) Decision boundaries for BPSK

(b) Equivalent receiver with just one decision boundary

Figure 4: Proposed simplification of the receiver

Probability of error expressions for this scenario are

$\displaystyle P_{1}(e)$ $\displaystyle=A_{1}+\frac{1}{2}A_{2}$
$\displaystyle=\frac{1}{2}A_{1}+\frac{1}{2}(A_{1}+A_{2}).$ (18)

Now, consider the following scenario with only one decision boundary as in
Fig. 4b. Distance between $s_{0}$ and $s_{1}$ is $d_{min}$ and distance
between $s_{1}$ and $D_{1}$ is $\beta$. Probability of error in this case is
given as,

$\displaystyle P_{2}(e)=\frac{1}{2}A_{3}+\frac{1}{2}A_{4}.$ (19)

As $A_{3}=A_{1}$ and $A_{4}=A_{1}+A_{2}$, we can write

$\displaystyle P_{2}(e)=\frac{1}{2}A_{1}+\frac{1}{2}(A_{1}+A_{2}).$ (20)

The motivation for writing the equations in this form is clear by observing
that $P_{1}(e)=P_{2}(e)$. Now, detection criteria for this detector can be
written as

$\displaystyle\hat{s}$
$\displaystyle=s_{1}\quad\text{if}\quad|y-s_{1}|\leq\beta,$
$\displaystyle=s_{0}\quad\text{elsewhere.}$ (21)

Given this criteria, we can now establish the relation between $\beta$ and
$\Theta_{d}$. Firstly, we fix the BER at $P$ and obtain the value of $\beta$
as described in the last subsection and to practically achieve that BER we
collect all the $d_{i}$’s such that $d_{i}\leq\beta$. These $d_{i}$’s form the
set $\Theta_{d}$. In other words, choosing any $d_{i}\in\Theta_{d}$ from
$\Theta$ will achieve the required BER.

##  IV Numerical Results

###  IV-A Probability of Hitting

For an ML detector, probability of hitting a solution in the search algorithm
for a MIMO system with $N_{t}$ transmit antennas each employed with M-QAM
modulation is:

$P_{ML}(H)=\frac{1}{M^{N_{t}}}$

For the proposed detector, we changed the search criterion to Eq (21) which
directly reduces the search space. As we are considering an AWGN environment
with mean 0 and variance $N_{0}$, all the search elements in $\Theta_{d}$ are
Gaussian distributed. Thus, probability of hitting becomes:

$P_{PD}(H)=\frac{M^{N_{t}}\times P(|y-s|\leq\beta)}{M^{N_{t}}}$
$P_{PD}(H)=\frac{M^{N_{t}}\times(1-Q(N_{0}\beta))}{M^{N_{t}}}$ $\Rightarrow
P_{PD}(H)=1-Q(N_{0}\beta)$

It can be see that, the proposed detector removes the exponential dependency
of probability of hitting on the transmitter and constellation sizes. Now it’s
just a function of $\beta$ that significantly helps in scaling the approach
for a massive MIMO system.

###  IV-B Complexity Analysis

To demonstrate how the proposed receiver compares to the ML receiver, we
performed the following simulation under the conditions.
$\bullet$ Modulation: 16-QAM with $2\times 2$ MIMO.
$\bullet$ Required BER: $P(\gamma)=2\times P_{min}(\gamma)$, where
$P_{min}(\gamma)$ is the BER obtained by ML receiver at the SNR ”$\gamma$”.
We assume perfect channel estimation and equalization here.

Procedure:
1) Obtain $\beta(\gamma_{0})$ at any $\gamma_{0}$ using Newton-Raphson.
2) Deploy linear search over the set $\Theta$.
3) Obtain the number of evaluations the search algorithm takes to detect every
symbol and average over all the symbols (10 million in this case).
4) Normalize the result to 1 by dividing it with $16^{2}$.
5) Repeat the same with the next SNR value.

Figure 5: Query Complexity vs SNR

ML receiver always takes the same number of evaluations, as it needs to find
the least value from the set, which is equal to the size of the search space.
Whereas the proposed receiver takes as many evaluations before it hits a value
that is less than $\beta(\gamma)$ (i.e., a value from the set $\Theta_{d}$).
This probability of hitting a solution increases with increasing SNR as the
$\beta(\gamma)$ value also increases with SNR provided $P(\gamma)$ is kept
constant. Hence, the number of evaluations required by the proposed receiver
decreases with increasing SNR as evident from the plot in Fig. 5.

###  IV-C BER vs SNR

We improved the computational complexity but we also do not want to trade this
off with performance. For this, we simulated the BER performance of the
proposed receiver on the same system as shown in Fig. 6.

Figure 6: BER vs SNR

##  V Conclusion

In this work, we have proposed a modified ML algorithm with application
dependent BER. Choosing a sub-optimal BER moves the decision boundaries such
that it increases the number of solutions for the cost function which in turn
reduces the size of search space for ML detection. We proposed a detector
design to exploit this property while eliminating the NULL regions. Then the
proposed detector is also extended to a MIMO scenario deriving the union bound
expression for probability of error calculation. This detector achieved a
significant computational complexity advantage depending on the target BER.
Simulation results also suggests the same.

## APPENDIX

### V-A

Let us consider the detector for 4-PAM.

Figure 7: 4-PAM constellation with new boundaries

We will calculate the average probability of error for this detector and
compare it with the proposed detector.

$P(e)=\frac{1}{4}(A_{1}+(A_{1}+A_{2})+(A_{1}+A_{2})+A_{2})$
$P(e)=\frac{3}{4}(A_{1}+A_{2})$

where,

$A_{1}=\frac{1}{2}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}$
$A_{2}=\frac{1}{2}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}$

Notice that $A_{2}$ is different from $A_{2}$ in Section II.

$P(e)=\frac{3}{8}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}+\frac{3}{8}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N}_{0}}\bigg{)}.$
(22)

### V-B

Consider 16-QAM as shown in the figure.

Figure 8: 16 QAM constellation with new boundaries

Consider the exterior points $S_{0}$, $S_{3}$, $S_{12}$, $S_{15}$. Probability
of detection and probability of error of these points are given as

$P(d/S_{0,15})=\bigg{[}1-\frac{1}{2}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}\bigg{[}1-\frac{1}{2}erfc\bigg{(}$
$\frac{\beta}{\sqrt{N}_{0}}\bigg{)}\bigg{]}$ $P(e/S_{0,15})=1-P(d/S_{0,15})$
$P(e/S_{3})=1-\bigg{[}1-\frac{1}{2}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}\bigg{]}^{2}$
$P(e/S_{12})=1-\bigg{[}1-\frac{1}{2}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N}_{0}}\bigg{)}\bigg{]}^{2}$

Similarly for the points
$S_{1}$,$S_{2}$,$S_{4}$,$S_{8}$,$S_{7}$,$S_{11}$,$S_{13}$,$S_{14}$ (neither
interior nor exterior points) probability of detection and probability of
error of these points are given by,

$\rightarrow
P(d/S_{1,2,7,11})=\bigg{[}1-\frac{1}{2}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}\bigg{]}\times$
$\bigg{[}1-\frac{1}{2}\bigg{[}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}+erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}\bigg{]}$
$P(e/S_{1,2,7,11})=1-P(d/S_{1,2,7,11})$

$\rightarrow
P(d/S_{4,8,13,14})=\bigg{[}1-\frac{1}{2}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N}_{0}}\bigg{)}\bigg{]}\times$
$\bigg{[}1-\frac{1}{2}\bigg{[}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}+erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}\bigg{]}$
$P(e/S_{4,8,3,14})=1-P(d/S_{4,8,13,14})$

For the interior points $S_{5}$,$S_{6}$,$S_{9}$,$S_{10}$, probability of
detection and probability of error of these points are given by,

$\rightarrow
P(d/S_{5,6,9,10})=\bigg{[}1-\frac{1}{2}\bigg{[}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}+$
$erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}\bigg{]}^{2}$
$P(e/S_{5,6,9,10})=1-P(d/S_{5,6,9,10})$

Total probability of error is given by,

$P(e)=\frac{1}{16}\bigg{[}\big{(}P(e/S_{3})+P(e/S_{13})+2P(e/S_{0})\big{)}+$
$\big{(}4P(e/S_{5})\big{)}+\big{(}4P(e/S_{1})+4P(e/S_{4})\big{)}\bigg{]}$

Substituting the values and simplifying it by neglecting the second order
terms gives,

$P(e)=\frac{3}{4}erfc\bigg{(}\frac{\beta}{\sqrt{N}_{0}}\bigg{)}+\frac{3}{4}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N}_{0}}\bigg{)}$
(23)

### V-C

We will derive the general expression for probability of error of an M-QAM.
Generalizing this to M-QAM we get $4$ exterior points, $(\sqrt{M}-2)^{2}$
interior points and $4(\sqrt{M}-2)$ points that are neither exterior nor
interior.

$P(e)=\frac{4}{M}P(e/S_{e})+\frac{(\sqrt{M}-2)^{2}}{M}P(e/S_{i})+$
$\frac{4(\sqrt{M}-2)}{M}P(e/S_{n})$ $\rightarrow
P(e)=\frac{1}{M}\bigg{[}\big{(}P(e/S_{3})+P(e/S_{13})+2P(e/S_{0})\big{)}$
$+\big{(}(\sqrt{M}-2)^{2}P(e/S_{5})\big{)}+$
$\big{(}2(\sqrt{M}-2)P(e/S_{1})+2(\sqrt{M}-2)P(e/S_{4})\big{)}\bigg{]}$

Simplification gives

$P(e)=\bigg{(}1-\frac{1}{\sqrt{M}}\bigg{)}\bigg{[}erfc\bigg{(}\frac{d_{min}-\beta}{\sqrt{N_{0}}}\bigg{)}+erfc\bigg{(}\frac{\beta}{\sqrt{N_{0}}}\bigg{)}\bigg{]}$

## References

  * [1] S. M. Kay., Fundamentals of Statistical Signal Processing: Estimation Theory.  Prentice Hall, 1993.
  * [2] J. Li, K. Letaief, and Z. Cao, “A reduced-complexity maximum-likelihood method for multiuser detection,” IEEE Transactions on Communications, vol. 52, no. 2, pp. 289–295, 2004.
  * [3] L. He and H. Ge, “Reduced complexity maximum likelihood detection for v-blast systems,” in IEEE Military Communications Conference, 2003. MILCOM 2003., vol. 2, pp. 1386–1391 Vol.2, 2003.
  * [4] P. Wolniansky, G. Foschini, G. Golden, and R. Valenzuela, “V-blast: an architecture for realizing very high data rates over the rich-scattering wireless channel,” in 1998 URSI International Symposium on Signals, Systems, and Electronics. Conference Proceedings (Cat. No.98EX167), pp. 295–300, 1998.
  * [5] H. Men and M. Jin, “A low-complexity ml detection algorithm for spatial modulation systems with $m$psk constellation,” IEEE Communications Letters, vol. 18, no. 8, pp. 1375–1378, 2014.
  * [6] R. Rajashekar, K. Hari, and L. Hanzo, “Reduced-complexity ml detection and capacity-optimized training for spatial modulation systems,” IEEE Transactions on Communications, vol. 62, no. 1, pp. 112–125, 2014.
  * [7] J.-S. Kim, S.-H. Moon, and I. Lee, “A new reduced complexity ml detection scheme for mimo systems,” IEEE Transactions on Communications, vol. 58, no. 4, pp. 1302–1310, 2010.
  * [8] G. Lindell, “Some exact union bound results for maximum likelihood detection in mimo systems,” Proceedings. International Symposium on Information Theory, 2005. ISIT 2005., pp. 2223–2227, 2005.