# Are strongly confined colloids good models for two dimensional liquids?

Jiting Tian  Institute of Nuclear Physics and Chemistry, China Academy of
Engineering Physics, 621999 Mianyang, China  Univ. Grenoble-Alpes, CNRS,
LIPhy, 38000 Grenoble, France  Walter Kob<EMAIL_ADDRESS>Laboratoire Charles Coulomb (L2C), University of Montpellier and CNRS, F-34095
Montpellier, France  Jean-Louis Barrat  Univ. Grenoble-Alpes, CNRS, LIPhy,
38000 Grenoble, France  Institut Laue Langevin, 38000 Grenoble, France

###### Abstract

Quasi-two-dimensional (quasi-2D) colloidal hard-sphere suspensions confined in
a slit geometry are widely used as two dimensional (2D) model systems in
experiments that probe the glassy relaxation dynamics of 2D systems. However,
the question to what extent these quasi-2D systems indeed represent 2D systems
is rarely brought up. Here, we use computer simulations that take into account
hydrodynamic interactions to show that dense quasi-2D colloidal bi-disperse
hard-sphere suspensions exhibit much more rapid diffusion and relaxation than
their 2D counterparts at the same area fraction. This difference is induced by
the additional vertical space in the quasi-2D samples in which the small
colloids can move out of the 2D plane, therefore allowing overlap between
particles in the projected trajectories. Surprisingly, this difference in the
dynamics can be accounted for if, instead of using the surface density, one
characterizes the systems by means of a suitable structural quantity related
to the radial distribution function. This implies that in the two geometries
the relevant physics for glass-formation is essentially identical. Our results
provide not only practical implications on 2D colloidal experiments but also
interesting insights into the 3D-to-2D crossover in glass-forming systems.

quasi-2D, glass transition, colloid suspensions, simulations

††preprint: APS/123-QED

##  I Introduction

When cooled down fast enough, a liquid can avoid crystallization and form a
solid with a disordered structure. This phenomenon, termed the glass
transition, has been widely observed in many natural and industrial systems
but still lacks a satisfactory fundamental understanding [1, 2, 3]. Direct
spatiotemporal information of particles in glass-forming liquids is highly
useful for investigating the physical mechanisms leading to the glass
transition. However, in common glass formers, such as metallic or molecular
liquids, the relevant dynamics of the individual atoms or molecules is on the
scale of nanometers and picoseconds (or even smaller) [4], thus very difficult
to access by current experimental techniques. In contrast to this, colloidal
suspensions not only share typical features with normal glass-forming liquids
but also have the big advantage that the trajectories of the constituent
colloid particles can be directly tracked by optical imaging, therefore
allowing a detailed analysis of the system’s dynamics on the particle level
[4]. As a result, in recent years colloidal suspensions have been widely used
as model glass-forming systems in experiments and have provided insightful
results on the physics of glass formation [5, 6, 7, 8, 9, 10, 11, 12, 13].

Compared to three-dimensional (3D) colloidal suspensions, 2D systems are more
popular due to the lower difficulty in optical imaging [7, 14, 15, 8, 9, 10,
11, 12, 13]. However, given that phase transitions (such as melting) in 2D and
in 3D systems are very different [16, 17, 18, 19, 20], one has to wonder
whether 2D and 3D glass transitions have indeed the same nature. Using
computer simulations, Flenner and Szamel showed that many typical features of
glassy dynamics in 3D systems (e.g., the transient cage formation, the
coupling between orientational and translational relaxations) are absent in 2D
samples, indicating that the glass transitions in 2D and in 3D differ
significantly [21]. However, later experiments [8, 9] and simulations [22, 23]
have clearly demonstrated that this difference is induced by Mermin-Wagner
fluctuations which refer to the phenomenon that particles in 2D can move
significantly without obviously changing their local environments (neighbors).
Once the influence of these long-length fluctuations has been removed by using
cage-relative quantities, the 2D systems are found to display all typical
features of glassy dynamics, just like their 3D counterparts [8, 9, 22, 23].
Hence this issue is resolved and 2D colloidal suspensions can be used as
valuable models of glass-forming systems [24].

Figure 1: Static structures of 2D and of quasi-2D systems. (a) Snapshots of
the simulation boxes; (b) Radial distribution function for two typical area
fractions. The inset is a zoom on the grey shaded region of the main plot; (c)
Radial distribution function for quasi-2D $\phi=0.70$ and its partial
components.

While in computer simulations of 2D systems one usually uses a strictly 2D
geometry, experimentally probed 2D samples are often quasi-2D rather than
exactly 2D, except for a few cases where colloids are placed at a water-air
interface to mimic an exactly 2D geometry [25, 26]. This kind of quasi-2D
geometry is usually realized by sedimentation of particles in a solvent or by
confinement in a narrow space between two parallel plates, i.e., in a slit
geometry [6, 8, 12, 13]. In both cases, bi-disperse or poly-disperse particles
are not exactly moving in a 2D plane (that’s why they are called “quasi-2D”)
and therefore the question arises whether or not this difference between 2D
and quasi-2D affects the analysis and interpretation of experimental 2D glass
transitions.

Using the first type of setup, Thorneywork et al. have carried out a series of
experiments of binary quasi-2D colloidal hard-sphere suspensions prepared by
strong sedimentation. Their results showed a remarkable agreement with 2D
theories regarding static structures and diffusion dynamics, implying that the
difference between 2D and quasi-2D can be ignored when colloidal particles
strongly sediment at the bottom of the container [27, 28]. In the second case
of confined colloidal suspensions in which one often uses density-matched
solvents, even if the slit height is approaching the diameter of the largest
particles to realize an extreme confinement, the smaller particles do have a
significant freedom in the vertical direction (see the illustration in Fig.
1(a)). How the out-of-plane fluctuations in these confined quasi-2D systems
affect the glassy dynamics is still an open question.

Effects of confinement on the dynamics of colloidal suspensions (or hard-
sphere liquids) have been intensively studied in the past decades [29, 30, 6,
31, 32, 33, 34]. Schmidt and Löwen have used Monte Carlo simulations to
investigate the phase diagram of hard spheres confined between two parallel
plates [29, 30], but they focused only on the phase transitions in
monodisperse systems and did not address the glass transition problem. Nugent
et al. [6] carried out experiments to investigate the colloidal glass
transition in confinement and found that the dynamics is slower than the one
of the bulk suspension at the same packing fraction, i.e., in confined samples
the glass transition is shifted to lower packing fractions. Subsequent
computer simulations showed that the relaxation times of glassy liquids with
different degrees of confinement collapse onto a master curve if plotted as a
function of excess entropy, implying that relaxation of bulk and confined
viscous liquids are not fundamentally different [31]. These studies only
involved moderate confinements, where plate separations were larger than 4
particle diameters. Several recent works [32, 33, 34] applied strong
confinement (with plate separations between 2 and 5 diameters) to hard-sphere
liquids, finding an interesting non-monotonic behavior of diffusion as a
function of wall separation. However, none of the above studies considered
extremely confined suspensions or exactly 2D systems and therefore the
difference and similarity between 2D and quasi-2D still need to be assessed.

Here we perform computer simulations of bi-disperse, nearly-hard sphere
particles to investigate the difference between 2D and quasi-2D geometries.
Solvent-mediated hydrodynamic interactions (HI) and non-penetrable hard walls
are included into our model to mimic real colloidal suspensions in extreme
confinement. In contrast to from most previous simulations in which HI have
been ignored, our model with an adequate consideration of HI should give a
more realistic description of the dynamics of real colloidal glass-formers,
especially the short-time diffusion where experiments have clearly show that
HI play an important role [35]. This is in contrast to the commonly used
Langevin dynamics (LD) which cannot predict accurately the diffusion constant
of colloidal suspensions, as shown in Ref. [36] and confirmed by our
simulations (see Appendix). As we will show below, our results demonstrate not
only a notable difference in the relaxation dynamics between 2D and quasi-2D
at the same high area fraction, but also their equivalence when they are
characterized by a suitable structural or dynamic quantity.

##  II Methods and Model

There are many numerical approaches to simulate colloidal suspensions, mainly
differing from each other by the level of accuracy in describing the HI
between colloids and hence the resulting computational cost. The simplest
method is LD, or its overdamped variant Brownian dynamics (BD), in which one
just adds to the Newtonian equation of each colloid a linear Stokes friction
term and its corresponding random term to simulate the effect of the
surrounding solvent. Since LD ignores inter-particle HI, it is only a
reasonable description for very dilute suspensions. This limitation is avoided
in the so-called Stokesian dynamics method which is based on the linearity of
the Stokes equations for low Reynolds number (Re) fluids in that it describes
HI with a combination of near-field pairwise lubrication forces and far-field
many-body terms [37, 38, 39, 40]. For colloidal suspensions, which usually
have a very small Re, Stokesian dynamics is a very suitable simulation method,
although the computational complexity is typically super-linear [37, 38] (only
recently reduced to $O(N)$ [39, 40]). More complex situations (such as non-
spherical particles, irregular boundaries, etc) and high-Re fluids require
more advanced methods [36], including explicit-solvent approaches like
dissipative particle dynamics [41, 42], stochastic rotation dynamics [43], or
numerical simulation methods like fluid particle dynamics [44] and smoothed
profile method [45]. These techniques are more general and powerful, but also
more difficult to code and way more expensive in CPU time.

Since the present work focuses on the difference between the glassy dynamics
in 2D and in quasi-2D suspensions rather than the precise effects of HI, we
adopt a simplified version of Stokesian dynamics, namely the so-called “fast
lubricate dynamics” (FLD) [46], to simulate the motion of the colloids. This
method is a minimal model to include in an approximate manner hydrodynamic
effects while keeping the computational cost at a reasonable level, as
requested by our interest in slow dynamics. In our work, which involves no
external shear of the samples, FLD consists in adding a lubrication force
$\mathbf{F}_{\rm lub}$ and its corresponding random force $\mathbf{F}_{\rm
lub}^{\rm ran}$ to the total force in the normal Langevin equation, and hence
one has:

$m\dot{\mathbf{v}}=\mathbf{F}_{\rm p}+\mathbf{F}_{\rm Sto}+\mathbf{F}_{\rm
Sto}^{\rm ran}+\mathbf{F}_{\rm lub}+\mathbf{F}_{\rm lub}^{\rm ran},$ (1)

where $\mathbf{F}_{\rm p}$, $\mathbf{F}_{\rm Sto}$, and $\mathbf{F}_{\rm
Sto}^{\rm ran}$ are the inter-particle conservative force, the Stokes damping
force and its corresponding fluctuation term (random force), respectively.

In this work $\mathbf{F}_{\rm p}$ is obtained from the pseudo-hard-sphere
(PHS) potential [47], which has shown a good ability to reproduce typical
properties of single- and binary-component hard-sphere liquids. We consider a
binary system of particles which we denote by A and B and whose size ratio is
1:1.4. The pairwise energy of the potential is expressed as

$u_{ij}(r)=\begin{cases}50(\frac{50}{49})^{49}\epsilon[(\frac{\sigma_{ij}}{r})^{50}-(\frac{\sigma_{ij}}{r})^{49}]+\epsilon,&r<(\frac{50}{49})\sigma_{ij}\\\
0,&r\geqslant(\frac{50}{49})\sigma_{ij}\end{cases}$ (2)

where we define $\sigma_{ij}=0.5(\sigma_{i}+\sigma_{j})$, with $\sigma_{i}$
the diameter of particle $i$. The large exponent in the power laws ensures
that the potential is very steep and close to a hard core interaction [47]. We
choose $\epsilon=10k_{\rm B}T$ so that the particles do indeed probe the
strongly repulsive part of the potential. $\mathbf{F}_{\rm Sto}$ and
$\mathbf{F}_{\rm Sto}^{\rm ran}$ are the same as in the standard LD algorithm
[48], i.e.,

$\mathbf{F}_{\rm Sto}(i)=-3\pi\eta\sigma_{i}\mathbf{v}_{i},$ (3)

$\mathbf{F}_{\rm Sto}^{\rm ran}(i)=\sqrt{2k_{\rm B}T\times
3\pi\eta\sigma_{i}/\delta t}\;\mathbf{G}_{i},$ (4)

where $\eta$ is the viscosity of the solvent, $\mathbf{v}_{i}$ is the velocity
of particle $i$, $\delta t$ is the timestep in the simulation, and each
component of $\mathbf{G}_{i}$, $G_{i\alpha}$, is an independent Gaussian
random number with zero mean and unit variance, i.e., $\langle
G_{i\alpha}\rangle=0$ and $\langle
G_{i\alpha}G_{j\beta}\rangle=\delta_{ij}\delta_{\alpha\beta}$.

Figure 2: Comparing diffusion and relaxation of 2D and quasi-2D systems: (a)
MSD and (b) SISF as a function of $t$ for six typical $\phi$; the wave-vector
$k$ is $k_{0}=5.0$, close to the the peak of the structure factor $S(k)$; (c)
$D_{\rm L}/D_{0}$ and $\tau_{\alpha}$ as a function of $\phi$. The inset in
(a) displays the log-scale slopes of the solid lines in the main plot.

The lubrication force on particle $i$ (position $\mathbf{r}_{i}$, velocity
$\mathbf{v}_{i}$) generated by a particle $j$ within a cutoff distance $r_{\rm
c}^{\rm lub}=0.75(\sigma_{i}+\sigma_{j})$ is given by

$\mathbf{F}_{\rm lub}(i)=-\mathbf{F}_{\rm lub}(j)=-a_{\rm
sq}(\mathbf{v}_{i}-\mathbf{v}_{j})\cdot\mathbf{n}\mathbf{n}$ (5)

with

$a_{\rm
sq}=\frac{3}{2}\pi\eta(\frac{\sigma_{i}\sigma_{j}}{\sigma_{i}+\sigma_{j}})^{2}/h,$
(6)

where
$\mathbf{n}=(\mathbf{r}_{j}-\mathbf{r}_{i})/|\mathbf{r}_{i}-\mathbf{r}_{j}|$
is the center-to-center unit vector from particle $i$ to $j$, and
$h=|\mathbf{r}_{i}-\mathbf{r}_{j}|-0.5\sigma_{i}-0.5\sigma_{j}$ is the
surface-to-surface distance between the two particles. Here we ignore the
${\rm log}(h)$ terms in the lubrication force [49], keeping only the leading
order term $a_{\rm sq}$ which scales as $1/h$ and thus diverges when two
neighbors touch each other. It can be inferred from Eq. (5) that
$\mathbf{F}_{\rm lub}$ is repulsive when two particles are approaching each
other and attractive when separating, therefore reflecting the hydrodynamic
effects mediated by the solvent. The corresponding random force is related to
$\mathbf{F}_{\rm lub}$ through the fluctuation-dissipation theorem, i.e.,

$\mathbf{F}_{\rm lub}^{\rm ran}(i)=-\mathbf{F}_{\rm lub}^{\rm
ran}(j)=\sqrt{2k_{\rm B}Ta_{sq}/\delta t}\;G_{ij}\mathbf{n},$ (7)

where $G_{ij}$ is an independent Gaussian random number with zero mean and
unit variance for the pair $ij$. Note that $\mathbf{F}_{\rm lub}$ and
$\mathbf{F}_{\rm lub}^{\rm ran}$ always satisfy Newton’s third law, ensuring
momentum conservation.

FLD includes the near-field pairwise lubrication forces of Stokesian dynamics
in an exact manner but significantly simplifies the far-field many-body part
by approximating it by an isotropic resistance tensor (linear damping force),
therefore trading accuracy against computational efficiency. Since for
confined dense suspensions it can be expected that long range HI are screened
and hence weaker than the near-field lubrication forces, this approximation
should be reasonable for our purposes. Further support for this choice is the
earlier finding that FLD is able to reproduce both short-time and long-time
diffusion constants of equilibrium 3D monodisperse hard-sphere suspensions,
while LD or BD cannot [36]. Our test simulations lead to similar results (see
Appendix and Fig. A1), confirming the ability of FLD to describe the dynamics
of equilibrium colloidal suspensions. Finally we mention that we have also
performed test simulations with LD and obtained similar conclusions (see
Appendix and Fig. A3), indicating that our main conclusions regarding the
relation between 2D and quasi-2D systems are not sensitive to the details of
the microscopic dynamics.

At the start of our simulations the 50:50 mixture of spheres (mass ratio
$m_{\rm A}/m_{\rm B}=1/1.4^{3}$) with a total number of $N=2500$ particles are
inserted into a simulation box. As shown in Fig. 1(a), two kinds of boxes are
used: In the strictly 2D systems, particles can move only in the plane, i.e.
$z=0$ (black dashed line), while in quasi-2D boxes, particles can move in 3D
but they are confined between two parallel hard walls that have a separation
of $H=1.1\sigma_{\rm B}$, i.e. a value which is close to the extreme
confinement case accessed by experiments. Whenever a particle hits a wall, it
undergoes a specular reflection, i.e., its velocity in the $z$ direction
changes sign. For both geometries, the box is a square with periodic boundary
conditions in the $x$ and $y$ directions. Denoting by $L$ the length of the
box, we define the area packing fraction as
$\phi=(N/2)\pi[(\sigma_{A}/2)^{2}+(\sigma_{B}/2)^{2}]/L^{2}$ and we consider
value between $0.40\leq\phi\leq 0.785$. In order to improve the statistics of
the results we average them for each $\phi$ and each geometry over 20
independent samples.

In our simulations, $\sigma_{A}$, $\tau$, and $k_{\rm B}T$ are used as reduced
units for length, time, and energy, respectively. Here
$\tau=\sigma_{A}^{2}/4D_{0}$ is the Brownian timescale, i.e. the average time
it takes a small particle to diffuse (in the dilute limit) over a distance
give by its size, and $D_{0}=k_{\rm B}T/3\pi\eta\sigma_{A}$ is the single-
particle diffusion constant of the small particles. For colloidal suspensions
one is in an overdamped regime, i.e. the mass of the particle should be
sufficiently small that inertial effects can be neglected. We have performed
test simulations of monodisperse 3D systems with a series of particle masses
and found that for FLD the diffusion dynamics already converges at $m=0.01$
(see Fig. A1 in Appendix), therefore we use in the following this value for
the mass.

From an initial configuration, we equilibrate the system for at least 10
$\tau_{\alpha}$ (defined below) and then perform a production run of length
1000 $\tau$. The timestep $\delta t$ for all simulations is $10^{-4}\ \tau$.
Our test simulations (results not shown here) demonstrate that a smaller
timestep ($10^{-5}\ \tau$) does not give different results and that $10^{-4}\
\tau$ gives an excellent energy conservation under $NVE$ condition with our
pseudo hard sphere system.

Previous studies have shown that in 2D systems the usual mean squared
displacement of a tagged particle is not able to characterize well the slowing
down of the dynamics of the system because of the presence of long wave-length
fluctuations discussed in the Introduction [21]. To quantify the motion of the
particles, we therefore use the cage-relative (CR) displacement
$\Delta\mathbf{r}_{\rm CR}(t)$, Refs. [8, 9, 22, 23], which is defined as the
displacement of a particle relative to its neighbors at time zero, i.e.

$\Delta\mathbf{r}_{\rm CR}(t)=\mathbf{r}(t)-\mathbf{r}(0)-\frac{1}{N_{\rm
nei}}\sum_{j}[\mathbf{r}_{j}(t)-\mathbf{r}_{j}(0)].$ (8)

Here $j$ runs over all the nearest neighbors of the particle in the initial
configuration (of number $N_{\rm nei}$), determined by the Voronoi cell
method. From $\Delta\mathbf{r}_{\rm CR}(t)$ one can obtain the mean squared
displacement (MSD), $\langle\Delta r^{2}_{\rm CR}(t)\rangle$, and the self
intermediate scattering function (SISF) for wave-vector $\mathbf{k}$, $F_{\rm
s,CR}(k,t)=\langle{\rm exp}({\rm i}\mathbf{k}\cdot\Delta\mathbf{r}_{\rm
CR}(t)\rangle$, to characterize the diffusion and relaxation dynamics,
respectively. Here $\langle.\rangle$ indicates an average over all particles
and over the different production runs. For $F_{\rm s,CR}(k,t)$, an additional
average is performed over all wave-vectors $\mathbf{k}$ that have modulus
$k_{0}=5.0\sigma^{-1}$. From $\langle\Delta r^{2}_{\rm CR}(t)\rangle$ and
$F_{\rm s,CR}(t)$ we can get, respectively, the long-time diffusion constant
$D_{\rm L}=\lim_{t\rightarrow\infty}\langle\Delta r^{2}_{\rm CR}(t)\rangle/4t$
and the $\alpha$-relaxation time $\tau_{\alpha}$ defined by $F_{\rm
s,CR}(k,\tau_{\alpha})=1/e$ (horizontal black dashed line in Figs. 2(b) and
Fig. 3(b)).

As in experiments, only the $x$ and $y$ coordinates are used in the data
analysis, $i.e.$ $\mathbf{r}_{i}(t)=(x_{i}(t),y_{i}(t),0)$ for any particle
$i$. Similarly, $\mathbf{k}=(k_{x},k_{y},0)$ is used when calculating the
static structure factor $S(k)$ and $F_{\rm s,CR}(k,t)$. Only for the dashed
lines in Fig. 6 where the real dynamics in 3D of quasi-2D samples are
calculated, we use $\mathbf{r}_{i}(t)=(x_{i}(t),y_{i}(t),z_{i}(t))$.

##  III Results and Discussions

To compare the structure of the two types of systems we show in Fig. 1(b)
their radial distribution function $g(r)$. For the moderate area fraction
$\phi$=0.60, $g(r)$ of the 2D and quasi-2D systems are basically identical,
except for $r<1.0$ where the allowed (projected) overlap between particles in
quasi-2D leads to a non-zero particle density (see the grey area and inset).
(The same positive values of $g(r)$ for $r<1.0$ have also been observed in
experimental quasi-2D colloidal suspensions [13].) As $\phi$ is increased to
0.76, the difference between the two structures becomes much more notable in
that the peaks in ${\rm g}(r)$ for 2D are higher than the ones in the quasi-2D
system, implying a denser packing. Considering that a denser packing
corresponds to slower dynamics in colloidal glass-forming systems, we can
expect that at high density, the quasi-2D will relax quicker than the 2D
systems at the same area fraction, and below we will see that this is indeed
the case. Figure 1(c) shows an example of decomposing g$(r)$ into its partial
components A-A, A-B, and B-B [2]. One recognizes that the first peak of $g(r)$
is solely due to the A-A pairs, the second is dominated by A-B pairs, while
the third one has contributions from the A-B as well as the B-B pairs. Below
we will come back to this result.

Figure 3: Comparing diffusion and relaxation of A and B particles in 2D and
quasi-2D systems: (a) MSD and (b) SISF as a function of $t$ for two values of
$\phi$; (c) Enhancement factor (defined in the text) of $D_{\rm L}/D_{0}$ and
of $\tau_{\alpha}$ as a function of $\phi$.

In Fig. 2 we compare the relaxation dynamics of the two systems by monitoring
the time dependence of the MSD and the SISF. Panel (a) shows the expected
slowing down of the diffusion as $\phi$ increases, independent of the geometry
considered. For low values of $\phi$ we find at short times a ballistic regime
that crosses over at around $t=0.1$ to a diffusive one, i.e. the MSD growth
linearly in time. For intermediate and high $\phi$ we see that the particles
are caged on intermediate time scales, i.e. the MSD is sub-diffusive. These
different regimes can be recognized better by monitoring the (double-
logarithmic) slope of the MSD, $\gamma(t)$, show in the inset of panel (a) for
the case of the quasi-2D system. At $\phi=0.40$ the system behaves like a
normal colloidal liquid and does not display caging. At $\phi=0.60$ it already
enters the glassy regime, with a minimum of $\gamma(t)$ somewhat lower than 1,
and with increasing $\phi$ this minimum becomes significantly deeper.

When we compare the diffusion dynamics of 2D and of quasi-2D systems at the
same area fraction (solid and dashed lines), it is clear that the former is
slower and that the difference increases gradually from insignificant at
$\phi=0.40$ to quite obvious at $\phi=0.78$, demonstrating that at high
packing fraction the dynamics of the two system at the same $\phi$ becomes
very different.

Figure 2(b) demonstrates that the intermediate scattering function sends the
same message as the MSD, i.e. that the dynamics of the 2D system becomes
quickly slower than the one of the quasi-2D system. Furthermore we recognize
that the shape of the correlator at long times is independent of the packing
fraction, in agreement with earlier results on 3D systems [50], and, more
importantly, is independent of the geometry considered.

To quantify the difference in the dynamics between the 2D and quasi-2D systems
we present in Fig. 2(c) the $\phi-$dependence of the normalized diffusion
constant $D_{\rm L}/D_{0}$, where the long and short time diffusion constants
$D_{\rm L}$ and $D_{0}$ have been defined above. From the graph one recognizes
that while at low $\phi$ this ratio is independent of the system, the relative
difference grows quickly with increasing $\phi$ in that the curve for the 2D
system decreases faster than the one for the quasi-2D system. The same
conclusion is reached by monitoring the $\alpha$-relaxation time
$\tau_{\alpha}$, introduced above, in that the curve for the 2D system
increases much quicker than the one for the quasi-2D system. (For the 2D case,
$\phi=0.78$ has a relaxation time larger than $10^{3}$ and therefore we do not
report the values of $D_{\rm L}/D_{0}$ and $\tau_{\rm\alpha}$). Note that,
both for 2D and quasi-2D, $\phi=0.76$ corresponds to a moderately supercooled
state, as evidenced by a relaxation time that is less than 500, and therefore
we expect that the difference in relaxation times will be much larger at
higher area fractions. Overall, these results clearly indicate that the
dynamics in 2D and in quasi-2D systems are very different at the same area
fraction when entering the supercooled, caged regime.

Figure 4: Unifying diffusion and relaxation dynamics in 2D and quasi-2D
systems: $D_{\rm L}/D_{0}$ and $\tau_{\alpha}$ as a function of (a) $g(r_{\rm
3rd})$ and (b) $\gamma_{\rm min}$; (c) MSD and SISF as a function of $t$ for
$g(r=r_{\rm 3rd})\approx 1.7$, 3.0 and 3.5, from normal to glassy (the
corresponding area fractions are given in the legend).

Since in the quasi-2D system the large (type B) particles are strongly
confined, one might expect that the differences in the dynamics of the two
types of systems are mainly due to the small (type A) particles because in the
quasi-2D system they can move vertically. To check whether this is indeed the
case, we present in Fig. 3(a) the mean squared displacement of the two types
of particles. For low packing fraction we see that the A particles are moving
faster than the B particles and that this dynamics is independent of the
geometry considered. More surprisingly is the observation that at high packing
fraction, $\phi=0.76$, there is no significant decoupling between the dynamics
of the big and small particles, even in the quasi-2D system. The inset shows
the ratio of the MSD of quasi-2D over that of 2D, and one recognizes that this
ratio is independent of the type of particle considered, i.e. the acceleration
of the dynamics in the quasi-2D system for the A and the B particles is the
same. The same conclusion is reached by probing the intermediate scattering
function for the two types of particles and geometries, see Fig. 3(b), showing
that not only the fast component of the dynamics (which is dominating the MSD)
is independent of the particle type, but also the slow component (which
dominates $F_{s,CR}(k,t)$).

To characterize the $\phi-$dependence of the dynamics for the two types of
systems on a quantitative level we define, in analogy with the inset in Fig.
3(a), an enhancement factor $\kappa_{D}$ ($\kappa_{\tau}$) as the ratio of
$D_{\rm L}/D_{0}$ ($\tau_{\rm\alpha}$) of the quasi-2D and 2D systems and show
the results for the A and B particles in Fig. 3(c). We see that all
enhancement factors increase from 1 to larger than 10 as $\phi$ increases, and
that the data for the A particles tracks closely the one for the B particles,
only showing a small difference for at $\phi=0.76$. This surprising
independence of the particle type can be understood from a simple argument:
When a A particle moves out of the $x-y$ plane, its effective diameter in this
plane decreases thus leaving more space for the neighboring B particles and
hence allowing the latter to move faster, even if their motion along the
$z$-axis is highly restricted.

The fact that the differences between the 2D and quasi-2D systems regarding
their structure and dynamics increase gradually from negligible at small
$\phi$ to significant at high $\phi$, hints the existence of a direct link
between structure and dynamics for both 2D and quasi-2D dynamics systems.
(Similar connections have been proposed for bulk systems at different specific
state points [51] but here we refer to relate systems with different
geometries to each other.)

Since the height of the peaks in $g(r)$ gives direct information about the
strength of the local structural order and hence of the cage, this height is a
good candidate to connect structure to dynamic properties. In practice we find
that the height of the third peak in g($r$), i.e. $g(r_{\rm 3rd})$, provides
such a link, as shown in Fig. 4(a): Systems with different geometries, but
with the same value of $g(r_{\rm 3rd})$, have the same diffusion constant and
relaxation time, thus demonstrating that there is an intimate connection
between the two types of systems. In contrast to this mapping involving
$g(r_{\rm 3rd})$, other simple structural quantities like the effective area
fraction $\Phi$ calculated by the effective particle diameters obtained from
the positions of the peaks in $g(r)$ [52, 12, 53], or the height of the other
two peaks in $g(r)$ (i.e. $g(r_{\rm 1st})$ and $g(r_{\rm 2nd}$)) do not allow
to obtain a good mapping (see Fig. A2 in the Appendix). Since $g(r_{\rm 3rd})$
is dominated by the contributions from the A-B and B-B pairs, see Fig. 1(c),
it contains most of the relevant structural information concerning the slow
component, i.e. the B particles. This is probably the reason why such a simple
quantity can be used to unify the relaxation dynamics for different
geometries. Finally we mention that we have also tested whether such a
unification can be obtained using the height of the main peak in the partial
radial distribution function for the B-B pairs, i.e., the maximum in $g_{\rm
BB}(r)$, but did not obtain a satisfactory result (see Fig. A2(d) in the
Appendix), implying that the structural information associated with the A-B
pairs is crucial for getting a reliable mapping.

Figure 5: ($\phi_{\rm 2D}$,$\phi_{\text{q-2D}}$) pairs whose corresponding
systems show the same dynamics. Figure 6: Comparing projected and real
diffusion dynamics in quasi-2D systems: MSD as a function of $t$ for (a) all,
(b) A, and (c) B particles.

Apart from the structural quantities discussed above, we find that also short-
time dynamic quantities can be used to predict the long-time diffusion and
relaxation behavior. Inspired by the recent experimental finding [8] regarding
the one-to-one relationship between the minimum value of the logarithmic MSD
slope $\gamma_{\rm min}$ and the relaxation time $\tau_{\alpha}$ for 2D and 3D
colloidal suspensions, we plot in Fig. 4(b) $D_{\rm L}/D_{0}$ and
$\tau_{\alpha}$ as a function of $\gamma_{\rm min}$. The graph demonstrates
that $\gamma_{\rm min}$ allows to make a map between the two systems, similar
to the case of $g(r_{\rm 3rd})$. Since $\gamma_{\rm min}$ reflects the extent
of the transient localization of colloids due to “caging”, the unification in
Fig. 4(b) implies that long-time diffusion and relaxation are determined by
the short-time caged state both for 2D and quasi-2D systems.

Combining the results from Figs. 4(a) and (b), one concludes that there exists
a one-to-one mapping between $g(r_{\rm 3rd})$ and $\gamma_{\rm min}$, and that
this mapping is independent of the geometry of the system. Since $g(r_{\rm
3rd})$ can predict the value of the short time quantity $\gamma_{\rm min}$ as
well as the long time quantities $D_{\rm L}/D_{0}$ and $\tau_{\rm\alpha}$, one
can envision the possibility that $g(r_{\rm 3rd}$) determines the entire
diffusion and relaxation processes. This is indeed the case, as shown in Fig.
4(c), where we present three sets of curves that correspond to different
values of $g(r_{\rm 3rd})$. One recognizes that, at a fixed $g(r_{\rm 3rd})$,
the MSD for the 2D and quasi-2D systems are very close to each other, and the
same is true for $F_{s,{\rm CR}}(k,t)$. (The small difference for g($r_{\rm
3rd})\approx 3.5$ is due to the high sensitivity of dynamics to g($r_{\rm
3rd})$ for dense suspensions and the difficulty to prepare a pair of 2D and
quasi-2D systems that have exactly the same g($r_{\rm 3rd})$.) From Fig. 4 one
thus concludes that the glassy dynamics in 2D and in quasi-2D liquids are
identical when probed at the same value of the relevant parameter, i.e.
$g(r_{\rm 3rd})$.

Figure 4 shows that for each packing fraction $\phi_{2D}$ of the 2D system
there is a fraction $\phi_{q-2D}$ for the quasi-2D system at which the two
dynamics are identical. Using a linear interpolation of the data points in
Fig. 4(a) one thus can obtain the correspondence between $\phi_{\rm 2D}$ and
$\phi_{\text{q-2D}}$ and the result is shown in Fig. 5. For packing fractions
corresponding to a normal liquid, e.g. $\phi_{\rm 2D}=0.40$, one has
$\phi_{\rm 2D}\approx\phi_{\text{q-2D}}$, in agreement with Fig. 4(c). For
systems that have a glassy dynamics, i.e. packing densities around 0.6 and
higher (grey region in Fig. 5), we find the data to be compatible with the
simple relation $\phi_{\text{q-2D}}-\phi_{\rm 2D}\approx 0.02$, i.e. a
constant shift in the packing fraction (dashed line). Note that such a
constant shift implies that with increasing $\phi$ the difference in
relaxation times of the two systems grows rapidly, because $\tau_{\alpha}$ is
a non-linear function of $\phi$, in agreement with Figs. 2 (a) and (b). We
also mention that the dynamics of quasi-2D systems will strongly depend on the
slit height $H$, i.e. the difference 0.02 reported above can be expected to
increase with increasing $H$. It will therefore be interesting to see in the
future whether or not a universal mapping for various $H$ exists, so that
experimental quasi-2D results can be easily compared to simulated or
theoretical (strictly) 2D results.

Although in experiments one measures projected (2D) coordinates, it should be
kept in mind that the real dynamics occurs in 3D (or more precisely, quasi-2D)
and hence one can wonder how the projected and the real dynamics differ from
each other. In Fig. 6(a) we thus compare projected and real MSD averaged over
all particles. It is clear that the two quantities coincide at long times,
i.e. once the horizontal displacements significantly exceed the vertical ones.
There are, however, obvious differences in the transient caged regime, i.e.
when $t\approx 1$, in particular for high area fractions, in that the real
dynamics shows a plateau in the MSD that is somewhat more extended (i.e. the
system is more glassy) and significantly higher (larger cage size) than the
MSD for the projected dynamics. This difference is more pronounced for A
particles, Fig. 6(b), than for the B particles, Fig. 6(c), due to the larger
available space. Hence Fig. 6 demonstrates that while the real and projected
dynamics at short and long times are very similar, there are significant
differences in the caged regime and therefore one should keep in mind these
differences when interpreting experimental data.

##  IV Summary and Conclusions

To answer the simple question of whether experimental quasi-2D colloidal
suspensions can represent faithfully 2D model systems, we have performed
extensive computer simulations to investigate to what extend the structural
and dynamic properties of the two types of systems differ from each other. For
this we have used a simple but realistic model for the colloidal particles and
a dynamics which includes the part of the hydrodynamic interactions that is
the most relevant for dense suspensions. We find, as expected, that the
quasi-2D samples exhibit faster dynamics than the 2D ones at the same area
fraction, reflecting the enhanced freedom associated with the displacement of
the small particles in the transverse direction. The acceleration of the
dynamics by more than a factor of 10 already for a moderate concentration of
$\phi=0.76$ is, however, a surprisingly large effect. A closer inspection of
the relaxation dynamics reveals that it is possible to map the dynamics of one
system onto the dynamics of the other one. This mapping can either be done by
using the third peak in $g(r$) or the minimum value of the double-logarithmic
slope of the MSD. The existence of such a mapping allows thus to establish a
connection between the dynamics of quasi-2D as they are experimentally
accessible to the dynamics of strictly 2D systems as they are considered in
simulations and analytical calculations. In the future it will be of interest
to probe whether such a mapping can also be done for other types of
interactions and how its details depend on the height of the confinement
space.

###### Acknowledgements.

We thank Bo Li and Eric Weeks for sharing with us their data and for useful
discussions. JT acknowledges the financial support by the China Scholarship
Council (CSC) (No. 201904890006) and his employer INPC for the strong support
which makes his stay in Grenoble possible. WK is a senior member of the
Institut Universitaire de France.

## Appendix

Figure A1: Results of test simulations of 3D monodisperse colloidal
suspensions: (a) The normalized short time diffusion coefficient $D_{\rm
S}/D_{0}$ and (b) the long time diffusion coefficient $D_{\rm L}/D_{0}$ as a
function of $\phi$ for LD and FLD and a series of masses. Also included are
the results from Ref. [36].

Real colloids in a solvent show an overdamped dynamics and hence inertia
effects should be absent. In order to check whether this is indeed the case we
perform test simulations with 256 monodisperse 3D spheres with the pseudo hard
sphere potential described in the main text. (Note that the original FLD
method uses a volume-fraction-dependent Stokes damping force, i.e.
$\mathbf{F}_{\rm Sto}=-3\pi\eta\sigma_{i}\mathbf{v}_{i}f(\phi)$ where
$f(\phi)=1+2.16\phi$ is fitted from Stokesian dynamics simulation results.
Since this factor was specifically obtained for 3D systems, we include it in
our test simulations of 3D systems, but not for 2D and for quasi-2D systems.)

Using LD or FLD and a series of masses we determine the short-time and the
long-time diffusion constants, i.e. $D_{\rm S}/D_{0}$ and $D_{\rm L}/D_{0}$,
where the former is defined as the maximum of $D(t)/D_{0}$, where
$D(t)=\langle\Delta r^{2}(t)\rangle/4t$ is the time-dependent self-diffusion
coefficient (see Ref. [36] for details). We compare our results for different
packing fractions with those in the literature in Fig. A1. From panel (a) one
concludes that for the FLD a decreasing mass results in an increasing $D_{\rm
S}/D_{0}$, finally converging for a mass around $m=0.01$. After a finite-size
correction [54],

$D_{\rm S}=D_{\rm S}(N)+D_{0}[1.76(\phi/N)^{1/3}-\phi/N]$ (9)

where $N=256$, the results are very close to the values from Ref. [36] which
agree well with experimental data. In contrast to this, the results for the LD
show no convergence and do not coincide with the reference values,
demonstrating that this type of microscopic dynamics is not suitable to
describe the dynamics of real colloidal systems.

In Fig. A1(b), we show the long time diffusion coefficient and we recognize
that also for this quantity the FLD can give a converged and accurate estimate
of $D_{\rm L}/D_{0}$ while this is not the case for the LD. Our results thus
confirm the findings of Ref. [36], i.e. FLD is better suited than LD to
describe diffusion in hard sphere-like colloidal suspensions.

Figure A2: Test whether various structural quantities can be used to establish
a mapping between the quasi-2D and the 2D systems. $D_{\rm L}/D_{0}$ as a
function of (a) effective area packing fraction $\Phi$, (b) $g(r_{\rm 1st}$),
(c) $g(r_{\rm 2nd}$), and (d) $g_{\rm BB}(r_{\rm peak})$.

As mentioned in the main text, we have also tested whether a different type of
structural quantity can be used to make the mapping from the quasi-2D system
to the 2D one. In Fig. A2 we show that this is not the case if one considers
the effective packing fraction $\Phi$, the height of the first or second peak
in $g(r)$, or the height of the sole peak in the partial distribution function
$g_{BB}(r)$. One thus concludes that for this mapping to work, one needs
information about the A-B correlation.

To test how the results presented in the main text depend on the microscopic
dynamics we carry out test simulations of the 2D and of quasi-2D systems with
LD, and show the main results in Fig. A3. We can see that the differences and
the equivalences between 2D and quasi-2D systems simulated by LD are similar
to those obtained for the FLD in the main text, except that for the LD the
mapping using $g(r_{\rm 3rd})$, shown in Fig. A3(b), is of inferior quality.

Figure A3: Results of test simulations of 2D and quasi-2D systems with LD:
$D_{\rm L}/D_{0}$ and $\tau_{\alpha}$ as a function of (a) $\phi$, (b)
$g(r_{\rm 3rd}$) and (c) $\gamma_{\rm min}$.

## References

  * Debenedetti and Stillinger [2001] P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001).
  * Binder and Kob [2011] K. Binder and W. Kob,  _Glassy materials and disordered solids: An introduction to their statistical mechanics_ (World Scientific, 2011).
  * Berthier and Biroli [2011] L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
  * Hunter and Weeks [2012] G. L. Hunter and E. R. Weeks, Rep. Prog. Phys. 75 (2012).
  * Weeks and Weitz [2002] E. R. Weeks and D. A. Weitz, Phys. Rev. Lett. 89 (2002).
  * Nugent  _et al._ [2007] C. R. Nugent, K. V. Edmond, H. N. Patel, and E. R. Weeks, Phys. Rev. Lett. 99 (2007).
  * Tanaka  _et al._ [2010] H. Tanaka, T. Kawasaki, H. Shintani, and K. Watanabe, Nat. Mater. 9, 324 (2010).
  * Vivek  _et al._ [2017] S. Vivek, C. P. Kelleher, P. M. Chaikin, and E. R. Weeks, Proc. Natl. Acad. Sci. U.S.A. 114, 1850 (2017).
  * Illing  _et al._ [2017] B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, Proc. Natl. Acad. Sci. U.S.A. 114, 1856 (2017).
  * Ma  _et al._ [2019] X. Ma, Z. S. Davidson, T. Still, R. J. Ivancic, S. S. Schoenholz, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett. 122 (2019).
  * Yang  _et al._ [2019] X. Yang, H. Tong, W. H. Wang, and K. Chen, Phys. Rev. E 99 (2019).
  * Lozano  _et al._ [2019] C. Lozano, J. R. Gomez-Solano, and C. Bechinger, Nat. Mater. 18, 1118 (2019).
  * Li  _et al._ [2020] B. Li, K. Lou, W. Kob, and S. Granick, Nature 587, 225 (2020).
  * Mishra and Ganapathy [2015] C. K. Mishra and R. Ganapathy, Phys. Rev. Lett. 114, 198302 (2015).
  * Nagamanasa  _et al._ [2015] K. H. Nagamanasa, S. Gokhale, A. Sood, and R. Ganapathy, Nature Physics 11, 403 (2015).
  * Kosterlitz and Thouless [1972] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 5, L124 (1972).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 6, 1181 (1973).
  * Halperin and Nelson [1978] B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
  * Nelson and Halperin [1979] D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979).
  * Young [1979] A. P. Young, Phys. Rev. B 19, 1855 (1979).
  * Flenner and Szamel [2015] E. Flenner and G. Szamel, Nat. Commun. 6 (2015).
  * Shiba  _et al._ [2016] H. Shiba, Y. Yamada, T. Kawasaki, and K. Kim, Phys. Rev. Lett. 117 (2016).
  * Shiba  _et al._ [2019] H. Shiba, T. Kawasaki, and K. Kim, Phys. Rev. Lett. 123 (2019).
  * Tarjus [2017] G. Tarjus, Proc. Natl. Acad. Sci. U.S.A. 114, 2440 (2017).
  * Chen  _et al._ [2006] W. Chen, S. Tan, Z. Huang, T.-K. Ng, W. T. Ford, and P. Tong, Phys. Rev. E 74, 021406 (2006).
  * Ma  _et al._ [2013] X. Ma, W. Chen, Z. Wang, Y. Peng, Y. Han, and P. Tong, Phys. Rev. Lett. 110, 078302 (2013).
  * Thorneywork  _et al._ [2014] A. L. Thorneywork, R. Roth, D. G. Aarts, and R. P. Dullens, J. Chem. Phys. 140 (2014).
  * Thorneywork  _et al._ [2018] A. L. Thorneywork, S. K. Schnyder, D. G. Aarts, J. Horbach, R. Roth, and R. P. Dullens, Mol. Phys. 116, 3245 (2018).
  * Schmidt and Löwen [1996] M. Schmidt and H. Löwen, Phys. Rev. Lett. 76, 4552 (1996).
  * Schmidt and Löwen [1997] M. Schmidt and H. Löwen, Phys. Rev. E 55, 7228 (1997).
  * Ingebrigtsen  _et al._ [2013] T. S. Ingebrigtsen, J. R. Errington, T. M. Truskett, and J. C. Dyre, Phys. Rev. Lett. 111 (2013).
  * Mittal  _et al._ [2008] J. Mittal, T. M. Truskett, J. R. Errington, and G. Hummer, Phys. Rev. Lett. 100 (2008).
  * Lang  _et al._ [2010] S. Lang, V. Boţan, M. Oettel, D. Hajnal, T. Franosch, and R. Schilling, Phys. Rev. Lett. 105 (2010).
  * Mandal  _et al._ [2014] S. Mandal, S. Lang, M. Gross, M. Oettel, D. Raabe, T. Franosch, and F. Varnik, Nat. Commun. 5 (2014).
  * Thorneywork  _et al._ [2015] A. L. Thorneywork, R. E. Rozas, R. P. Dullens, and J. Horbach, Phys. Rev. Lett. 115, 268301 (2015).
  * Bolintineanu  _et al._ [2014] D. S. Bolintineanu, G. S. Grest, J. B. Lechman, F. Pierce, S. J. Plimpton, and P. R. Schunk, Comput. Part. Mech. 1, 321 (2014).
  * Bossis and Brady [1984] G. Bossis and J. F. Brady, J. Chem. Phys. 80, 5141 (1984).
  * Sierou and Brady [2001] A. Sierou and J. F. Brady, J. Fluid Mech. 448, 115 (2001).
  * Fiore and Swan [2019] A. M. Fiore and J. W. Swan, J. Fluid Mech. 878, 544 (2019).
  * Ouaknin  _et al._ [2021] G. Y. Ouaknin, Y. Su, and R. N. Zia, J. Comput. Phys. 442 (2021).
  * Hoogerbrugge and Koelman [1992] P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992).
  * Soddemann  _et al._ [2003] T. Soddemann, B. Dünweg, and K. Kremer, Phys. Rev. E 68 (2003).
  * Malevanets and Kapral [1999] A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605 (1999).
  * Tanaka and Araki [2000] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
  * Nakayama and Yamamoto [2005] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71 (2005).
  * Kumar and Higdon [2010] A. Kumar and J. J. Higdon, Phys. Rev. E 82 (2010).
  * Jover  _et al._ [2012] J. Jover, A. J. Haslam, A. Galindo, G. Jackson, and E. A. Müller, J. Chem. Phys. 137 (2012).
  * Allen and Tildesley [2017] M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ (Oxford University Press, 2017).
  * Ball and Melrose [1997] R. C. Ball and J. R. Melrose, Physica A 247, 444 (1997).
  * Kob [1999] W. Kob, Journal of Physics: Condensed Matter 11, R85 (1999).
  * Bacher  _et al._ [2014] A. K. Bacher, T. B. Schrøder, and J. C. Dyre, Nat. Comm. 5, 1 (2014).
  * de J. Guevara-Rodríguez and Medina-Noyola [2003] F. de J. Guevara-Rodríguez and M. Medina-Noyola, Phys. Rev. E 68, 8 (2003).
  * Abaurrea-Velasco  _et al._ [2020] C. Abaurrea-Velasco, C. Lozano, C. Bechinger, and J. D. Graaf, Phys. Rev. Lett. 125 (2020).
  * Ladd [1990] A. J. Ladd, J. Chem. Phys. 93, 3484 (1990).


.gifpng.pngconvert gif:#1 png: .gif

††thanks: Both authors contributed equally††thanks: Both authors contributed
equally

# Adaptive atomic basis sets

Danish Khan  Chemical Physics Theory Group, Department of Chemistry,
University of Toronto, St. George Campus, Toronto, ON, Canada  Vector
Institute for Artificial Intelligence, Toronto, ON, M5S 1M1, Canada
Maximillian L. Ach  Department of Physics, University of Toronto, St. George
Campus, Toronto, ON, Canada  Department of Physics, Ludwig-Maximilians-
Universität München (LMU), Munich, Germany  O. Anatole von Lilienfeld
<EMAIL_ADDRESS>Chemical Physics Theory Group, Department of
Chemistry, University of Toronto, St. George Campus, Toronto, ON, Canada
Vector Institute for Artificial Intelligence, Toronto, ON, M5S 1M1, Canada
Acceleration Consortium, University of Toronto, Toronto, ON, Canada
Department of Physics, University of Toronto, St. George Campus, Toronto, ON,
Canada  Department of Materials Science and Engineering, University of
Toronto, St. George Campus, Toronto, ON, Canada  Machine Learning Group,
Technische Universität Berlin and Institute for the Foundations of Learning
and Data, 10587 Berlin, Germany  Berlin Institute for the Foundations of
Learning and Data, 10587 Berlin, Germany

###### Abstract

## Abstract

Atomic basis sets are widely employed within quantum mechanics based
simulations of matter. We introduce a machine learning model that adapts the
basis set to the local chemical environment of each atom, prior to the start
of self consistent field (SCF) calculations. In particular, as a proof of
principle and because of their historic popularity, we have studied the
Gaussian type orbitals from the Pople basis set, i.e. the STO-3G, 3-21G, 6-31G
and 6-31G*. We adapt the basis by scaling the variance of the radial Gaussian
functions leading to contraction or expansion of the atomic orbitals. A data
set of optimal scaling factors for C, H, O, N and F were obtained by
variational minimization of the Hartree-Fock (HF) energy of the smallest 2500
organic molecules from the QM9 database. Kernel ridge regression based machine
learning (ML) prediction errors of the change in scaling decay rapidly with
training set size, typically reaching less than 1% for training set size 2000.
Overall, we find systematically lower variance, and consequently the larger
training efficiencies, when going from hydrogen to carbon to nitrogen to
oxygen. Using the scaled basis functions obtained from the ML model, we
conducted HF calculations for the subsequent 30’000 molecules in QM9. In
comparison to the corresponding default Pople basis set results we observed
improved energetics in up to $99\%$ of all cases. With respect to the larger
basis set 6-311G(2df,2pd), atomization energy errors are lowered on average by
$\sim$31, 107, 11, and 11 kcal/mol for STO-3G, 3-21G, 6-31G and 6-31G*,
respectively — with negligible computational overhead. We illustrate the high
transferability of adaptive basis sets for larger out-of-domain molecules
relevant to addiction, diabetes, pain, aging. Generally, deviation from cc-
pvQZ based HF results for counter-poise corrected atomization energies, HOMO-
LUMO gaps, and dipole moments also decreases — in line with the variational
principle.

##  I Introduction

The choice of atomic basis set used for expanding molecular orbitals is one of
the central components of any quantum chemical calculation due to its large
influence on both computational cost and accuracy1. Since their popularization
in the early 1970s, Gaussian Type Orbitals (GTOs) have been by far the most
widely used type2; 3. The primary reason for this is the more efficient
evaluation of four-center integrals due to the Gaussian product rule. Since
the quality of the wavefunction and the cost of a quantum mechanical
calculation directly depend on the size and quality of the basis set, several
strategies and basis set families have been developed to improve this trade-
off. A detailed documentation can be found in numerous, comprehensive
reviews1; 4; 5. Since optimizing the basis set individually for each molecule
is impractical, the free parameters, namely the contraction coefficients and
the radial exponents, are in the majority of cases obtained by fitting to
atomic calculations and kept constant for each element6; 3; 7; 8; 9. Despite
the large success of this approach, it is known that optimal values can differ
significantly for molecular systems3; 10. While using such immutable basis
sets generally does not pose a problem when large (and balanced) basis sets of
triple zeta quality or better are employed, for more incomplete (and
frequently used due to tractability) basis sets there is still room for
improvement.
In analogy to the recently introduced machine learning (ML) assisted adaptive
hybrid density functional approximations11, we introduce ML models that
provide improved basis sets (STO-3G2, 3-21G6, 6-31G7, 6-31G*12) by adapting
them, on-the-fly, to the system of interest. This is done by scaling the
radial exponents of contracted GTOs (CGTOs) by a constant factor depending on
the local chemical environment of each atom. A training dataset of optimal
scaling factors is obtained via variational minimization of Hartree-Fock (HF)
energies of 2,000 smallest molecules from the QM913 dataset. The data is
subsequently used to train a ML model of the optimal scaling factors for each
atomic basis set. For out of sample molecules, this leads to improved HF total
and atomization energies in up to $\sim$99.9$\%$ of cases at no additional
computational cost. Since the optimal scaling factors are obtained via
application of the variational principle, the improved basis sets lead to
improved approximations of the exact ground state HF wavefunctions. As a
result of this, simlutaneous improvements are also obtained for other physical
properties like HOMO-LUMO gaps and dipole moments indicating a foundational ML
model for quantum chemistry. Due to the strong locality of the scaling
factors, the model is shown to be easily transferable to larger systems of
biologically relevant molecules suggesting modest training data needs to be
generally applicable.

##  II Theory and methods

###  II.1 Basis function scaling

Gaussian type basis functions are composed of linear combinations of primitive
Gaussian type orbitals (PGTOs) of the form

$\displaystyle g_{ik}^{\text{PGTO}}(\alpha_{ik})$
$\displaystyle=N_{ik,n}Y_{lm}(\theta,\phi)|\mathbf{r}-\mathbf{R}_{A}|^{n-1}\mathrm{e}^{(-\alpha_{ik}|\mathbf{r}-\mathbf{R}_{A}|^{2})},$
(1) $\displaystyle N_{ik,n}$
$\displaystyle=\sqrt{\frac{2(2\alpha_{ik})^{n+\frac{1}{2}}}{\Gamma(n+\frac{1}{2})}}$

which form contracted orbitals (CGTOs)

$\displaystyle\phi_{k}^{\text{CGTO}}=\sum_{i}c_{ik}~{}g_{ik}(\alpha_{ik})$ (2)

where $Y_{lm}$ are spherical harmonics, $\mathbf{R}_{A}$ denotes the nuclear
coordinate on which the PGTO/CGTO is centered, $c_{ik}$ are fixed (pre-
optimized) linear combination coefficients (called contraction coefficients),
$\Gamma$ denotes the gamma function and the sum runs over the length of
contraction. To make these basis functions adaptive to the local atomic
environment, the radial exponents of all PGTOs within a CGTO are multiplied by
a single constant factor $\zeta_{k}$, i.e.

$\displaystyle g_{ik,\mathrm{scaled}}^{\mathrm{PGTO}}(\zeta_{k}\alpha_{ik})=$
$\displaystyle N_{ik,nlm}Y_{lm}(\theta,\phi)R_{l}(\mathbf{r},\mathbf{R}_{A})$
(3)
$\displaystyle\cdot\mathrm{e}^{(-\zeta_{k}\alpha_{ik}|\mathbf{r}-\mathbf{R}_{A}|^{2})}$

and

$\displaystyle\phi_{k,\mathrm{scaled}}^{\mathrm{CGTO}}=\zeta_{k}^{\frac{3}{2}}\sum_{i}c_{ik}~{}g_{ik}(\zeta_{k}\alpha_{ik})$
(4)

where the prefactor is for normalization. $\zeta_{k}$ are unique for each CGTO
and assumed to depend on the atom’s local chemical environment. This is
similar to the non-adaptive scaling suggested for the 4-31G basis set3. The
local atomic environment is defined as the sphere
$S_{A}(\mathbf{R}_{A},r_{\mathrm{cut}})$ of radius $r_{\mathrm{cut}}$ centered
at $\mathbf{R}_{A}$. The $\zeta_{k}$ are then assumed to be functions of the
atomic species present within the sphere $S_{A}$,
$\zeta_{k}=f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon S_{A}})$.
The function $f$ in our work is estimated using ML. In our work we modify
valence CGTOs and polarization functions and use the same scaling factors for
s and p orbitals, since the exponents between these are shared in the default
basis sets as well. We note that including more scaling factors could however
further improve the method.
The optimal set of scaling factors for a molecule,
$\boldsymbol{\zeta}_{\text{opt}}$, is obtained by minimizing the HF energy as
a function of these $\zeta$

$\displaystyle\boldsymbol{\zeta}_{\text{opt}}=\arg\min_{\boldsymbol{\zeta}}E_{\text{HF}}(\boldsymbol{\zeta})$
(5)

Gradients of the HF energy with respect to $\zeta_{k}$ can be calculated using
the general HF response expression

$\displaystyle\frac{\partial E_{\text{HF}}}{\partial\phi_{k}}$
$\displaystyle=\sum_{\mu\nu}\biggl{(}P_{\mu\nu}\frac{\partial
h_{\mu\nu}}{\partial\phi_{k}}-W_{\mu\nu}\frac{\partial
S_{\mu\nu}}{\partial\phi_{k}}\biggr{)}$
$\displaystyle+\frac{1}{2}\sum_{\mu\nu\lambda\sigma}P_{\mu\nu}P_{\lambda\sigma}\frac{\partial(\mu\nu||\lambda\sigma)}{\partial\phi_{k}}$
$\displaystyle=2\sum_{\nu}\biggl{(}P_{k\nu}\frac{\partial
h_{k\nu}}{\partial\phi_{k}}-W_{k\nu}\frac{\partial
S_{k\nu}}{\partial\phi_{k}}\biggr{)}$ (6)
$\displaystyle+\sum_{\nu\lambda\sigma}P_{k\nu}P_{\lambda\sigma}\frac{\partial(k\nu||\lambda\sigma)}{\partial\phi_{k}}$
(7)

and

$\displaystyle\frac{\partial\phi_{k}}{\partial\zeta_{k}}=\frac{3}{2}\frac{\phi_{k}}{\zeta_{k}}-\zeta_{k}^{\frac{3}{2}}|\mathbf{r}-\mathbf{R}_{A}|^{2}\sum_{i}c_{ik}\alpha_{ik}~{}g_{ik}(\zeta_{k}\alpha_{ik})$
(8)

where $P_{\mu\nu}$, $S_{\mu\nu}$, $W_{\mu\nu}$, $h_{\mu\nu}$ denote elements
of the density, overlap, energy weigthed density and core Hamiltonian matrices
respectively9. It should be noted that an additional term due to $\zeta_{k}$
will appear when gradients such as atomic forces are evaluated. To ensure
smoothness of the forces, the following form can be used for changing
geometries:

$\displaystyle\zeta_{k}$
$\displaystyle=1+f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})\sum_{j=1}^{n}f_{\mathrm{cut}}(r_{jA})$ (9) $\displaystyle
f_{\mathrm{cut}}(r_{jA})$
$\displaystyle=\begin{cases}\frac{1}{2}\left(\cos(\frac{\pi
r_{jA}}{r_{\mathrm{cut}}})+1\right),~{}r_{jA}<r_{\mathrm{cut}}\\\
0,~{}r_{jA}>r_{\mathrm{cut}}\end{cases}$ (10)

where $f_{\mathrm{cut}}$ is a soft cut-off function that smoothly decays to 0
at $r_{\mathrm{cut}}$14 and $n$ denotes the number of atoms in the system.
Nuclear gradients can be evaluated as
$\nabla_{\mathbf{R}_{A}}\zeta_{k}\partial_{\zeta_{k}}\phi_{k}\partial_{\phi_{k}}E_{HF}$
using eq. 7 and 8. along with

$\displaystyle\nabla_{\mathbf{R}_{A}}\zeta_{k}=\nabla_{\mathbf{R}_{A}}f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})\sum_{j=1}^{n}f_{\mathrm{cut}}(r_{jA})$
$\displaystyle-\frac{1}{2}f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})\sum_{j:\mathbf{R}_{j}\epsilon S_{A}}\frac{\pi\sin(\frac{\pi
r_{jA}}{r_{\mathrm{cut}}})}{r_{\mathrm{cut}}r_{jA}}(\mathbf{R}_{A}-\mathbf{R}_{j})$
(11)

where the term
$\nabla_{\mathbf{R}_{A}}f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})$ needs to be evaluated by propagating the gradient through the ML
model15.
In the current work we have employed numerical gradients (central finite
differences with step size $10^{-8}$) for simplicity. The L-BFGS16 optimizer,
as implemented in Scipy17, was used to obtain optimal $\zeta_{k}$ of the
training set molecules by minimizing HF total energies. HF calculations were
performed using version 2.3.0 of the PySCF package18; 19. All basis sets were
obtained from basis set exchange20.

Figure 1: Optimal scaling factor, $\zeta_{\mathrm{opt}}$, distributions for
Hydrogen, Carbon, Nitrogen and Oxygen valence shell orbitals from the first
2000 smallest QM913 molecules found via minimization of total HF energies
employing numerical gradients. Data used for training, see learning curves in
Fig 2.

Using this methodology, we calculated optimal scaling factors of the first
2000 molecules from QM913, a dataset of 134k small organic molecules with up
to nine heavy atoms out of C,O,N,F. For STO-3G, this was additionally done for
30k more molecules. The distribution of $\zeta_{\mathrm{opt}}$ can be seen in
Fig. 1 for each chemical element and contracted GTO.

###  II.2 Machine Learning

The ML model used throughout this work is Kernel Ridge Regression21; 22; 23
(KRR) due to its high robustness and simplicity, which has been covered
extensively24; 25; 26; 27; 11. Using this method, the function $f$ is
approximated as a weighted sum of similarity measures to all atoms (of the
same type) in the training set

$\displaystyle f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})=\sum_{J}^{N_{\mathrm{train}}}\alpha_{J}k(\mathbf{x}(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}}),\mathbf{x}_{J})$ (12)

where $\alpha_{J}$ are the regression weights, $\mathbf{x}$ is the atomic
representation feature vector which depends only on the set of nuclear charges
and coordinates within the atom’s local environment $S_{A}$
($\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon S_{A}}$), and $k(.,.)$
denotes a kernel function acting as a similarity measure. The kernel function
used in our work is the atomic Gaussian kernel

$\displaystyle
k(\mathbf{x}_{I},\mathbf{x}_{J})=\exp{\left(-\frac{||\mathbf{x}_{I}-\mathbf{x}_{J}||^{2}_{2}}{2\sigma^{2}}\right)}$
(13)

The gradient
$\nabla_{\mathbf{R}_{A}}f(\\{Z_{i},\mathbf{R}_{i}\\}_{i:\mathbf{R}_{i}\epsilon
S_{A}})$ in eq. 11 is then straightforward to evaluate as

$\displaystyle\nabla_{\mathbf{R}_{A}}f=-\nabla_{\mathbf{R}_{A}}\mathbf{x}_{I}\sum_{J}^{N_{\mathrm{train}}}\alpha_{J}\frac{||\mathbf{x}_{I}-\mathbf{x}_{J}||_{2}}{\sigma^{2}}k(\mathbf{x}_{I},\mathbf{x}_{J})$
(14)

The regression weights $\mathbf{\alpha}$ are obtained from the set of training
labels $\mathbf{y}^{\mathrm{train}}$ via the following equation

$\mathbf{\alpha}=(\mathbf{K}+\lambda\cdot\mathbf{I})^{-1}\mathbf{y}^{\mathrm{train}}$
(15)

where $\mathbf{K}$ is the kernel matrix of the training set and $\lambda$ is a
regularization parameter.
In our work the training and predictions are performed separately for each
unique chemical element. This leads to improved predictions since the basis
set scalings have different variances for each chemical element (cf. Fig. 1).
Moreover, this also leads to a reduction of the size of the kernel matrix
$\mathbf{K}$ from $(\sum_{i}n_{i})^{2}$ to $(\max n_{i})^{2}$, where $n_{i}$
denotes number of atoms belonging to chemical element $i$ in the training set.
This results in more efficient training and predictions since the
computational cost and memory requirement of the method scale with the size of
$\mathbf{K}$.

In this work we have relied upon an extended version of the Many-Body
Distribution Functionals (MBDF) representation for generating atomic feature
vector mappings $\mathbf{x}$28. The reason for this choice is its compact
size, along with high predictive power, which makes the impact on the quantum
chemistry calculation cost negligible. Briefly, MBDF generates the feature
vector components as suitable functionals of two- and three-body
distributions:

$\displaystyle F_{\textrm{2}}^{nm}[i]$
$\displaystyle=\int_{0}^{\infty}dr~{}g_{n}(r)~{}\partial_{r}^{m}\rho_{i}(r),$
(16) $\displaystyle\rho_{i}(r)$ $\displaystyle=\sum_{j:\mathbf{R}_{j}\epsilon
S_{i}}^{M}~{}Z_{j}\mathcal{N}{\left(R_{ij},\sigma_{r}\right)}$ $\displaystyle
F_{\textrm{3}}^{nm}[i]$
$\displaystyle=\int_{0}^{\pi}d\theta~{}g_{n}(\theta)~{}\partial_{\theta}^{m}\rho_{i}(\theta),$
(17) $\displaystyle\rho_{i}(\theta)$
$\displaystyle=\sum_{j:\mathbf{R}_{j}\epsilon S_{i},k:\mathbf{R}_{k}\epsilon
S_{i}}^{M}~{}\frac{(Z_{j}Z_{k})^{\frac{1}{2}}\mathcal{N}\left(\theta_{ijk},\sigma_{\theta}\right)}{(R_{ij}R_{jk}R_{ik})^{2}}$

where $\mathcal{N}$ denotes a normalized Gaussian, $Z_{j}$ corresponds to the
nuclear charge of atom $j$, $R_{ij}$ and $\theta_{ijk}$ denote interatomic
distances and angles respectively and $g_{n}$ are suitable weighting functions
similar to 2, 3-body interaction potentials. Gradients
$\nabla_{\mathbf{R}_{A}}\mathbf{x}_{I}$ of the generated feature vectors (eq.
14) are implemented in the MBDF code28 (cf. Data and code availability).

##  III Results and Discussion

###  III.1 Numerical results

As outlined above, optimal scaling factors were calculated for the first 2500
molecules from the QM913 dataset, and then randomly drew 500 for testing and
up to 2000 of the remaining for training. The data set comprises molecules
with a 14669 hydrogen, 8860 carbon, 2356 oxygen, 1713 nitrogen and 22 fluorine
atoms. Scaling factors for the latter varied only negligibly and were hence
neglected for the rest of the study. During the optimization, we generally
found $\zeta$ to be well behaved and convex, making it suitable for machine
learning. Further, we note that a larger and better sampled training set is
likely to make the model even more transferable through the use of atomic
fragments based on amons29.

Figure 2: Mean absolute errors (MAE) of predicted scaling factors $\zeta$ for
STO-3G (valence orbitals), 3-21G, 6-31G (both inner and outer valence) and
6-31G* (inner and outer valence and polarization functions) as a function of
training set size (number of molecules) using the KRR-MBDF model on a
validation set of 500 out-of-sample QM9 molecules. Note that the average
number of atoms per molecule is different for each chemical element (7.33
Hydrogen, 4.43 Carbon and $\sim$1 Nitogen/Oxygen atoms per molecule). Training
sets were obtained from the $\zeta_{\mathrm{opt}}$ values calculated for the
first 2000 smallest QM9 molecules. Refer to Fig. 1 for training data
distributions.

The scaling factors were found to be machine-learnable quite efficiently using
our very simple KRR-MBDF model as shown in the learning curves in Fig. 2 for a
set of 500 out-of-sample molecules from the QM9 dataset. Lower offsets for
Nitrogen and Oxygen were generally observed throughout likely due to the small
variance of $\zeta_{\mathrm{opt}}$ values (Fig. 1) for these elements. Larger
offset and more erratic behaviour for the 6-31G* polarization function
learning is likely because there are 3 scaling factors in this basis set per
atom which are coupled to each other whereas the ML models are separate for
each. Distributions in Fig. 1 indicate that variance of the optimal scaling
factor is inversely proportional to the electronegativity of the element which
is in line with chemical intuition.

Figure 3: Change in Hartree Fock total energy for a validation set of 500
QM913 molecules compared to the default unscaled form as a function of
training set size.

This indicates that the method should be quite transferable with modest
training data needs for such elements. For Carbon and Hydrogen as well we
found the prediction errors to be sufficiently small so as to lead to
considerable improvements for out-of-sample molecules. This was verified by
using the ML model to predict basis sets adapted to 30,000 larger molecules
drawn at random from the QM9 dataset. Fig. 3 shows changes in total energies
obtained from HF calculations when using the adaptive basis sets as a function
of training set size. Remarkably, it can be seen that the adaptive basis sets
already start improving upon the default Pople basis sets with $\sim$100
training molecules and only $\sim$500 molecules are sufficient to recover the
vast majority of the possible accuracy increase.

| Total Energies | Atomization Energies | Average CPU time
---|---|---|---
Basis Set |  Improvement rate |  Average improvement [kcal/mol] |  Improvement rate |  Average improvement [kcal/mol] |  Default (adaptive) [sec/molecule]
aSTO-3G | 98.6% | 8.56 (9.53) | 82.9% | 31.2 | 0.590 (0.590)
a3-21G | 99.8% | 39.9 | 96.7% | 107.8 | 1.167 (1.167)
a6-31G | 99.2% | 3.29 | 99.9% | 11.9 | 1.204 (1.205)
a6-31G* | 99.2% | 2.37 | 93.7% | 11.0 | 1.464 (1.465)

Table 1: Hartree-Fock total energy decrease (Left), atomization energy
improvement vs 6-311G(2df,2pd) (Middle) and average computational cost
comparison when employing adaptive vs default basis sets for 30,000 QM9
molecules (Right). For aSTO-3G, the improvement in the total energy was
additionally computed using optimal (non-ML) scaling factors and are given in
brackets. All timings evaluated on a compute node equipped with a 16-core AMD
Ryzen 9 7950X @ 5.7GHz processor and 128 GB DDR5 RAM.

These metrics are quantified in Table 1 which shows changes in the energetics
obtained from using adaptive basis sets predicted by the ML model trained on
the first 2000 smallest QM9 molecules. For all 4 basis sets, the adaptive
versions lead to lowering of the total HF energies for almost all of the
30,000 tested molecules. Lowering of the total HF energy indicates a better
approximation to the optimal single Slater-determinant ground state
wavefunction due to the variational principle. The success rate is
particularly high for double-zeta basis sets, which is likely due to the
higher flexibility offered by having two scaling factors compared to STO-3G,
where all valence orbitals are scaled equally. For STO-3G, we computed optimal
scaling factors $\boldsymbol{\zeta}_{\text{opt}}$ for all 30k validation
molecules and found that our ML model on average recovered 90% of the energy
improvement (8.56/9.53 kcal/mol) that can be achieved using
$\boldsymbol{\zeta}_{\text{opt}}$ obtained via optimization. The corresponding
changes in atomization energies are slightly less consistent likely due to
basis set superposition errors (BSSE) which were not present in the training
data since the $\boldsymbol{\zeta}_{\text{opt}}$ were obtained by minimizing
molecular total energies. Nevertheless, the magnitude of the improvements can
be seen to be even larger than for total energies. Throughout our work, we
found particularly large improvements in energy for the 3-21G basis set. A
likely explanation for this is the small size of the basis set along with the
previously mentioned flexibility offered by two separate scaling factors. We
note here that, since smaller basis sets yield the largest accuracy increase,
adaptive basis sets should be beneficial to Hamiltonians used in quantum
computers where only small basis sets are employed due to hardware
constraints30. We also note again that the optimal scalings for polarization
functions (Fig. 2) were found to be harder to learn which can likely be
improved with a more balanced training set selection or a ML model that
accounts for the coupling between the scaling factors into account by learning
them together. Larger improvements from the HF calculations should be possible
with more accurate polarization function scalings since these are known to
have a drastic impact 1. In Table 1, we also include average CPU timings for
the HF calculations using adaptive vs default basis sets. For adaptive basis
sets, this contains the scaling factor prediction time by the ML model and the
HF calculation time. As can be seen, the increased accuracy using adaptive
basis sets can be obtained at virtually no additional cost.

| Atomization energy error [kcal/mol]
---|---
Molecule |  STO-3G |  aSTO-3G |  3-21G |  a3-21G |  6-31G |  a6-31G |  6-31G* |  a6-31G* |  6-31G**
Aspirin | 221.1 | 218.6 | 229.5 | 190.4 | 182.2 | 172.7 | 23.4 | 17.2 | 11.5
Glucose | 207.1 | 113.3 | 253.6 | 207.6 | 198.4 | 188.9 | 37.5 | 8.7 | 10.4
MDMA | 7.6 | 68.6 | 224.5 | 159.3 | 177.2 | 168.2 | 26.0 | 15.0 | 9.4
Mescaline | 37.5 | 109.7 | 248.9 | 181.6 | 197.9 | 185.6 | 31.3 | 16.8 | 11.2
Metformin | 114.9 | 99.3 | 188.4 | 165.2 | 150.3 | 144.5 | 29.6 | 24.8 | 11.3
Methylglyoxal conformer 1 | 83.6 | 84.3 | 99.7 | 84.8 | 86.6 | 83.2 | 8.7 | 4.6 | 4.4
Methylglyoxal conformer 2 | 78.5 | 82.5 | 100.1 | 83.7 | 82.1 | 78.2 | 8.6 | 4.1 | 4.3
Nicotine | 39.7 | 13.6 | 181.7 | 135.9 | 139.0 | 133.4 | 20.2 | 13.2 | 5.9
Paracetamol | 123.2 | 118.0 | 181.8 | 155.4 | 136.5 | 132.8 | 23.4 | 17.7 | 9.1
Resveratrol | 172.6 | 196.4 | 265.5 | 226.4 | 190.1 | 186.6 | 36.1 | 31.3 | 13.6
Salicylic acid | 175.2 | 153.6 | 176.3 | 158.2 | 136.3 | 132.4 | 21.8 | 19.3 | 9.2
Uric acid | 372.0 | 258.3 | 245.6 | 227.5 | 201.7 | 191.3 | 30.5 | 33.6 | 20.3
Mean Error | 136.1 | 126.4 | 199.6 | 164.7 | 156.1 | 149.4 | 24.8 | 17.2 | 10.0

Table 2: Counterpoise corrected Hartree-Fock atomization energy errors
compared to 6-311G(2df,2pd) for a test set of 12 out-of-domain
biologically/societally relevant molecules using adaptive and default Pople-
style basis sets. The Methylglyoxal conformers correspond to the aldehyde
group being on the same (1) and opposite side (2) of the keto group. Figure 4:
Hartree-Fock atomization energy mean absolute errors vs average computational
cost for a test set of 12 biologically/societally relevant molecules (Table 2)
using adaptive and default Pople-style basis sets. All timings were evaluated
on a compute node equipped with a 16-core AMD Ryzen 9 7950X @ 5.7GHz processor
and 128 GB DDR5 RAM.

To further test the transferability of the method we designed a test set of 12
larger, out-of-domain, biologically and societally relevant molecules. All
adaptive basis sets for these molecules were obtained through predictions from
the ML model trained on the 2000 smallest QM9 molecules. Table 1 in the SI
shows the change in total HF energy achieved using adaptive basis sets for
these molecules while atomization energy errors are shown in table 2.
Atomization energies of these molecules are compared directly to
6-311G(2df,2pd) and in order to mitigate BSSE, a counterpoise correction was
applied to all basis sets. As expected, the error declines with increasing
basis set size with the exception of STO-3G and for nearly all molecules, it
consistently improves when adaptive basis sets are employed. The better
performance of STO-3G compared to 3-21G and 6-31G was verified using the Orca
quantum chemistry package (version 5.0.3)31 and is likely linked to error
cancellation between the total and atomic energies. As mentioned earlier, the
large impact of polarization functions can be seen as the a6-31G* basis set
nearly reduces the error 10-fold compared to its upolarized form 6-31G. The
improvements obtained from each adaptive basis set in this case as well comes
at virtually no added computational cost as shown in the computational cost vs
accuracy plot in Fig. 4. Similar to the observation from QM9 testing,
particularly large improvement can be seen with the adaptive a3-21G basis but
also with the inclusion of polarization functions in the a6-31G* basis set. We
expect to see a similar shift in the Pareto front for even larger and well
balanced basis sets with more polarization functions being adaptable.

| HOMO-LUMO Gap [eV] | Dipole Moment Norm [D]
---|---|---
Molecule |  3-21G/a3-21G |  6-31G/a6-31G |  6-31G*/a6-31G* |  cc-pVQZ |  3-21G/a3-21G |  6-31G/a6-31G |  6-31G*/a6-31G* |  cc-pVQZ
Aspirin | 11.65/11.8 | 11.46/11.48 | 11.53/11.53 | 11.31 | 1.97/1.94 | 2.32/2.34 | 2.24/2.16 | 2.23
Glucose | 15.07/15.4 | 14.92/14.97 | 15.15/14.42 | 11.31 | 4.21/4.17 | 4.49/4.43 | 3.85/3.83 | 3.46
MDMA | 11.99/12.3 | 11.93/11.85 | 11.81/11.8 | 11.03 | 1.56/1.54 | 1.79/1.66 | 1.22/1.24 | 1.06
Uric Acid | 11.54/11.6 | 11.39/11.42 | 11.58/11.6 | 10.89 | 3.36/3.43 | 3.50/3.56 | 3.42/3.52 | 3.55
Mescaline | 12.17/12.5 | 12.11/12.08 | 12.1/12.1 | 10.8 | 5.53/5.35 | 6.03/5.47 | 4.71/4.52 | 4.24
Metformin | 13.32/13.55 | 13.06/13.1 | 13.44/13.45 | 11.56 | 4.47/4.32 | 4.56/4.52 | 4.45/4.45 | 4.38
Nicotine | 12.68/12.99 | 12.46/12.43 | 12.32/12.3 | 11.86 | 3.32/2.96 | 3.63/3.42 | 3.08/3.04 | 3.04
Paracetamol | 11.59/11.92 | 11.56/11.45 | 11.49/11.43 | 10.54 | 2.63/2.62 | 2.88/2.89 | 2.57/2.57 | 2.63
Resveratrol | 9.47/9.5 | 9.4/9.36 | 9.41/9.4 | 9.18 | 3.96/3.82 | 4.11/4.08 | 3.69/3.67 | 3.52
Salicylic Acid | 11.43/11.62 | 11.28/11.26 | 11.29/11.25 | 11.04 | 3.32/3.27 | 3.83/3.81 | 3.38/3.31 | 3.35
MAE | 0.87/1.02 | 0.74/0.71 | 0.79/0.71 | - | 0.38/0.32 | 0.58/0.47 | 0.15/0.12 | -

Table 3: HOMO-LUMO gaps and dipole moment norms using default and adaptive
Pople-style basis sets. Mean absolute errors (MAE) are compared to the cc-pVQZ
values in all cases.

Since an improvement in the wavefunction should correspond to an improvement
in accuracy for all quantum mechanical properties, we further tested the
changes in HOMO-LUMO gaps and dipole moments. This is displayed in Tab. 3 for
the same set of out-of-domain molecules. All values are compared to the larger
and more balanced cc-pVQZ basis set8 since dipole moment values are known to
be highly sensitive to basis set size32. All adaptive basis sets here as well
were predicted by the ML model trained on the first 2000 smallest QM9
molecules. As extensively discussed in prior studies, the form and choice of
the optimal basis set can vary significantly for different molecular
properties33; 34. For both HOMO-LUMO gaps and dipole moments the errors do not
consistently decrease with increasing basis set size, however the adaptive
version improves the errors in all cases except the HOMO-LUMO gaps for a3-21G.
The dipole moment norms with a3-21G however do show an improvement indicating
more accurate electron densities which is seen with all 3 adaptive versions.
Improvement of multiple properties that were not part of training highlights
the ”foundational model” aspect35 of this method for basis sets. This can also
be exploited for the generation of accurate quantum chemistry datasets for ML
while retaining the cost of a cheaper basis set similar to how the revQM9
dataset was generated using an adaptive exchange-correlation density
functional35.

###  III.2 Transferability

Figure 5: Exemplary molecules from the QM9 dataset with diverse functional
groups and their optimal (HF energy-minimized) STO-3G valence orbital scaling
factors $\zeta_{\text{opt}}$.

The ease with which the orbital scaling factors can be learned, as indicated
by our numerical results, suggests a strong mapping between them and an atom’s
local environment. This is desirable since the ML model can be trained on a
small dataset of diverse local chemical environments and still remain
transferable and applicable to larger compounds. Fig. 5, shows scaling factors
for 14 molecules from the QM9 dataset with diverse functional groups and
different orbital hybridizations for the STO-3G basis set. Upon examination,
$\zeta$ appears to be a very local quantity, for the most part influenced only
by its direct bonding partners. This is consistent with the local and short-
ranged nature of atomic orbitals and underlines the flexibility of this
approach. This strong locality allows partitioning large molecules with
multiple functional groups into smaller fragments (or amons 29) with, for the
most part, independent scaling factors. This is likely the reason why our ML
model seems to be transferable to larger molecules.
In previous work, $\Delta$-36 and multilevel37; 26 machine learning methods
have been used to predict corrections to calculations performed in small basis
sets with significant improvements. These improvements are even larger since
molecular properties are learned from reference data using a larger basis set
however separate models and labels are required for each property. Adaptive
basis sets, on the other hand, are a foundational ML model thereby providing
improved wavefunctions and multiple improved properties/observables at once
since the method remains fully ab initio. The training data acquisition cost
is also lower since labels are generated via application of the variational
principle to the same small basis set without need for more expensive, larger
basis set reference calculations. Furthermore, the method might still be more
transferable than $\Delta$-/multilevel methods due to the strong locality of
the scaling factors depending on the functional groups. Such comparisons will
be analyzed in a future work employing a larger and more diverse training
dataset while extending the method to basis sets from other families.

##  IV Conclusion

In this work have we have presented a ML model for adaptive basis sets which
improves Hartree-Fock wavefunctions at negligible additional cost. The method
predicts optimally scaled radial exponents of contracted Gaussian type
orbitals for each atom within a molecule depending on its local chemical
environment. Optimal scaling factors for the STO-3G, 3-21G, 6-31G and 6-31G*
Pople basis sets were generated via variational minimization of Hartree-Fock
total energies. Application of the variational principle provides training
data without need for higher level reference calculations and guarantees
simultaneous improvement of all observables due to more accurate ground state
wavefunctions. Testing on 30,000 random molecules from the QM9 dataset13
showed that the default basis sets could be improved upon already after
training on only $\sim$100 molecules by using a local, atomic regression model
(cf. Fig.3). After increasing the training set size to 2000 molecules, HF
total and atomization energies systematically improve and the accuracy
increase was found to be particularly large for 3-21G (cf. Tab 1).

Making use of the high locality of $\zeta$, we investigated a variety of
biologically/societally relevant molecules, most of which are significantly
larger than the molecules used during training. The model is found to be
transferable, as shown by the observed improvements in energetics, HOMO-LUMO
gaps and dipole moment norms for these larger molecules further highlighting a
foundational model aspect (cf. Tabs. 2,3). Analysis of the scaling factors for
STO-3G indicates their strongly local character as similar functional groups
and atoms with the same bonding partners often evoke similar optimal scaling
factors (cf. Fig. 5). This supports the locality hypothesis which underpins
the transferability of the method which is desirable in order to efficiently
cover diverse chemistries through the use of suitable partitioning schemes,
e.g. amons29 or GEMS 38. Alternative methods for parameterizing $\zeta$ could
be possible, however due to the increasing complexity of the problem for
larger basis sets, we have opted for simple kernel ridge regression based
machine learning with a compact local representation, such as many-body
distribution functionals (MBDF).

While this work has focused on HF in order to exploit its variational nature,
scaled basis sets are expected to offer improvements to post-HF methods (such
as CCSD(T) which is known to converge slowly with basis set size) building
upon the improved ground state HF wavefunction, and also Kohn-Sham density
functional theory based methods employing the same self-consistent-field
framework as HF. Since the size of the basis set remains constant within the
adaptive basis set framework, improvements in wavefunction and properties of
test molecules are obtained for negligible computational overhead, improving
the computational cost vs accuracy trade-off for practically any post-Hartree-
Fock quantum chemistry calculation (cf. Tab. 4). Regarding density
functionals, note that adaptive atomic basis sets can easily be paired with
other adaptive parameters in the exchange correlation potential, such as the
adaptive hybrid density functional approximation35.

Finally, while this work has been limited to the smallest Pople basis set
family, the approach is expected to be equally applicable to other popular
basis set families in quantum chemistry, since they are also based on similar
contracted Gaussian type orbitals. It is also straightforward to incorporate
other chemical elements and to atomic orbital based treatment of condensed
phase systems through programs, such as FHI-AIMS39.

##  V Data and code availability

Generated training data, scripts and a trained model for obtaining optimal
scaling factors for the basis sets discuessd in this work are available at
https://github.com/dkhan42/aBasis. Script for generating MBDF along with
gradients is available at https://github.com/dkhan42/MBDF.

##  VI Acknowledgments

D.K. acknowledges discussions with G. Domenichini. We acknowledge the support
of the Natural Sciences and Engineering Research Council of Canada (NSERC),
[funding reference number RGPIN-2023-04853]. Cette recherche a été financée
par le Conseil de recherches en sciences naturelles et en génie du Canada
(CRSNG), [numéro de référence RGPIN-2023-04853]. This research was undertaken
thanks in part to funding provided to the University of Toronto’s Acceleration
Consortium from the Canada First Research Excellence Fund, grant number:
CFREF-2022-00042. O.A.v.L. has received support as the Ed Clark Chair of
Advanced Materials and as a Canada CIFAR AI Chair. O.A.v.L. has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No. 772834).

## References

  * Jensen (2013) F. Jensen, WIREs Computational Molecular Science 3, 273 (2013), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1123 .
  * Hehre  _et al._ (2003a) W. J. Hehre, R. Ditchfield, R. F. Stewart, and J. A. Pople, The Journal of Chemical Physics 52, 2769 (2003a), https://pubs.aip.org/aip/jcp/article-pdf/52/5/2769/11347023/2769_1_online.pdf .
  * Ditchfield  _et al._ (2003) R. Ditchfield, W. J. Hehre, and J. A. Pople, The Journal of Chemical Physics 54, 724 (2003), https://pubs.aip.org/aip/jcp/article-pdf/54/2/724/11217135/724_1_online.pdf .
  * Hill (2012) J. G. Hill, International Journal of Quantum Chemistry 113, 21–34 (2012).
  * Davidson and Feller (1986) E. R. Davidson and D. Feller, Chemical Reviews 86, 681 (1986), https://doi.org/10.1021/cr00074a002 .
  * Binkley  _et al._ (1980) J. S. Binkley, J. A. Pople, and W. J. Hehre, Journal of the American Chemical Society 102, 939–947 (1980).
  * Hehre  _et al._ (2003b) W. J. Hehre, R. Ditchfield, and J. A. Pople, The Journal of Chemical Physics 56, 2257 (2003b), https://pubs.aip.org/aip/jcp/article-pdf/56/5/2257/11046672/2257_1_online.pdf .
  * Dunning (1989) J. Dunning, Thom H., The Journal of Chemical Physics 90, 1007 (1989), https://pubs.aip.org/aip/jcp/article-pdf/90/2/1007/15358102/1007_1_online.pdf .
  * Schäfer  _et al._ (1992) A. Schäfer, H. Horn, and R. Ahlrichs, The Journal of Chemical Physics 97, 2571 (1992), https://pubs.aip.org/aip/jcp/article-pdf/97/4/2571/11359554/2571_1_online.pdf .
  * Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund,  _Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory_ , 1st ed. (Dover Publications, Inc., Mineola, 1996).
  * Khan  _et al._ (2024a) D. Khan, A. J. A. Price, M. L. Ach, and O. A. von Lilienfeld, “Adaptive hybrid density functionals,” (2024a), arXiv:2402.14793 [physics.chem-ph] .
  * Petersson  _et al._ (1988) G. A. Petersson, A. Bennett, T. G. Tensfeldt, M. A. Al‐Laham, W. A. Shirley, and J. Mantzaris, The Journal of Chemical Physics 89, 2193 (1988), https://pubs.aip.org/aip/jcp/article-pdf/89/4/2193/11255285/2193_1_online.pdf .
  * Ramakrishnan  _et al._ (2014) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, Scientific Data 1 (2014), https://doi.org/10.1038/sdata.2014.22.
  * Behler (2011) J. Behler, The Journal of Chemical Physics 134, 074106 (2011), https://doi.org/10.1063/1.3553717 .
  * Christensen  _et al._ (2019) A. S. Christensen, F. A. Faber, and O. A. von Lilienfeld, The Journal of Chemical Physics 150, 064105 (2019), https://doi.org/10.1063/1.5053562 .
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal, Mathematical Programming 45, 503 (1989).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Nature methods 17, 261 (2020).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Computational Molecular Science 8, e1340 (2018), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1340 .
  * Sun  _et al._ (2020) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, The Journal of Chemical Physics 153, 024109 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0006074/16722275/024109_1_online.pdf .
  * Pritchard  _et al._ (2019) B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson, and T. L. Windus, Journal of Chemical Information and Modeling 59, 4814 (2019), pMID: 31600445, https://doi.org/10.1021/acs.jcim.9b00725 .
  * Vapnik (1999) V. Vapnik,  _The nature of statistical learning theory_ (Springer science & business media, 1999).
  * (22) C. E. Rasmussen _et al._ ,  _Gaussian processes for machine learning_ , Vol. 1 (Springer).
  * Deringer  _et al._ (2021) V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, Chemical Reviews 121, 10073 (2021), pMID: 34398616, https://doi.org/10.1021/acs.chemrev.1c00022 .
  * Faber  _et al._ (2017) F. A. Faber, L. Hutchison, B. Huang, J. Gilmer, S. S. Schoenholz, G. E. Dahl, O. Vinyals, S. Kearnes, P. F. Riley, and O. A. von Lilienfeld, Journal of Chemical Theory and Computation 13, 5255 (2017), pMID: 28926232, https://doi.org/10.1021/acs.jctc.7b00577 .
  * Huang  _et al._ (2023) B. Huang, O. A. von Lilienfeld, J. T. Krogel, and A. Benali, Journal of Chemical Theory and Computation 19, 1711 (2023), pMID: 36857531, https://doi.org/10.1021/acs.jctc.2c01058 .
  * Heinen  _et al._ (2023) S. Heinen, D. Khan, G. F. von Rudorff, K. Karandashev, D. J. A. Arrieta, A. J. Price, S. Nandi, A. Bhowmik, K. Hermansson, and O. A. von Lilienfeld, arXiv preprint arXiv:2308.11196 (2023), 2308.11196 .
  * Lee  _et al._ (2023) S. Lee, S. Heinen, D. Khan, and O. A. von Lilienfeld, arXiv preprint arXiv:2308.00389 (2023), 2308.00389 .
  * Khan  _et al._ (2023) D. Khan, S. Heinen, and O. A. von Lilienfeld, The Journal of Chemical Physics 159, 034106 (2023), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0152215/18048072/034106_1_5.0152215.pdf .
  * Huang and von Lilienfeld (2020) B. Huang and O. A. von Lilienfeld, Nature Chemistry 12, 945 (2020).
  * Kwon  _et al._ (2023) H.-Y. Kwon, G. M. Curtin, Z. Morrow, C. T. Kelley, and E. Jakubikova, International Journal of Quantum Chemistry 123, e27123 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.27123 .
  * Neese (2022) F. Neese, WIREs Computational Molecular Science 12, e1606 (2022).
  * Hickey and Rowley (2014) A. L. Hickey and C. N. Rowley, The Journal of Physical Chemistry A 118, 3678 (2014).
  * García  _et al._ (2021) V. García, D. Zorrilla, M. Fernández, and J. Sánchez-Márquez, Journal of Computational Chemistry 42, 2154 (2021), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.26745 .
  * Rappoport and Furche (2010) D. Rappoport and F. Furche, The Journal of Chemical Physics 133, 134105 (2010), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3484283/15431477/134105_1_online.pdf .
  * Khan  _et al._ (2024b) D. Khan, A. J. A. Price, M. L. Ach, O. Trottier, and O. A. von Lilienfeld, “Adaptive hybrid density functionals,” (2024b), arXiv:2402.14793 [physics.chem-ph] .
  * Ramakrishnan  _et al._ (2015) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, Journal of Chemical Theory and Computation 11, 2087 (2015), pMID: 26574412, https://doi.org/10.1021/acs.jctc.5b00099 .
  * Zaspel  _et al._ (2019) P. Zaspel, B. Huang, H. Harbrecht, and O. A. von Lilienfeld, Journal of Chemical Theory and Computation 15, 1546 (2019), pMID: 30516999, https://doi.org/10.1021/acs.jctc.8b00832 .
  * Unke  _et al._ (2024) O. T. Unke, M. Stöhr, S. Ganscha, T. Unterthiner, H. Maennel, S. Kashubin, D. Ahlin, M. Gastegger, L. Medrano Sandonas, J. T. Berryman, _et al._ , Science Advances 10, eadn4397 (2024).
  * Blum  _et al._ (2009) V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009).