# Using Metadynamics to Reveal Extractant Conformational Free Energy
Landscapes

Xiaoyu Wang  Chemical Sciences and Engineering Division, Argonne National
Laboratory, 9700 S Cass Ave, Lemont, IL 60439<EMAIL_ADDRESS>Michael J.
Servis  Chemical Sciences and Engineering Division, Argonne National
Laboratory, 9700 S Cass Ave, Lemont, IL 60439<EMAIL_ADDRESS>
###### Abstract

Understanding the impact of extractant functionalization on metal binding
energetics in liquid-liquid extraction (LLE) is essential to guide the
development of better separations processes. Traditionally, computational
extractant design uses electronic structure calculations on metal-ligand
clusters to determine the metal binding energy of the lowest energy state.
Although highly accurate, this approach does not account for all the relevant
physics encountered under experimental conditions. Such methodologies often
neglect entropic contributions, such as temperature effects and ligand
flexibility, in addition to approximating solvent-extractant interactions with
implicit solvent models. In this study, we use classical molecular dynamics
simulations with an advanced sampling method, metadynamics, to map out
extractant molecule conformational free energies in the condensed phase. We
generate the complete conformational landscape in solution for a family of
bidentate malonamide-based extractants with different functionalizations of
the head group and the side chains. In particular, we show how such alkyl
functionalization reshapes the free energy landscape, affecting the free
energy penalty of organizing the extractant into the cis-like metal binding
conformation from the trans-like conformation of the free extractant in
solution. Specifically, functionalizing alkyl tails to the center of the head
group has a greater influence on increasing molecular rigidity and disfavoring
the binding conformation than functionalizing side chains. These findings are
consistent with trends in metal binding energetics based on experimentally
reported distribution ratios. We also consider a different bidentate
extractant molecule, CMPO, and show how the choice of solvent can further
reshape the conformational energetic landscape. This study demonstrates the
feasibility of using molecular dynamics simulations with advance sampling
techniques to investigate extractant conformational energetics in solution,
which, more broadly, will enable extractant design that accounts for entropic
effects and explicit solvation.

###### keywords:

American Chemical Society, LaTeX

LLE,MD,Metadynamics

##  1 Introduction

Liquid liquid extraction (LLE), also called solvent extraction, is a
commonplace hydrometallurgical process used to separate and recover critical
metals, such as rare earth, transuranic, and platinum group elements. 1 In
LLE, metal ions are transported from an aqueous feed phase, often acidic, to
an immiscible organic phase that typically contains an amphiphilic extractant
molecule, which binds the targeted metal ions, and a nonpolar solvent, also
called the diluent. LLE is a thermodynamic process driven by small free energy
differences in metal solvation between phases, in which large energetic terms
cancel out, resulting in an overall energy on the order of $-k_{B}T$ to ensure
reversibility for further reprocessing. 2 While an LLE process involves many
energetic terms that are often challenging to isolate and quantify, their net
contributions are reflected in the overall separations performance.
Specifically, the distribution ratio is the ratio of concentration of the
target metal ion in the extractant phase to the aqueous phase. While metal-
extractant binding is typically a strong, enthalpically favorable interaction,
for multidentate extractants, the reorganization of the extractant into its
binding conformation is often highly unfavorable. 3, 4 Thus, traditional
ligand design work has focused on utilizing state-of-the-art density
functional theory (DFT) to suggest candidate extractants by constraining the
flexibility of functional head groups. 5, 6, 7, 8 In these works, a handful of
molecular structures are investigated by geometry optimizations at 0 K
followed by single-point energy calculations. While capable of computing
highly accurate energies, such a strategy is insufficient to provide a full
picture of molecular energetics at temperature in the solution phase, as the
computational cost of DFT calculations prohibits consideration of the ensemble
of structures potentially relevant to metal extraction under experimentally
realistic conditions. Moreover, since geometry optimization is carried out at
0 K, neither the flexibility of molecular structure nor the contribution of
entropy are considered. In addition, solvents are typically treated implicitly
by mean-filed theories, which neglect that, at the molecular level, solvents
can interact with extractants in more complex ways.

Herein, we present classical molecular dynamics (CMD) simulations and profile
the full energetic landscape of extractant conformations, with solvent
molecules being modeled explicitly. To enhance the sampling of all extractant
conformations, especially for unfavorable high-energetic ones, a free energy
estimation technique, metadynamics (MTD) 9, 10, 11, is carried out by adding
repulsive bias potentials as functions of observable collective variables.
Unlike regular unbiased CMD, these time-evolving bias potentials discourage
the system from visiting low-energetic configurations, which provides a way of
mapping out the whole free-energy surface defined by the collective variables
with reasonable computational cost. In contrast to DFT techniques, which
necessitate the simplification of solvents to avoid expensive calculations,
CMD uses semiempirical force fields to explicitly account for solvent
molecules. While empirical CMD potentials are less accurate than those
calculated from DFT, we are able to benchmark the force field performance
against DFT, and then explore the full range of conformations present in the
solution phase. By using DFT to identify the binding geometry of the
extractant molecule in the presence of the metal, we can then compute the free
energy cost of the molecule organizing into that conformation. In a practical
system, this energetic cost competes with the favorable metal-extractant
binding. Here, by isolating this conformation energy in the absence of the
metal, but at temperature and in the presence of explicit solvent, we can
isolate how changes in the extractant structure reshape the conformational
energetic landscape. This will quantify how the extractant structure favors or
disfavors the binding conformation compared to the lowest energy conformation
that is typically assumed by the nonbinding extractant. Overall, this will
provide a complete picture of the inherent conformational energetics of an
extractant molecule, allowing us to connect changes to this landscape to
separations performance.

In this work, we consider two major types of extractants: malonamides and
carbamoylmethylphosphine oxides (CMPO). For the malonamide type, we choose
$N,N,N^{\prime},N^{\prime}$-tetrahexyl malonamide (THMA),
$N,N,N^{\prime},N^{\prime}$-tetrahexyl-2-methyl malonamide (MeTHMA),
$N,N,N^{\prime},N^{\prime}$-tetrahexyl-2,2-dimethyl malonamide (DiMeTHMA),
$N,N^{\prime}$-dimethyl-$N,N^{\prime}$-dibutyl-2-tetradecyl malonamide
(DMDBTDMA), and 3,9-diaza-3,9-dioctylbicyclo[4.4.0]decane-2,10-dione (BMA),
which have shown significantly different separation performances with slight
modifications in molecular structures. 12, 13, 14, 15 For the type of CMPO, we
select the most famous derivative, n-octyl(phenyl)-$N,N^{\prime}$-diisobutyl
CMPO (we will simply refer to it as CMPO in the following context), which has
been implemented in the TRansUranic EXtraction (TRUEX) process. 16 Figure 1
shows the molecular structures of these extractants covered in this work. By
constructing the two adjacent torsion angles in the extractants’ head group as
collective variables (shown as $\psi$ and $\phi$ in Figure 1), we are able to
sample and generate the full conformational free-energy landscapes by varying
$\psi$ and $\phi$. Such free energy landscapes are able to represent full
energetic information of these multidentate extractant conformations,
including the free unbound trans, transitional gauche and bound cis
conformations. The dramatically different energetic landscapes across
different extractant functionalizations provides different perspectives on how
LLE separations can be influenced by the molecular structure and
conformational energetics.

Figure 1: Schematic molecular structures of extractants considered in this
study.

##  2 Simulation methodology

###  2.1 DFT calculations

For all DFT calculations, we use a complex consisting of one Eu(III), with one
bidentate malonamide or CMPO and three bindentate nitrate anions to fill the
coordination vacancies and neutralize the complex. The Eu(III) cation has 7 –
9 coordination sites, and we choose 8 as the coordination number in all DFT
calculations for simplicity. Note that the coordination chemistry of
lanthanides and extraction stoichemetry are beyond the scope of this study;
the metal-extractant complexes are selected only to provide a preliminary
understanding of the extractant geometry as it binds with the metal ion,
without deliberately considering variation in stoichimetry of the complex or
nitrate denticity. The geometry of the Eu(III)-extractant complex is optimized
with M06 17 hybrid functional theory. The Stuttgart/Dresden relativistic
small-core (RSC) effective core potential (ECP) basis set, obtained from Basis
Set Exchange 18, is used for Eu(III), while 6-31G(d) is used for other
elements in the complex. All DFT calculations are performed in Gaussian 16 19.

###  2.2 CMD and MTD

CMD uses semi-empirical force fields to describe the interaction between
bonded and nonbonded atoms in the simulation system. By solving the Newtonian
equation of motion, atomic positions are propagated with time, from which
ensemble quantities are averaged and compared with experiments. One of the
challenges of CMD is that the simulated system is often “trapped” in local
minima on the free energy surface (FES), which prevents the CMD from exploring
the entire configurational space. The detailed underlying theory of how MTD
works is beyond the scope of this work; instead, we will only briefly
introduce the concept of MTD. MTD is a type of advanced sampling method that
operates by spawning history-dependent Gaussian bias potentials to the
original FES, which discourages the system from re-visiting low-energetic
configurations and eventually helps to “escape” local energy minima. The MTD
simulation will stop when the whole configurational space is explored, and the
negative of the summation of all spawned Gaussian bias potentials will be
removed from the original FES.

For all simulations, we model the extractants using the General AMBER Force
Field (GAFF2) 20 with partial charges using the Austin Model 1 with bond
charge correction (AM1-BCC) method. 21, 22 We also test the Optimized
Potentials for Liquid Simulations (OPLS) force field 23 (see below). The
parameter set for dodecane is taken from an optimized GAFF model to prevent
unphysical freezing at room temperature. 24 The Lorentz-Berthelot mixing rule
is used for all Lennard-Jones cross-term interactions. A trans extractant
molecule is packed with 150 dodecane molecules using PACKMOL 25. The
simulation boxe volumes are equilibrated for 1 ns in NPT ensembles with a
Nose-Hoover style thermostat and barostat 26, 27, 28 at 303 K and 1
atmosphere. Then, simulations are switched to NVT ensemble with Nose-Hoover
style thermostat at 303 K for 100 ns MTD sampling. All simulations are carried
out with a 1 fs timestep. For MTD simulations, every 1000 steps, a Gaussian
bias potential is added to the dihedral angles $\psi$ and $\phi$ defined in
Figure 1. The height of the potential is 0.5 kJ/mol, and the width is 0.1 rad
for both $\psi$ and $\phi$. We use LAMMPS 29 as the CMD engine, and MTD
simulations are performed with the open source, community-developed PLUMED
software package. 30

##  3 Results and discussions

###  3.1 Force field comparison

To evaluate whether the general purposed classical MD force fields are able to
capture the conformation energetics of the extractant molecules accurately, we
compare a single gas-phase molecule modeled using GAFF2 and OPLS force fields
to DFT calculations. We consider the CMPO extractant because the molecule is
not symmetric and has a phosphine oxide group which may be poorly defined in
these general-purposed classical force fields. Figure 2 shows the Ramachandran
plots for molecular mechanics and quantum mechanics. The Ramamchandran plot 31
is a 2D heat map showing the properties of interest in different combinations
of two adjacent dihedrals: here, we use it to show the potential energy
surfaces (PES) by varying the dihedral angles of O=P-C-C ($\phi$) and P-C-C=O
($\psi$) in CMPO (shown in Figure 1), which control the head group
conformations of the extractant. The PES for GAFF2 shows that, compared to DFT
calculations, the classical force field is able to capture all the fine
features, given that we have not performed any reparameterizations for the
CMPO molecule. Major deviations from DFT in GAFF2 are observed only in some
high-energy areas at extreme values of $\phi$ and $\psi$. Such deviations are
acceptable considering that classical force fields, which use harmonic bond
and angle potentials, often miss-treat intramolecular potentials when
molecules are severely distorted. Compared to GAFF2 and DFT, OPLS does not
accurately describe molecular mechanics: OPLS captures neither the shape nor
the peak height of many features of the PES. Therefore, we proceed with the
CMD simulations using GAFF2.

Figure 2: PESs of CMPO in gas-phase DFT calculations (left), GAFF2 (middle),
and OPLS (right). Axes are defined by the corresponding dihedral angles in
Figure 1. All energy contour maps are sharing one color-coding in the right
panel.

###  3.2 Effects of alkyl substitution on malonamide $\alpha$-carbon

To test the impact of functionalization on the head group of the extractant,
we investigate three malonamides where the substitutions of $\alpha$-carbon
are varied. The addition of alkyl functional groups to the $\alpha$-carbon
that joins the amide groups in the malonamide molecule is known to affect
separations performance. 13, 12 Specifically, when MeTHMA and THMA are
compared, the middle methyl group on the head group is known to result in
lower distribution ratios for both uranyl and Eu(III) extraction.13, 12 A more
significant drop in the distribution ratios for DiMeTHMA has been observed by
grafting the second methyl group, where nearly negligible extraction was
reported for both uranyl and Eu(III) metal ions. 13, 12 We expect that this
effect will be reflected in the conformational energetics of the molecule,
which we explore through DFT and CMD supplemented with MTD. In particular, we
hypothesize that such drops in distribution ratios can potentially be
explained by the increased free energy penalty for assuming the binding
conformation with sequential $\alpha$-carbon functionalization based on
comparing DFT calculations of the bound-state molecular geometry to the MTD-
determined conformational free energy landscape.

In order to identify the bound-state extractant conformation, we first carry
out DFT calculations in which an extractant molecule forms a complexant with
one Eu(III) and three nitrates. The molecular structures of these clusters are
shown in the top panels of Figure 3. As we discussed above, compared with
static DFT calculations, CMD can provide more details on the extractant
confomrational energetics, especially when temperature and solvent effect are
taken into consideration. The bottom panels in Figure 3 show Ramamchandran
plots, in which the conformational free energy landscapes are mapped from CMD
simulations using the MTD sampling technique. We note that, unlike the
potential energy surface calculated from single gas-phase configurations
described in the force field validation and shown in Figure 2, the MTD
simulations provide free energy surfaces. The optimized DFT structures of the
extractant bound with the Eu(III) provide the combinations of $\phi$ and
$\psi$ assumed to define the binding state, which we project onto the 2D
Ramamchandran free energy surface to show where the binding conformation falls
in the total energetic landscape of the extractant. We calculate one DFT
structure per malonamide; the multiple bound states in Figure 3 are inferred
from molecular symmetry. Note that our DFT calculations are only used to
determine $\phi$ and $\psi$ at the bound state conformation of a single
extractant, and we are not trying to replicate any particular experimental
stoichiometry of the metal-extractant complex. Nevertheless, all our DFT
optimized structures suggest that cis-malonamide binds with the metal ion in a
bidentate fashion, consistent with experimental studies. 32, 33, 34

With the binding conformation of each malonamide identified, we proceed with
CMD using MTD to determine how changes to the molecular structure affects the
cis conformation in the absence of the metal. From inspection of the binding
geometries and conformational free energy surface, the fundamental competition
between metal binding and molecular deformation is apparent. The enthalpically
favorable metal binding process favors the cis conformation, while the
intramolecular conformational energetics prefer gauche and, even more so,
trans states. Therefore, changes in the extractant that move the accessible
cis-like regions further from the true cis region near $\phi=0$ and $\psi=0$
presumably weaken the metal-extractant binding. Thus, when the conformational
energetic penalty for metal binding is added to the favorable binding process,
the overall free energy of extraction decreases.

Figure 3: Top: bound state structures optimized by DFT. Bottom: Ramachandran
plots of a) THMA, b) MeTHMA, c) DiMeTHMA in dodecane from MTD simulations.
Axes are defined by the corresponding dihedral angles in Figure 1. All energy
contour maps are sharing one color-coding given in panel c). Stars in each map
represent the DFT-optimized conformations shown at the top, when extractants
bound with an Eu(III) ion.

When a methyl group is grafted on the $\alpha$-carbon from THMA to MeTHMA,
free energy Ramachandran plots have shown that the flexibility of the
extractant decreases: low-energetic areas are much more confined for MeTHMA,
showing conformations that become inaccessible. The free energy of the
$\phi=0$ and $\psi=0$ region increases significantly for MeTHMA, making the
binding conformations less accessible, which potentially explains the
observation that the distribution ratio drops from THMA to MeTHMA. For THMA,
since there is additional symmetry around the $\alpha$-carbon, they do not
possess chirality and should also be symmetric across $y=x$ as well. For
MeTHMA, the middle alkyl group gives chirality to the molecule, so we do not
expect this symmetry. (In fact, this symmetry helps to demonstrate when the
MTD trajectory has adequately sampled the conformational landscape, as the
features and energy scales should be consistent across the applicable symmetry
planes.) The unfavorable cis region also shifts to the bottom right of the
Ramachandran plot. However, the binding conformation responds to the higher
energetic cost of the cis conformation and the asymmetry of the free energy
surface across $y=x$. The optimized bound state $\phi$ and $\psi$ identified
by DFT shifts to the upper left corner of the cis region for MeTHMA compared
to THMA, reflecting these changes to the underlying free energy surface. The
DFT optimized bound-state structures also show that although the cis
conformation is most relevant for both THMA and MeTHMA, these structures are
slightly different, where carbonyl groups are more parallel for MeTHMA.

By adding the second methyl group on the $\alpha$-carbon, going from MeTHMA to
DiMeTHMA, the most significant change in Figure 3 is that the free energy
landscape becomes symmetric across $y=x$ again. The DFT optimized bound state
$\phi$ and $\psi$ also shift back to positions more similar to THMA. Overall,
compared with the other malonamides, the low-energetic areas in the free-
energy landscape become even more confined. The central (cis) high-energy area
shifts back closer to $\phi=0$ and $\psi=0$, resulting in the disappearance of
a low-energy area that is easily accessible for the bound state. These effects
potentially explain why DiMeTHMA only shows negligible extraction performance
for Eu(III) and uranyl metal ions, where the high energetic barrier near the
cis conformation makes the binding state for DiMeTHMA much less energetically
accessible.

###  3.3 Effects of alkyl tail length on malonamides

The functionalization of longer alkyl chains provides better solubility in the
organic phase,35, 36 improving the metal loading capacity. However, the
potential effect of alkyl length on conformational flexibility is not
understood. To test the effect of side chain lengths, we first compare
MeDMDBMA, an analogue of MeTHMA but with shorter and different alkyl side
chains. The left and center panels in Figure 4 shows the free energy
landscapes of MeDMDBMA and MeTHMA, respectively. The free energy maps of these
two molecules are similar, with DFT-determined bound-state conformations at
similar locations on the Ramachandran plot. We also consider DMDBTDMA, an
extractant used in the DIAMEX (DIAMide Extraction) process for the recovery of
trivalent actinide ions,37 which has a much longer alkyl tail, replacing the
methyl group on the $\alpha$-carbon of MeDMDBMA. Noting that, compared to
THMA, slightly better separation performance has been reported with DMDBTDMA
for uranium, 38 it is worth mentioning that other interaction contributions,
such as the assembly of multiple extractants and the complexation with metal
ions, will also play important roles in driving separation. It is likely that
the additional sides chains of THMA introduce steric hindrance that disfavors
formation of extracted metal complexes, which we are not considering in this
study. Further simulation studies using MTD on metal-extractants solvates
could identify differences in steric effects between extractants with
different side chain lengths. Overall, the unfavroable regions for DMDBTDMA
have energies higher than those of MeTHMA and MeDMDBMA, indicating that the
length of the alkyl group on the $\alpha$-carbon has a greater influence on
increasing conformational energetics and lowering the head group’s flexibility
than on the side chain. While the addition of an alkyl tail length to the
molecule is necessary to improve its solubility, these results suggest that
different locations on the malonamide affect its conformation flexibility and,
presumably, its separation performance.

Figure 4: Top: bound state DMDBTDMA optimized by DFT; bottom: Ramachandran
plot of DMDBTDMA in dodecane from MTD simulations. X and Y axes are defined by
the dihedral angles in Figure 1. Stars represent the DFT-optimized
conformations shown at the top, when the extractant bound with an Eu(III) ion.

###  3.4 Bicyclic malonamide

A common strategy in ligand design to improve selectivity and separation
performance is to rigidify the backbone of the molecule, which preorganizes
the molecule into the preferred binding conformation. 5, 39 To design a better
malonamide for rare earth separation, Lumetta et al. synthesized a bicyclic
malonamide, finding 7 orders of magnitude improvement in the distribution
ratio of Eu(III).15 Such a significant increase in the distribution ratio was
achieved through restricting the head group to the binding conformation by
introducing a bicyclic carbon group. We use the MTD to map out the free energy
landscape of this bicyclic malonamide, with the Ramachandran plot shown in
Figure 5.

Figure 5: Top: bound state BMA optimized by DFT; bottom: Ramachandran plot of
BMA in dodecane from MTD simulations. X and Y axes are defined by the dihedral
angles in Figure 1. Stars represent the DFT-optimized conformations shown at
the top, when the extractant bound with an Eu(III) ion. We use the white color
to represent the region, where the CMD simulation cannot sample even with MTD.

The conformational free-energy landscape of the bicyclic malonamide is
dramatically different from that of all of the malonamides investigated above.
The trans conformation, typically where the global mininum resides, completely
disappears. Instead, the extrantant molecule is now restricted in the binding,
cis conformation, shown in blue colors in the top left quadrant. There is a
wide region (colored by white) where our CMD cannot explore this
configurational space even with the MTD, due to the high-energy barrier. The
bound-state structures calculated with DFT are much closer to the minima than
with the other malonamides. The energetic difference between the bound
structure and the minimum free energy conformation is around 15 kJ/mol, which
is significantly less than those of the malonamides investigated above,
although not sufficiently less to explain the seven orders of magnitude
increase in the Eu(III) distribution compared to THMA reported by Lumetta et
al..

###  3.5 CMPO: Asymmetric binding moieties

Compared with traditional PUREX (Plutonium Uranium Reduction EXtraction), CMPO
was designed to improve the functionality and separation of actinides in LLE.
40, 41, 16 Compared to the malonamides that we have covered above, CMPO not
only has a highly asymmetric molecular structure, but also has an entirely
different functional head group consisting of the phosphine oxide group in
addition to the carbonyl group. Here, we will use CMD with MTD to investigate
its conformational behavior and show its differences in comparison to that of
the malonamides. The left panel in Figure 6 shows the energetic landscape of
CMPO in dodecane, with its asymmetry readily apparent in the conformation free
energy surface. In comparison to the malonamides we investigated, the
conformational free energies of CMPO are generally much shallower, not only
for the bound conformation but also for other nearby conformations. Such a
lower energy landscape indicates that more conformations are easily accessible
at room temperature, which gives CMPO greater flexibility in solution. This
lower overall energy of the binding state of CMPO could explain its increased
performance compared with that of malonamides, in addition to the greater
basicity of the phosphine oxide binding site.

Figure 6: Top: bound state CMPO optimized by DFT; bottom: Ramachandran plot of
CMPO in dodecane and TBP from MTD simulations. X and Y axes are defined by the
dihedral angles in Figure 1. Stars represent the DFT-optimized conformations
shown at the top, when the extractant bound with an Eu(III) ion.

A primary benefit of studying conformational energetics with CMD using MTD is
the ability to include explicit solvation. To explore this, we also
investigate CMPO in a solution of TBP, which is a more polar solvent than
dodecane with its dielectric constant of 8.178 compared to 2.006 for dodecane.
Our DFT calculations with the PCM implicit solvent model show that the
potential energy (not free energy) of the cis conformation ($\Delta E_{cis-
trans}$) is stabilized by 9 kJ/mol in TBP compared to dodecane. However, it is
still challenging to use DFT calculations to map out the whole free energy
landscape for all the possible combinations of $\phi$ and $\psi$ that define
the configuration of the head group. Our MTD simulations for CMPO in TBP
(shown in the Figure 6 right panel) agree with the above DFT calculations in
that the bound cis conformation is more strongly stabilized by TBP compared to
dodecane, with the difference in free energy between the DFT-calculated bound
conformation decreasing by 12 kJ/mol going from dodecane to TBP. The
differences between this value and that obtained from DFT could result from
some combination of differences in the energetic description (empirical force
field versus quantum mechanical), the explicit versus implicit solvent, and
being a free energy rather than a potential energy. Another observation from
the free energy surface is that not only is the central “hot spot”
corresponding to the cis conformation greatly lowered in TBP solvent, the
overall shape of the free energy landscape is altered. Local minima in the
bottom left quadrant ($\phi=-90^{\circ}$, $\psi=-90^{\circ}$) in dodecane
shift to ($\phi=-50^{\circ}$, $\psi=-60^{\circ}$) in TBP. In addition, the
local minima in the top right quadrant ($\phi=60^{\circ}$, $\psi=90^{\circ}$)
expand in TBP. These small shifts indicate that certain CMPO conformations are
stabilized through molecular-level solvent-extractant interactions that may be
difficult to capture with mean-field implicit solvent methods. Overall, our
methodology of using CMD with MTD has shown its capability in investigating
not only extractant conformation energetics with different head groups but
also the solvent effect. Additionally, the conformational behavior of CMPO
potentially explains its superior separation performance over those of
malonamides.

###  3.6 Minimum free energy pathway

Although all energetic information should be contained in these 2D
Ramachandran plots, the interpretation of these maps is key to elucidate the
reaction rate and the relative stability between the start and end states. 42
To quantify the energetic difference between the bound and unbound states, we
identify the minimum free energy pathway (MFEP) from the 2D Ramachandran
plots. Briefly, the MFEP is a simplified reaction coordinate that connects two
configurations of interest on the free energy surface. Although the
identification of MFEP is non-trivial and remains research-active, here, we
apply the methodology developed by Fu et al., called MULE. Here, we give an
example of THMA to illustrate how we generate the MFEP from our energetic
landscape. As shown in Figure 7 left panel, using MULE, we identify four MFEPs
from the THMA bound state to the lowest energetic, free, and unbound
conformation. Each MFEP is associated with a free energy transition profile in
the right panel. Since these four MFEPs should be identical to each other, we
condense them into one MFEP per extractant molecule by taking the average of
these profiles. All other MFEPs on free-energy landscapes for different
extractant molecules are shown in SI.

Figure 7: Left: bound-to-unbound MFEPs identified on the free energy
landscape. Right: free energies associated with these four MFEPs.

Figure 8 summarizes all MFEPs from bound to unbound conformations for all
extractant molecules investigated in this work. For the bound state THMA to
MeTHMA, we can see that there is a 9 kJ/mol difference in free energy
(equivalent to 3.6 kT at room temperature). Compared to the experimental
distribution ratios for Eu(III) 12, this difference in energetics potentially
accounts for the one order of magnitude higher distribution ratio for THMA.
From MeTHMA to DiMeTHMA, the free-energy difference at the bound state
increases by 2 kJ/mol, consistent with the lower distribution ratio associated
with DiMeTHMA, although the magnitude of the decrease is either underestimated
here or partially attributable to another effect.12, 15 When we look at the
overall shape of MFEPs, we can see that local minima are formed for MeTHMA,
which is relevant for the transitional gauche conformation. Such a local
minimum does not exist for either THMA or DiMeTHMA, which indicates that the
middle alkyl substitution is influencing the conformational energetics and
free-energy transition pathways. From MeTHMA to MeDMDBMA, the lengths of the
side chains are reduced, which has an effect on the position of the local
minimum in the MFEP. For DMDBTDMA, its MFEP shape is similar to that of
MeDMDBMA, with an enhanced free energy barrier in the transition pathway,
suggesting that the length of the alkyl chain in the head group will influence
the molecular flexibility and kinetics of reorganization.

Figure 8: Bound-to-unbound MFEPs of all systems investigated.

For CMPO, a fundamentally different extractant molecule, the differences in
energy between the bound and free states in the MFEP profiles are generally
lower than those of all malonamides. The solvent effect on the MFEP is also
pronounced, where TBP not only lowers the bound-state energy but also biases
the overall MFEP. Much of the MFEP profile is higher in TBP than in dodecane,
resulting in a lower population of other conformations along the transition
pathway, suggesting more sluggish interconversion kinetics despite the overall
lower energy of the bound conformation. While this work mainly focuses on
exploring the application of the MTD advanced free energy sampling method in
quantifying extractant conformational energetics, more investigation of the
solvent effects and solvent-solute interactions will be explored in future
studies.

##  4 Conclusion

Traditional computational ligand design studies in liquid-liquid extraction of
metals heavily rely on cluster DFT calculations, in which the metal-extractant
binding energy is typically used as the primary means to evaluate ligand
performance. Such a strategy neglects the degrees of freedom of the extractant
structure in solution in favor of a single, lowest-energy structure. To
address this limitation, we implemented a combined approach with DFT and CMD
to study the impact of the extractant molecular structure on the energetic
penalty associated with organizing the extractant molecule into the metal-
binding conformation. An advanced sampling method, metadynamics, enabled
sampling of the entire conformational free energy landscape of a single metal-
free extractant molecule in solution. These CMD results provide entropic
information, directly reporting free energies instead of potential energies,
in addition to considering solvent effects explicitly. By using DFT
calculations to project bound-state structures onto the free-energy landscapes
sampled from CMD, we show how functionalization on the extractants’ head
groups reshapes the conformational energetics. These changes in the
conformational free-energy landscape potentially explain the relationship
between the functionalization of the extractant molecule and the observed
separation performance as measured by the metal-ion distribution ratios. As
the unfavorable conformational free energy associated with assuming the
binding conformation is added to the competing enthalpic favorability of the
metal-ligand binding, changes to the conformational free energy in the absence
of metal ions still provide important information on the separation behavior
of the extractant. Overall, this study highlights the opportunity to use CMD
simulations with the MTD technique to explore extractant behavior, learning
how to tune intramolecular conformational energetics to enhance separation
performance, while also explicitly accounting for solvent-solute interactions
that are present in practical LLE processes.

This work was supported by U.S. Department of Energy (DOE), Office of Science,
Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and
Biosciences Division, Separation Science Program, under Contract DE-
AC02-06CH11357 to UChicago Argonne, LLC, Operator of Argonne National
Laboratory. We gratefully acknowledge the computing resources provided on
Bebop, a high-performance computing cluster operated by the Laboratory
Computing Resource Center at Argonne National Laboratory.

MFEPs of all extractant molecules studied in this work.

## References

  * Hidayah and Abidin 2018 Hidayah, N. N.; Abidin, S. Z. The evolution of mineral processing in extraction of rare earth elements using liquid-liquid extraction: A review. _Minerals Engineering_ 2018, _121_ , 146–157
  * Špadina et al. 2021 Špadina, M.; Dufrêche, J.-F.; Pellet-Rostaing, S.; Marčelja, S.; Zemb, T. Molecular forces in liquid–liquid extraction. _Langmuir_ 2021, _37_ , 10637–10656
  * Hay et al. 2004 Hay, B. P.; Gutowski, M.; Dixon, D. A.; Garza, J.; Vargas, R.; Moyer, B. A. Structural criteria for the rational design of selective ligands: convergent hydrogen bonding sites for the nitrate anion. _Journal of the American Chemical Society_ 2004, _126_ , 7925–7934
  * Boehme and Wipff 2002 Boehme, C.; Wipff, G. Carbamoylphosphine oxide complexes of trivalent lanthanide cations: Role of counterions, ligand binding mode, and protonation investigated by quantum mechanical calculations. _Inorganic Chemistry_ 2002, _41_ , 727–737
  * Lewis et al. 2011 Lewis, F. W.; Harwood, L. M.; Hudson, M. J.; Drew, M. G.; Desreux, J. F.; Vidick, G.; Bouslimani, N.; Modolo, G.; Wilden, A.; Sypula, M., et al. Highly efficient separation of actinides from lanthanides by a phenanthroline-derived bis-triazine ligand. _Journal of the American Chemical Society_ 2011, _133_ , 13093–13102
  * McCann et al. 2016 McCann, B. W.; Silva, N. D.; Windus, T. L.; Gordon, M. S.; Moyer, B. A.; Bryantsev, V. S.; Hay, B. P. Computer-aided molecular design of bis-phosphine oxide lanthanide extractants. _Inorganic Chemistry_ 2016, _55_ , 5787–5803
  * Ellis et al. 2017 Ellis, R. J.; Brigham, D. M.; Delmau, L.; Ivanov, A. S.; Williams, N. J.; Vo, M. N.; Reinhart, B.; Moyer, B. A.; Bryantsev, V. S. “Straining” to separate the rare earths: how the lanthanide contraction impacts chelation by diglycolamide ligands. _Inorganic Chemistry_ 2017, _56_ , 1152–1160
  * Zhang et al. 2022 Zhang, X.; Adelman, S. L.; Arko, B. T.; De Silva, C. R.; Su, J.; Kozimor, S. A.; Mocko, V.; Shafer, J. C.; Stein, B. W.; Schreckenbach, G.; Batista, E. R.; Yang, P. Advancing the Am Extractant Design through the Interplay among Planarity, Preorganization, and Substitution Effects. _Inorganic Chemistry_ 2022, _61_ , 11556–11570
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proceedings of the national academy of sciences_ 2002, _99_ , 12562–12566
  * Barducci et al. 2011 Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2011, _1_ , 826–843
  * Valsson et al. 2016 Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint. _Annual review of physical chemistry_ 2016, _67_ , 159–184
  * McNamara et al. 1999 McNamara, B. K.; Lumetta, G. J.; Rapko, B. M. Extraction of europium (III) ion with tetrahexylmalonamides. _Solvent Extraction and Ion Exchange_ 1999, _17_ , 1403–1421
  * Lumetta et al. 1999 Lumetta, G. J.; McNamara, B. K.; Rapko, B. M.; Hutchison, J. E. Complexation of uranyl ion by tetrahexylmalonamides: an equilibrium modeling and infrared spectroscopic study. _Inorganica Chimica Acta_ 1999, _293_ , 195–205
  * Nigond et al. 1995 Nigond, L.; Condamines, N.; Cordier, P.; Livet, J.; Madic, C.; Cuillerdier, C.; Musikas, C.; Hudson, M. Recent advances in the treatment of nuclear wastes by the use of diamide and picolinamide extractants. _Separation Science and Technology_ 1995, _30_ , 2075–2099
  * Lumetta et al. 2002 Lumetta, G. J.; Rapko, B. M.; Garza, P. A.; Hay, B. P.; Gilbertson, R. D.; Weakley, T. J.; Hutchison, J. E. Deliberate design of ligand architecture yields dramatic enhancement of metal ion affinity. _Journal of the American Chemical Society_ 2002, _124_ , 5644–5645
  * Philip Horwitz et al. 1985 Philip Horwitz, E.; Kalina, D. C.; Diamond, H.; Vandegrift, G. F.; Schulz, W. W. The truex process-a process for the extraction of the tkansuranic elements erom nitric AC in wastes utilizing modified purex solvent. _Solvent Extraction and Ion Exchange_ 1985, _3_ , 75–109
  * Zhao and Truhlar 2008 Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. _Theoretical Chemistry Accounts_ 2008, _120_ , 215–241
  * Pritchard et al. 2019 Pritchard, B. P.; Altarawy, D.; Didier, B.; Gibson, T. D.; Windus, T. L. New basis set exchange: An open, up-to-date resource for the molecular sciences community. _Journal of Chemical Information and Modeling_ 2019, _59_ , 4814–4820
  * Frisch et al. 2016 Frisch, M. J. et al. Gaussian˜16 Revision C.01. 2016; Gaussian Inc. Wallingford CT
  * He et al. 2020 He, X.; Man, V. H.; Yang, W.; Lee, T.-S.; Wang, J. A fast and high-quality charge model for the next generation general AMBER force field. _The Journal of Chemical Physics_ 2020, _153_ , 114502
  * Jakalian et al. 2000 Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. _Journal of computational chemistry_ 2000, _21_ , 132–146
  * Jakalian et al. 2002 Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. _Journal of computational chemistry_ 2002, _23_ , 1623–1641
  * Jorgensen et al. 1996 Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. _Journal of the American Chemical Society_ 1996, _118_ , 11225–11236
  * Vo et al. 2015 Vo, Q. N.; Hawkins, C. A.; Dang, L. X.; Nilsson, M.; Nguyen, H. D. Computational study of molecular structure and self-association of tri-n-butyl phosphates in n-dodecane. _The Journal of Physical Chemistry B_ 2015, _119_ , 1588–1597
  * Martínez et al. 2009 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. _Journal of Computational Chemistry_ 2009, _30_ , 2157–2164
  * Martyna et al. 1994 Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. _The Journal of Chemical Physics_ 1994, _101_ , 4177–4189
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _Journal of Applied physics_ 1981, _52_ , 7182–7190
  * Shinoda et al. 2004 Shinoda, W.; Shiga, M.; Mikami, M. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. _Physical Review B_ 2004, _69_ , 134103
  * Thompson et al. 2022 Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D., et al. LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. _Computer Physics Communications_ 2022, _271_ , 108171
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Computer physics communications_ 2014, _185_ , 604–613
  * Ramachandran et al. 1963 Ramachandran, G.; Ramakrishnan, C.; Sasisekharan, V. Stereochemistry of polypeptide chain configurations. _Journal of Molecular Biology_ 1963, _7_ , 95–99
  * Ellis and Antonio 2012 Ellis, R. J.; Antonio, M. R. Coordination structures and supramolecular architectures in a cerium (III)–malonamide solvent extraction system. _Langmuir_ 2012, _28_ , 5987–5998
  * Ellis et al. 2013 Ellis, R. J.; Meridiano, Y.; Chiarizia, R.; Berthon, L.; Muller, J.; Couston, L.; Antonio, M. R. Periodic behavior of lanthanide coordination within reverse micelles. _Chemistry–A European Journal_ 2013, _19_ , 2663–2675
  * Ellis et al. 2014 Ellis, R. J.; Meridiano, Y.; Muller, J.; Berthon, L.; Guilbaud, P.; Zorz, N.; Antonio, M. R.; Demars, T.; Zemb, T. Complexation-Induced Supramolecular Assembly Drives Metal-Ion Extraction. _Chemistry–A European Journal_ 2014, _20_ , 12796–12807
  * Berthon et al. 2007 Berthon, L.; Martinet, L.; Testard, F.; Madic, C.; Zemb, T. Solvent penetration and sterical stabilization of reverse aggregates based on the DIAMEX process extracting molecules: Consequences for the third phase formation. _Solvent extraction and ion exchange_ 2007, _25_ , 545–576
  * Bonnett et al. 2023 Bonnett, B. L.; Sheyfer, D.; Wimalasiri, P. N.; Nayak, S.; Lal, J.; Zhang, Q.; Seifert, S.; Stephenson, G. B.; Servis, M. J. Critical fluctuations in liquid–liquid extraction organic phases controlled by extractant and diluent molecular structure. _Physical Chemistry Chemical Physics_ 2023, _25_ , 16389–16403
  * Erlinger et al. 1998 Erlinger, C.; Gazeau, D.; Zemb, T.; Madic, C.; Lefrancois, L.; Hebrant, M.; Tondre, C. Effect of nitric acid extraction on phase behavior, microstructure and interactions between primary aggregates in the system dimethyldibutyltetradecylmalonamide (DMDBTDMA)/n-dodecane/water: A phase analysis and small angle X-ray scattering (SAXS) characterisation study. _Solvent extraction and ion exchange_ 1998, _16_ , 707–738
  * Makombe et al. 2022 Makombe, E.; Bourgeois, D.; Berthon, L.; Meyer, D. Uranium (VI) and thorium (IV) extraction by malonamides: Impact of ligand molecular topology on selectivity. _Journal of Molecular Liquids_ 2022, _368_ , 120701
  * Johnson et al. 2023 Johnson, K. R.; Driscoll, D. M.; Damron, J. T.; Ivanov, A. S.; Jansone-Popova, S. Size Selective Ligand Tug of War Strategy to Separate Rare Earth Elements. _JACS Au_ 2023, _3_ , 584–591
  * Kalina et al. 1981 Kalina, D.; Horwitz, E.; Kaplan, L.; Muscatello, A. The extraction of Am (III) and Fe (III) by selected dihexyl N, N-dialkylcarbamoylmethyl-phosphonates,-phosphinates and-phosphine oxides from nitrate media. _Separation Science and Technology_ 1981, _16_ , 1127–1145
  * Horwitz et al. 1982 Horwitz, E.; Kalina, D.; Kaplan, L.; Mason, G.; Diamond, H. Selected alkyl (phenyl)-N, N-dialkylcarbamoylmethylphosphine oxides as extractants for Am (III) from nitric acid media. _Separation Science and Technology_ 1982, _17_ , 1261–1279
  * 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. _The journal of physical chemistry B_ 2005, _109_ , 6676–6687
  * Fu et al. 2020 Fu, H.; Chen, H.; Wang, X.; Chai, H.; Shao, X.; Cai, W.; Chipot, C. Finding an optimal pathway on a multidimensional free-energy landscape. _Journal of Chemical Information and Modeling_ 2020, _60_ , 5366–5374