RESEARCH ARTICLE | JUNE 21 2024

A polarizable valence electron density based force field forhigh-energy interactions between atoms and molecules
Special Collection: 2024 JCP Emerging Investigators Special Collection

José Romero; Paulo Limão-Vieira; Thana Maihom; Kersti Hermansson; Michael Probst 

J. Chem. Phys. 160, 235101 (2024)
https://doi.org/10.1063/5.0210949


View
OnlineExportCitation

09 August 2024 10:12:59

The Journalof Chemical Physics

ARTICLEpubs.aip.org/aip/jcp

A polarizable valence electron density based force
field for high-energy interactions between atoms
and molecules
Cite as: J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949
Submitted: 28 March 2024 • Accepted: 29 May 2024 •Published Online: 21 June 2024
José Romero,1 , 2,a)Paulo Limão-Vieira,2
1 , 3,b)and Michael Probst

Thana Maihom,3 , 4

Kersti Hermansson,5

AFFILIATIONS
1Institute of Ion Physics and Applied Physics, University of Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria
Atomic and Molecular Collisions Laboratory, CEFITEC, Department of Physics, Universidade NOVA de Lisboa,
2829-516 Caparica, Portugal
3
School of Molecular Science and Engineering, Vidyasirimedhi Institute of Science and Technology, Rayong 21210, Thailand
4
Department of Chemistry, Faculty of Liberal Arts and Science, Kasetsart University, Kamphaeng Saen Campus,
Nakhon Pathom 73140, Thailand
5
Department of Chemistry-Ångström, Uppsala University, Box 538, SE-75121 Uppsala, Sweden
2

09 August 2024 10:12:59

Note: This paper is part of the 2024 JCP Emerging Investigators Special Collection.
a)Electronic mail: j.romero@alumni.fct.unl.ptb)Author to whom correspondence should be addressed: michael.probst@uibk.ac.at

ABSTRACT

High-accuracy molecular force field models suited for hot gases and plasmas are not as abundant as those geared toward ambient pressure andtemperature conditions. Here, we present an improved version of our previous electron-density based force field model that can now accountfor polarization effects by adjusting the atomic valence electron contributions to match ab initio calculated Mulliken partial charges. Using aslightly modified version of the Hohenberg–Kohn theorem, we also include an improved theoretical formulation of our model when applied
to systems with degenerate ground states. We present two variants of our polarizable model, fitted from ab initio reference data calculated at
CCSD(T)/cc-pVTZ and CCSD(T)/CEP-31G levels of theory, that both accurately model water dimer interaction energies. Further improvements include the additional interaction components with fictitious non-spherically symmetric, yet atom-centered, electron densities andfitting the exchange and correlation coefficients against analytical expressions. The latter removes all unphysical oscillations that are observedin the previous non-polarizable variant of our force field.
© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial 4.0
International (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/). https://doi.org/10.1063/5.0210949

I. INTRODUCTION
In computational chemistry and physics, the rationale behind
employing force fields stems from the ability to “integrate out” theelectronic degrees of freedom. One can thereby avoid solving thehigh-dimensional partial differential Schrödinger equation in favorof a tractable many-body problem of coupled second-order ordinary differential equations (ODEs) embodied in Newton’s laws ofmotion.

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

Force fields are typically formulated as sums of inverse powerlaws, with the Lennard-Jones 12-6 potential form standing out asa prominent example. Some examples of force field models thatare often used in Molecular Dynamics (MD) simulations and that
rely on such inverse power laws include Jorgensen’s Transferable Intermolecular Potential (TIP) family of water models (e.g.,
TIP3P and TIP4P),1–4 Optimized Potential for Liquid Simulations(OPLS),5–7 and Assisted Model Building and Energy Refinement(AMBER)8,9 family of force field models. Many other examples of160, 235101-1

The Journalof Chemical Physics

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

pubs.aip.org/aip/jcp

calculated with Kohn-Sham DFT (KS-DFT), these electron densities
are described as fictitious as they solely are convenient mathematical abstractions that lack a rigorous physical interpretation. Using
Hohenberg and Kohn’s theorem,32 it is demonstrated that for nondegenerate cases, mapping the non-interacting electron densities ofthe molecules to the energy of the interacting system is possible. This
argument extends to systems with degenerate ground states in their
non-interacting components through Levy’s constrained search formulation.33 The main advantage of this approach is that besides
requiring the usual molecular specification needed by any, full atom,
force field model, e.g., the Z-matrix, it only needs partial atomic
charges, e.g., Mulliken or Löwdin charges, of every molecule in thesystem (when isolated) once a set of exchange and correlation coefficients are fitted against experimental or accurate ab initio reference
data. However, in our previous paper, these partial charges weretreated as constants, which, therefore, made the previous version of
our model non-polarizable, and potentially too inaccurate for somesystems.
In Sec. II, we highlight three main limitations in our previous
model. Those limitations are mostly a consequence of the nonpolarizability of the model, and we provide theoretical means to
overcome them. In particular, as a theoretical foundation in thisimproved version of our model, we use a slightly modified, yetequivalent, version of a proof that deals with quantum mechanicalsystems with degenerate ground states34 and was previously derived
by Katriel et al. We favor our version of this proof as it is consistent
with the style of notation adopted in the theoretical derivation of ourprevious model. Here in Sec. II, we also provide a detailed description of a minimization problem that allows us to calculate atomic
contributions of the atomic electron density contributions that taketheir surroundings into account. This minimization problem, likesimilar others found in the literature,35,36 is interpreted as a selfenergy minimization of the electron clouds. Moreover, like many of
those other methods, the minimization problem we propose simply
involves a quadratic function, which translates to having to solve asystem of linear equations, therefore eliminating any convergence
related issues. Conservation of charge is guaranteed with a Lagrange
multiplier.
In Sec. III describes the systems used as reference and validation
datasets and the software used to calculate both datasets. We also
describe all the other well-established methods and the respective
software that we use to benchmark the accuracy of our model.
In Sec. IV presents preliminary findings for both calibration and validation datasets, employing our model with an orderten approximation alongside several other established force fieldmodels.
Finally, in Sec. V, we summarize our work, add final remarks,
and mention potential avenues for future developments of ourmodel.
II. THEORY AND PROBLEM FORMULATIONThis section is divided into three parts. Each of these parts
focuses on improving different limitations of our previous force
field model, limitations that, despite the model’s merits, remained
either unresolved or incompletely answered. In summary, in thepresent work, we target the following weaknesses of our previous
model:

160, 235101-2

09 August 2024 10:12:59

different empirical formulas can be found in the literature, such as
Halgren’s 14-7 buffered potential,10 featured in the Merck Molecular Force Field (MMFF).11–15 A standard issue associated with force
fields is the lack of transferability of optimized parameters between
models. In contrast to ab initio calculations, where the option
exists to combine different levels of theory for different molecules
and/or for parts of a same molecule(s) of a quantum mechanicalsetup,16–22 such flexibility is not always feasible for force fields. Furthermore, unless readily available, the force field parameterization
for a set of molecules must be fitted against reliable experimental
data and/or accurate ab initio calculations. This task is often very
time-consuming and can be discouraging for scientists with limited
experience in parameterizing molecular force field models.
Efforts have been made in the literature to connect quantum
chemical calculations and force field expressions on a more fundamental level. For example, Bond Order Potentials (BOPs) canbe viewed as expressions closely related to the realm of density
functional (or wave function based) calculations.23,24 In addition
to providing physical clarity, these approaches avoid a proliferation of fitting parameters. However, despite their appealing nature,these approaches are not widely used, except for the Tight Binding (TB) method, which essentially also falls in this category. Much
more widely used at present are machine learning approaches. Thesemethods typically assume that the total energy of a system is thesum of the energy of all atoms in their respective environments. Asa first step, the Cartesian coordinates of this environment are converted into symmetry-adapted descriptors. In the two most popular
approaches, these descriptors are inputs to a neural network25 or to a
Gaussian approximation potential.26 The parameters of these models are obtained by fitting to the results of Density Functional Theory(DFT) calculations. The main advantage of the machine-learning
approach is the flexibility of the potential energy expression and its
learning-on-the-fly potential.
Another strategy for constructing force fields that has gained
popularity involves approximating electronic densities inside integral expressions that resemble Coulomb interaction integrals. Suchan approach offers several advantages. It overcomes the problems
arising from interaction energy formulas directly based on internuclear distances by instead using electron density distributions,which are more physical in nature. Such an approach also exploits
the experience gained from density functional calculations over theyears, and, consequently, requires fewer fitting parameters. Thiskind of approach often relies on Gaussian and/or Slater types ofbasis-functions to approximate the electron clouds. Exact expressions for these integrals exist when such basis functions are used.
Furthermore, compared to machine learned potentials, empiricalintuition can be built in more easily since their resulting expressions
strongly resemble the solution(s) of simple classical and/or quantum
mechanical problems. Notable examples include the Gaussian Electrostatic Model (GEM),27–29 the Monomer Electron Density Force
Field (MEDFF),30 and our own force field model,31 which we intend
to improve here.
One unique aspect of our model is how it calculates interaction energies within a system. We achieve this by repeatedly
applying the Laplace operator to the non-interacting electron densities of the system’s molecules. This process results in a finite
sum of fictitious non-interacting electron densities, for which exact
expressions can be derived. In the same fashion as the orbitals

ARTICLEThe Journalof Chemical Physics

1.

2.

3.

ARTICLEIts theoretical foundation could only establish unique maps
between non-interacting and interacting electron densities if all systems (interacting and non-interacting) are
non-degenerate.
All the contributions for the fictitious electron density inside
the integrals of the exchange and correlation functionals
were limited to atom-centered functions that have sphericalsymmetry.
The net atomic contributions to the electron densities
were scaled to match the Mulliken charges of the isolatedmonomers regardless of their surroundings in the system(hence non-polarizable).

In the continuation of this section, we will address, point-by-point,
each of those weaknesses by providing explicit formulae and/or
proofs wherever needed and/or applicable. Thus, subsection A below
refers to point 1 above, B to point 2, and C to point 3.
A. Hohenberg–Kohn theorem and degenerate
ground states∣ψk (r, r2 , . . . , rN )∣2
dr2 ⋅ ⋅ ⋅ drN .
K/N
R k=1
K

ρε0 (r) = ∫ 3 ⋅ ⋅ ⋅ ∫ 3 ∑R(1)From Lemma 1, we establish that for every eigenspace of any
Hamiltonian operator, there is a unique eigendensity associated with
it. This eigendensity can be easily computed when an orthonormal
basis of the corresponding eigenspace is known. Consequently, theexistence of the map f ∶= Ĥ ↦ ρε0 is affirmed. In addition, the reciprocal mappings of any Hamiltonian operator to its external potential
and vice versa are facilitated by the following functions:g(Ĥ) = Ĥ − T̂ − Û EE ,(2a)h(Û Ext ) = T̂ + Û EE + Û Ext ,(2b)

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

FIG. 1. Commutative diagram depicting the maps we have shown to exist of any
molecular system between its Hamiltonian, its external potential, and its groundstate eigendensity.

and finally, from Theorem 1, we know that the map v ∶= ρε0 ↦ Û Ext
exists. All maps just described are summarized in the diagram shownin Fig. 1.
We assume a pair of molecules A and B, with respective groundstate energies εA0 and εB0 and respective non-interacting groundstate eigendensities ρεA0 and ρεB0 . Each eigendensity can be mapped
to its respective non-interacting Hamiltonian (Ĥ A and Ĥ B ), and,
therefore, the map
ρεA0 ⊗ ρεB0 ↦ Ĥ A ⊗ Ĥ B

(3)

must also exist. Since the Hamiltonian of each molecule can bemapped to the position and atomic number of all atoms, the noninteracting Hamiltonians can be mapped to the corresponding
interacting Hamiltonian, that is, the map
ρεA0 ⊗ ρεB0 ↦ Ĥ AB(4)

must exist; therefore, the non-interacting ground state eigendensities of molecules A and B can be mapped to the eigendensity of theinteracting system; in other words, the map
ρεA0 ⊗ ρεB0 ↦ ρεAB0(5)must also exist. This argument can be used on a per atom basis,
which serves to justify separating the entire interaction energy on a
per atom pair sum of contributions. Mapping the non-interacting
ground state eigendensities of each atom in our system to theground state eigendensity of the whole is advantageous since itallows us to exploit the spherical symmetry of their ground state
eigendensities. However, just like when Hohenberg and Kohn’s original proof is used instead, this generalization only holds for the
ground state eigendensities of both interacting and non-interacting
systems.
B. Atom-centered non-spherical electron
density contributions
Going back to our prior study, we introduced a straightforward
model that, through the adjustment of a finite set of parameters,
aspired to serve as a universal model for the interaction energies
between two interacting molecules. Although those aspirations were
not exactly realized, our previous model exhibited promising results,

160, 235101-3

09 August 2024 10:12:59

Using the Hohenberg–Kohn theorem, one can prove that there
exists a map between the non-interacting electron densities of two
molecules to the energy of the interacting system. With this knowledge, in our prior model, we aimed to, empirically, find an expression for this map in the form of a sum of fictitious Coulomb-like
integral expressions. However, this argument fails in justifying theexistence of a map between the non-interacting electron densities
and the interacting electron density due to ambiguity issues should
the interacting system allow degenerate ground states. Such kind of
map is useful to serve as a theoretical foundation of a polarizable
model, that is, a model whose interacting electron density, fictitious
or not, depends on the surroundings of the molecules of a system.
In this subsection, we provide theoretical means to overcome thisissue.
Let us begin by considering a molecule with N electrons in its
neutral state. The Hamiltonian operator governing this molecular
is denoted as Ĥ : L2 (R3N ) → L2 (R3N ), which we assume to be selfadjoint. The ground state energy, represented by the smallest eigenvalue of Ĥ, is denoted as ε0 ∈ R. The degeneracy of the ground stateis denoted as K, and let the wavefunctions ψ1 , . . . , ψK ∈ L2 (R3N ) collectively form an orthonormal basis for the lowest energy eigenspaceof Ĥ. We define the eigendensity, denoted as ρε0 ∈ L1 (R3 ), for thiseigenspace as the following integral:pubs.aip.org/aip/jcp

The Journalof Chemical Physics

ARTICLE

particularly for neutral systems.31 However, it lacks any atom-wise
contributions in the fictitious electron densities, refer to Eqs. (8a)and (8b), with no spherical symmetry, e.g., dipole-like contributions.
The interaction energy between two molecules, A and B, as definedby our simple model, is given by
EE
EN
NN
ΔEAB = ΔEAB
+ ΔEAB+ ΔEAB
K

EE(k)
EN(k)
+ ∑ (ak XCAB + bk XCAB ),
k=0

(6)

in which the naïve model contributions are
EE
ΔEAB=∫ 3∫ 3R

ρεA0 (r)ρεB0 (r′ )∥r − r′ ∥R

dr′ dr,A
B⎛ MA Zm ρεB0 (r) MB Zm ρεA0 (r) ⎞EN
ΔEAB= −∫ 3 ∑
dr,
+
∑
A
B
R ⎝m=1 ∥r − Rm ∥m=1 ∥r − Rm ∥ ⎠(7a)

To address this, we add to our model the following pair of family of functionals of the non-interacting ground state eigendensities
of A and B:ED(k): W 2k+1,1 (R3 ) × W 2k+1,1 (R3 ) → R,

(9a)

ND(k)

: W 2k+1,1 (R3 ) × W 2k+1,1 (R3 ) → R,

(9b)

XCAB

XCAB

wherein the electron clouds and nuclei are made to interact with fictitious atom-centered electron density contributions with cylindricalsymmetry, which are also conveniently parameterized by the groundstate eigendensities of A and B, whose exact expressions we assume
to be
ED(k)

XCAB
(7b)

2k+1
′
⎛ ρεA0 (r)rAB ∇r′ ρεB0 (r )1⋅
= ∫ 3∫ 3
2 R R ⎝ ∥rAB ∥∥r − r′ ∥

+
A B
Zm
Zm′
NN
ΔEAB
=∑ ∑
A
B ,
m=1 m′ =1 ∥Rm − Rm′ ∥
MA

MB

ND(k)

XCAB

R

⎝

∥r − r ∥

EN(k)

XCAB

(8a)

(10b)

⎠

B 2k
A 2k
⎛ MA Zm ∇ ρεB0 (r) MB Zm ∇ ρεA0 (r) ⎞+dr, (8b)=∫ 3 ∑∑
R ⎝m=1∥r − RAm ∥∥r − RBm ∥ ⎠m=1

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949

2k+1
BρεA0 (r) ⎞
Zm
rBA ∇
⋅
dr,
∥r
∥
∥r
− RBm ∥ ⎠
BA
m=1
MB

+ ∑

A
where RAm and Zmlabel the coordinates and atomic number of theB
mth atom of A, RBm and Zmlabel the mth atom of B, and the operator ∇2k is shorthand notation for k consecutive applications of the
Laplace operator. Despite its simplicity, in our previous study, we
found this model to reproduce water dimer interactions to a fairdegree of accuracy, albeit with a rather high approximation order(six or higher).
One can easily see from the exchange and correlation functionals above that the heuristic philosophy driving this model isthe introduction of fictitious electron clouds by means of consecutive applications of the Laplace operator on the non-interacting
atomic density contributions, thus making each of these contributions spherically symmetric by construction. This is likely a limiting
factor in our model, since capturing asymmetries in the interacting
densities is likely to be more successfully achieved by the inclusion of
atom-centered contributions that have no spherical symmetry in thefictitious interacting densities. An example where this plays a significant role is a water molecule with its lone electron pair, which canclearly not be described by spherical densities residing on the oxygen
atom. This may also have contributed to the requirement of having to use rather high order expansions (six or higher) to accurately
model water dimer interactions in our previous model.

© Author(s) 2024

2k+1
ρεB0 (r)⎛ MA ZA rAB ∇⋅=∫ 3 ∑ m
R ⎝m=1 ∥r AB ∥∥r − RAm ∥∥r − r ∥′ρεB (r)∇2k
r′ ρεA0 (r ) ⎞ ′+ 0
dr dr,′(10a)in which ∇2k+1 is shorthand notation for one and k consecutive
applications of the gradient and Laplacian operators, respectively,
and the distance vectors rBA , rAB ∈ R3 arerAB = ∫ 3rρεB0 (r)NB

R

rBA = ∫ 3
RrρεA0 (r)NAdr − ∫ 3rρεA0 (r)dr − ∫ 3rρεB0 (r)
NB

R

dr,

(11a)

dr = −rAB.

(11b)

NA

R

One can, therefore, see that for the zeroth-order contributionsof Eqs. (10a) and (10b), that is, when k = 0, the fictitious electron densities obtained through differentiation are atom-centereddipoles, which resemble p-type atomic orbitals, oriented along the
distance vectors between molecules A and B. Finally, the interactionenergy between A and B in our improved model is
K

EE(k)

EE
EN
NN
ΔEAB = ΔEAB
+ ΔEAB+ ΔEAB+ ∑ (ak XCAB
k=0
EN(k)

+ bk XCABED(k)+ ck XCABND(k)+ dk XCAB

),

(12)

where ak , bk , ck , and dk denote the multiplying coefficients, 4(K + 1)
in total, which need to be fitted against reliable ab initio and/or
experimental data.

160, 235101-4

09 August 2024 10:12:59

2k′⎛ ρεA0 (r)∇r′ ρεB0 (r )1
EE(k)XCAB = ∫ 3 ∫ 3′R

ρεB0 (r)rBA ∇2k+1
ρεA0 (r′ ) ⎞ ′r′⋅
dr dr,∥rBA ∥∥r − r′ ∥ ⎠(7c)and the exchange and correlation contributions are2

pubs.aip.org/aip/jcp

The Journalof Chemical Physics

ARTICLE

C. Modeling Mulliken partial charges
The final limitation in our previous model that we intend tocorrect is its inability to polarize the molecules in a form that allows
for charge redistribution. With this, we mean that the atomic contributions to the total electron density can shift depending on theenvironment of the atoms. It is noted that our previous modelaccounts for dipole moments of the molecules to an extent, since
the atomic density contributions are multiplied by scaling a scaling
factor that ensures that these contributions agree with the Mulliken
partial charges obtained from high accuracy ab initio calculations.
This has the disadvantage of requiring ab initio calculations on every
molecular species in the system that one aims to model.
In this subsection, we improve our model by lifting this requirement by not only modeling interaction energies, but also the Mulliken charges. Even though Eqs. (10a) and (10b) can account forasymmetries in the interacting ground state eigendensities, they do
not allow the net atomic contributions to the density to change, akinto how in our previous version of our model the atomic density contributions are simply scaling factors that match the Mulliken chargesof an isolated monomer ab initio calculation.
Let us begin by considering a simple system with two molecules
A and B. We approximate the interacting ground state eigendensitiesof these molecules asMAρ̃εA0 (r) = ∑ ∑nMBρ̃εB0 (r) = ∑ ∑
m=1

n

3/2

B B
χm
cm,n

3/2

λAm,n
)
(
2
exp (λAm,n ∥r − RAm ∥ ) π
λB( m,n )B
B 2
exp (λm,n ∥r − Rm ∥ ) π

,

(13a)

,

(13b)

A B
where χm
, χm ∈ R serve as labels for unknown scalar variables
that allow the net atomic contributions to the eigendensity of themolecules to match the ideal ab initio calculated Mulliken partialA
B
charges, whereas the constants labeled by cm,n, cm,n, λAm,n , and λBm,n aretaken from our previously tabulated atom-specific values.31 A new
set of atom-specific coefficients, for the full electron variant of ourmodel, is calculated and presented in Table S1 of the supplementarymaterial, for all elements in the first three rows of the Periodic Table.
This new set of coefficients are optimized so that the last threeGaussian contributions are the same as their Effective Core Potential (ECP) counterparts, thus allowing us to separate between coreand valence electrons. We obtain the polarization variables from thefollowing minimization problem:⎧⎪
minimize⎪⎪⎪⎨⎪subject to
⎪
⎪
⎪
⎩

Pol A
A
B
(χ1 , . . . , χM)
UAB
, χ1B , . . . , χM
A
B
MA

A A

MB

A B

∑ ∑ χm cm,n + ∑ ∑ χm cm,n = NA + NB ,
m=1

n

m=1

(14)

n

where N A and N B are the number of electrons of molecules A and
Pol
B and UAB
: RMA +MB → R is a function, that is yet to be defined, ofthe polarization scalars associated with every atom in the system,
which we treat as optimization variables. Naturally, the equality constraint in Eq. (14) is imposed by using a Lagrange multiplier. Foranionic systems, N A + N B , in the equality constraint of Eq. (14), is
replaced with N A + N B + 1. In addition, for the full electron variantof our model, only the valence electron contributions are allowed toJ. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

change, that is, the inner core shell electron density contributionsof the atoms of each molecule are treated as frozen electron shells.
Mathematically, this translates to their, the core shell contributions,associated polarization scalars being constrained to always being
equal to one regardless of what surrounds them; hence, we do nottreat these as optimization variables but rather as mere constants.
Contrary to our previous model, where computing all
intramolecular interactions can be avoided since they cancel out
when calculating interaction energies, this is no longer the case if
we allow the molecules in our system to become polarized. Thisis a direct consequence of the intramolecular interaction energy
components being dependent on what surrounds every molecule inour system. Therefore, the energies of each molecule when isolated
need to be calculated, only once, and then subtracted from the totalenergy, that is, intra- and intermolecular interaction energies where
the latter is disregarded in our previous model due to it canceling asa consequence of its, inaccurate, non-polarizability.
While it can be proved that there exists a functional thatmaps the non-interacting neutral ground state eigendensity of all
the atoms in a polyatomic system to its interacting, not necessarily neutral, ground state eigendensity, which can be done using themodified Hohenberg–Kohn theorem, it does not provide an explicit
Pol
formula for the map. As an ansatz, we model UAB
as being a functional whose mathematical expression is identical to our exchange
and correlation corrected interaction energy [Eq. (12)], and whosepolarization coefficients minimize it, akin to a self-energy minimization. However, the physical interpretation of this functional
is dubious; thus, a set of different exchange and correlation coefficients are used to model the Mulliken charges. This separate set of
exchange and correlation coefficients are obtained from least squares
minimization against the Mulliken charges of the same dataset used
to calibrate the interaction energy model.
This method extends beyond simply modeling the Mulliken
charges of diatomic systems (Fig. 2), which is the typical focus instudies using molecular dynamics simulations. These simulations
often involve thousands of molecules.
III. METHODS
The current, improved version of our model is calibrated withthe same reference data that were used to calibrate the previous nonpolarizable version.31 In our previous study, the valence electron
densities of all the atoms, in the first three rows of the Periodic Table,
are fitted, and tabulated, via least-squares minimization against theelectronic density obtained from Density Functional Theory (DFT)
calculations using the B3LYP37 hybrid functional with Grimme’s
empirical dispersion formula GD3.38 The Effective Core Potential(ECP) basis-set CEP-31G39–41 were used for all atoms except hydrogen and helium where the full-electron basis-set cc-pVDZ42–45 isused instead. The water dimer interaction calculations used as reference data consist of a 25 × 200 uniform grid of calculations in
the O⋅ ⋅ ⋅O distance coordinate and a 25 × 200 uniform grid of calculations in the H⋅ ⋅ ⋅O distance coordinate. This means that theinteratomic distances are sampled with 200 points, whereas theangular component at which one of the molecules of water is rotatedin an axis perpendicular to the plane of the molecule is sampledwith 25 points, using CCSD(T)46–50 level of theory with the basis-set
cc-pVTZ.51

160, 235101-5

09 August 2024 10:12:59

m=1

A A
cm,n
χm

pubs.aip.org/aip/jcp

The Journalof Chemical Physics

ARTICLEpubs.aip.org/aip/jcp

To validate our model, comparisons between different energies
predicted by our model and accurate ab initio calculations are madeon a set of different complexes with small molecules, namely H2 ,N2 , O2 , and CO2 . These comparisons include the binding energies
between two atoms, for the diatomic molecules mentioned above.
To determine its ability to predict covalent bonds. Comparisons
between Vertical Electron Affinities (VEAs) of four small molecules
and interaction energies between dimer complexes of two identicalmolecules in their neutral and anion states are also provided. The
comparisons of the dimer complexes consist in displacing one of themolecules by a distance ΔL ∈ [1, 5], in Å, and rotating it by an angleθ ∈ {0, π/4, π/2}.
The bond lengths and angles used for all the small molecules
considered for this study, for both the calibration and validation, aresummarized in Table I.
All the calculations used to validate our model are obtained
from ROCCSD(T) and CCSD(T) calculations (for non-singlet andsinglet states, respectively), which use the basis-set cc-pVTZ. Moreover, we also compare the results in the validation dataset withthe predicted values from the already established force field models: General AMBER Force Field (GAFF),52 Merck Molecular Force
Field (MMFF94S),11–15 both of them using the open-source software Open Babel,53 Effective Fragment Potential (EFP)54–56 method
using the open-source software LIBEFP,57 and MB-nrg58,59 using
the open-source software MBX.60 All ab initio calculations that areused and presented for this study are obtained from the Gaussian

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

TABLE I. Bond lengths and angles used for all the small molecules studied in thisarticle. This list of molecules includes all the molecules in the calibration and validation
datasets.

Molecule
H2 O
CO2
H2
N2
O2

Type

Length (Å)

O–H
H–O–H
C–O
O–C–O
H–H
N–N
O–O

0.983

Angle (○ )
108.9

1.163
180.0
0.74
1.09
1.16

16 software61 using the LEO High Performance Computing (HPC)
facilities in the University of Innsbruck.
The fitting of the exchange and correlation coefficients used
to polarize the molecules as described in Eq. (14) is performed
with the following combination of methods: An initial optimization using the Nelder–Mead simplex algorithm62–64 using only
220 randomly picked water dimer configurations of our calibration dataset, followed by another minimization using a low memory
Broyden–Fletcher–Goldfarb–Shanno (LBFGS)65–71 algorithm, followed by a final minimization using the same LBFGS algorithm
against all water dimer configurations in our calibration dataset.

160, 235101-6

09 August 2024 10:12:59

FIG. 2. Diagram illustrating the energy calculation of a hypothetical system comprising three atoms, denoted as a, b, and c, whose respective non-interacting neutral atomic
ground state eigendensities are ρa , ρb , and ρc . Contrasting these, the atomic contributions to the modeled interacting eigendensity for the system are represented as ρ̃a , ρ̃b ,and ρ̃c . The overall charge of the system is denoted by Q, where for a neutral system, Q = 0. See text for explanation of the labels.

The Journalof Chemical Physics

ARTICLE

Exact values for the gradients for the minimization problem tocalibrate, and just the calibration, the XC coefficients for the polarization model are obtained using automatic (forward mode) differentiation in Julia using the open-source package ForwardDiff.72 For
the full electron variant of our method, only the Nelder–Mead algorithm is employed for the fitting of its parameters, due to numericalnoise issues stemming from huge numbers that arise when using thisvariant of our molecular force field model.
Unlike our non-polarizable model, where the XC coefficientsare obtained through an unconstrained least squares minimization
problem, in this study, we constrain them to follow exact alternating expressions to avoid both overfitting and unphysical oscillations
when extrapolating to shorter interatomic distances than those existing in our calibration dataset. The exact expressions used for every
XC coefficient, that is not zeroth order, are
A
,
(2k)!

(15a)

bk = (−1)k

B
,
(2k)!

(15b)

ck = (−1)kC,(1 + 2k)!

(15c)

dk = (−1)kD,(1 + 2k)!

(15d)where A, B, C, D ∈ R are optimized parameters (distinct for thepolarization and interaction energy models, and, hence, a totalof eight fitting parameters are used regardless of the order ofthe model) ∀n ∈ N+ /{0}. All zeroth-order coefficients are set tozero. Equations (15a)–(15d) are chosen due to these bearing astrong resemblance the coefficients previously obtained in a numerical experiment, whose values are not presented in this study, inwhich we perform least squares minimization against our calibration dataset with every coefficient being treated as an unconstrained
optimization variable.

IV. RESULTS AND DISCUSSION
Allowing any of the zeroth-order XC coefficients a0 , b0 , c0 , andd0 [refer to Eq. (12)] to be different than zero significantly compromises the accuracy of our model, especially for low-energy interaction configurations. However, setting these zeroth-order coefficientsto zero ensures that our model closely resembles the naïve model forlarge interatomic distances. This is a desirable feature since, at considerable distances between interacting molecules, the naïve model
approximates the exact answer due to the negligible correlation
between the electrons of each molecule. Moreover, by constraining
the other XC coefficients to exact expression as shown in Eq. (15a)through Eq. (15d), and then minimizing the coefficients of these
expressions with respect to our reference data, rather than having each XC coefficients as a minimization variable, eliminates
undesired oscillations at short interatomic distances.
Despite this, the ECP variant of our model noticeably underestimates the interaction energies, compared to accurate ab initio
calculations of a water dimer in the H⋅ ⋅ ⋅O distance coordinate
shown in the left panels of Fig. 3. This is likely due to the core electron contributions of the oxygen atom to the electron density not
being modeled. In contrast, the full electron variant of our model
accurately predicts these interaction energies; however, it displays
discernible numerical noise, despite 64-bit floating-point precision
being used throughout our entire code. Similar numerical instability

FIG. 3. Comparison between the interaction energies of full electron and ECP-fitted models, both at order ten, and the reference ab initio energies for a water dimer inhydrogen bonding (H⋅ ⋅ ⋅O; left panels) and oxygen facing (O⋅ ⋅ ⋅O; right panel) conformations. In the left panels, ΔL is the distance between O and H atoms, and in the rightpanel, ΔL is the distance between the O atoms.

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

160, 235101-7

09 August 2024 10:12:59

ak = (−1)kpubs.aip.org/aip/jcp

The Journalof Chemical Physics

pubs.aip.org/aip/jcp

of this study. Moreover, the inclusion of post-Born–Oppenheimer
corrections for high-energy configurations is another factor toconsider.
We analyzed the atomic partial charges by comparing the ones
obtained from our model and the reference ab initio dataset. Looking at the hydrogen bonded dimer configuration (left panels of
Fig. 4), a noticeable decline in accuracy is always observed when theH⋅ ⋅ ⋅O distances become shorter. All errors of the partial charges
in the all-electron variant are smaller than 0.15 e and only theerror of hydrogen bond acceptor oxygen number 4 is larger than
0.1 e, where e is the elementary charge. The partial charge errors
for the ECP variant are even much smaller, for every atom less than
0.07 e. The right panels of Fig. 4 provide an identical examination
for the “repulsive” water dimer conformation. Notwithstanding thecompletely different water–water orientation with the oxygen atoms
facing each other, the same conclusion concerning the superiority of
the ECP basis set variant and the overall accuracy is drawn.
In a broader analysis of the predicted partial charges from both
our ECP and full electron fitted variants of our model, comparedto their respective Mulliken charges from the SCF electron densities
computed at CCSD(T)/CEP-31G and CCSD(T)/cc-pVTZ levels oftheory, a consistent trend emerges. Both versions of our fitted model
showcase acceptable accuracy in estimating Mulliken charges, particularly for configurations with lower energy (see Fig. 5). Notably,
the ECP variant slightly outperforms in terms of accuracy, boasting a slightly larger coefficient of determination against its reference
ab initio data. This marginal improvement in accuracy is likely also
influenced by the inherent numerical error noise introduced due tothe use of limited floating-point precision representation and thenecessity to handle exponentiation of relatively large numbers, asevident in the decaying coefficients of the Gaussian basis functions

FIG. 4. Comparison of the errors in the partial charges from the two variants of our model (full-electron and ECP) for different water dimer geometries. The bars indicate theabsolute value of the differences with respect to Mulliken charges from ab initio calculations on the CCSD(T)/cc-pVTZ and CCSD(T)/CEP-31G levels of theory, respectively.
The two water dimer configurations used and the numbering of the atoms are shown in the sketches on top of the error bar plots.

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

160, 235101-8

09 August 2024 10:12:59

issues have also been reported in other studies that use comparable polarizable force field models.36 This numerical noise is due tothe high-order powers of the extinction/decay coefficients of theGaussian s-type functions, which are up to 20 for an order ten
approximation.
A close inspection of the water dimer displaced along theO⋅ ⋅ ⋅O distance coordinate shown in the right panel of Fig. 3 showsthat both the full electron and ECP variants of our model reproduce its respective interaction energies with a fair degree of accuracy.
Neither of the two variants of our model display unphysical oscillations through this distance coordinate. Like the H⋅ ⋅ ⋅O distancecoordinate, the full electron variant of our model displays the same
numerical noise in the O⋅ ⋅ ⋅O distance coordinate. Interestingly, forinteratomic distances between 1.5 and 2.0 Å, compared to the reference ab initio reference data, our full electron fitted model underestimates energies. Nonetheless, for interatomic distances shorter
than those illustrated in Fig. 3, the accuracy of our ECP fitted model
may diminish. This discrepancy arises from its inaccurate depiction of inner shell electrons, rendering it unsuited for extremely
high-energy molecular dynamics simulations where such configurations might be more prevalent. To enhance the accuracy of thefull electron variant of our model, for high-energy interactions, theincorporation of similar polarization coefficients into the core shell
contributions is needed, that is, no longer treating the inner shell
contributions as frozen electron clouds. Incorporating these innershell coefficients would require significant changes to the fitting of
the parameters used to model partial charges since relying solely
on Mulliken charges as reference would be insufficient. Instead,this would require the consideration of shell-specific contributions
to the electron density for each configuration within our calibration dataset. This more nuanced approach is beyond the scope

ARTICLEThe Journalof Chemical Physics

ARTICLEpubs.aip.org/aip/jcp

FIG. 5. Modeled atomic partial charges plotted against reference Mulliken charges, distinguished by the interaction energy of the respective geometry. Left panel: All-electron
variant; right panel: ECP variant. The orange and green data points refer to charges from configurations with interaction energies above or below +50 kJ/mol, respectively.
For clarity, the distribution of the data points along each axis is plotted above and to the right of the 2D plot panels.

In this fictitious scenario, a single electron cloud, characterized
by an amplitude equal to one and with decay coefficient λ = 3, itscenter is separated from the atom by a distance d. Figure 6 illustrates the convergence of the series, demonstrating that as the orderincreases, the sum converges to a defined limit. While we do not provide exact, nor approximate, expressions for the limit when K → ∞in this study, this topic still warrants to be mentioned given thatobtaining such an expression would contribute not only to reduce
the overall computational cost, but also to alleviate numerical noise,like it is observed in Fig. 3. Floating-point precision errors are prevented in Fig. 6 by obtaining exact values using the open-source
software Maxima,75 using its default arbitrary-precision arithmetic,
and only casting the result of each sum to a floating-point precision
representation after a factorization of the exact result.

A. Limit toward an infinite order expansion

To validate the order ten ECP variant of our model, for highinteraction energies, a set of different dimer complexes of two identical molecules, none of them being part of the calibration dataset,
is analyzed (Fig. 7). Specifically, the complexes H2 + H2 , N2 + N2 ,O2 + O2 , and CO2 + CO2 are considered, due to these being smallenough to allow high accuracy ab initio calculations to be made,and since the molecules that are part of these complexes are of considerable interest in various other scientific fields, such as biology,
atmospheric sciences, and materials science. Refer to Fig. S2 throughFig. S5 of the supplementary material for diagrams illustrating thegeometric configuration of the dimer complexes used in Fig. 7. For
the ab initio calculations of the interacting dimer complexes, all thesystems are assumed to be in their singlet states, except for all theO2 + O2 complexes where the systems are assumed to be quintets.
While these results are not shown in this article, identical results
are obtained for the same geometric configurations of O2 + O2 complexes in their triplet states if the interatomic separation between the
dimers is kept at not too long distances. For larger separations, theab initio calculations fail to converge with equal or less than 200 SCF
iterations, possibly due to basis-set limitations.
The ECP order ten variant of our model replicated the ab initio
calculated energies with a fair degree of accuracy, except for theN2 + N2 complexes. For the H2 + H2 complexes, the EFP results
displayed higher accuracy compared to our model, even at higher

Expressing our exchange and correlation contributions withcoefficients following the exact expressions in Eq. (15a) throughEq. (15d) establishes the legitimacy of treating this expansion asa series. This is evident in Fig. 6, where a normalized sum ofelectron–nuclei contributions is expanded up to orders 10 and
100 for a hypothetical system comprising an atom with an atomic
number of one.

FIG. 6. Comparison between normalized sums of electron–nuclei exchange andcorrelation functionals with 10 and 100 contributions of a fictitious system ofjust one atom and an electron cloud modeled by a single s-type Gaussianbasis-function with unit amplitude and decay coefficient λ = 3.

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

B. Validation against neutral complexes

160, 235101-9

09 August 2024 10:12:59

used to approximate atomic electron contributions (refer to Table S1of the supplementary material), namely for fitted full electron fitted
densities.
While the accuracy of our results is qualitatively comparable
to findings reported using other similar and well-established methods found in the literature, such as Atom-Bond Electronegativity
Equalization Method (ABEEM)73 and Electronegativity Equalization Method (EEM),74 it is noted that these results are confined toour calibration dataset. As such, validation against configurations
outside this dataset is crucial to affirm the broader applicability of
our model. Given the numerical instabilities inherent in our fullelectron model, arising from the challenges of handling exponentiation of large numbers, our validation is solely directed toward the
ECP variant of our model. Specifically, we validate both the interaction energy and polarization models against dimer configurations ofmolecules that are small enough to allow CCSD(T)/CEP-31G levelof theory to be used.

The Journalof Chemical Physics

ARTICLEpubs.aip.org/aip/jcp

FIG. 7. Comparison between the interaction energies of various neutral complexes obtained with an order ten ECP variant of our model against its ideal ab initio results. For
reference, comparisons with other established molecule force field models often found in the literature are also included. None of the dimer configurations shown in thesecomparisons are part of the calibration dataset.

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

configurations. Consequently, it proves inadequate for intramolecular geometry optimization. Therefore, to ensure accurate moleculargeometries, reliance on either a dependable database or theoretical calculations capable of optimizing molecular geometry is stillneeded.
Nonetheless, for systems comprising multiple molecules thatare treated as rigid objects, conducting geometry optimizations
remains feasible. Although explicit calculation or comparison offorces is lacking, in the dimer complexes, our data align with aunit slope in most cases. Thus, no scaling factor for the energies is required and it seems likely that also the forces derived
from our model are in good agreement with those obtained at the
CCSD(T)/cc-pVTZ level theory, provided that the molecules aretreated as rigid objects.
C. Validation against anionic complexes
For the same dimer complexes in their anion states, we findresults that somewhat resemble a linear trend between the interaction energies calculated with our order ten ECP variant of our modeland the equivalent total energies obtained from ROCCSD(T)/ccpVTZ level of theory as shown in Fig. 8. However, unlike the interaction energies of Fig. 7, the linear regressions obtained for all systems
show large intersects. No clear trend is observed between molecularsize nor atomic numbers, and the numerical value of these intersects.
Furthermore, the slope of these linear regressions varied roughly
between 0.4 and up to 2.0 depending on the system of dimers. As it
has already been mentioned, two likely reasons for the discrepancies
between our polarizable model and their ideal ab initio calculation counterparts are the narrowness of our calibration dataset,and our model’s inability to discern between atomic-shell specificcontributions to the electron density.

160, 235101-10

09 August 2024 10:12:59

energy regimes. Nevertheless, the ECP variant of our model displays a fair degree of accuracy for this complex for energies roughly
up to 50 kJ/mol. For the N2 + N2 complexes, both our model andMMFF94S displayed similar levels of accuracy, both differing to theideal ab initio result by a factor of ∼3.0 (Fig. 7). For the O2 + O2
and CO2 + CO2 complexes, our model displayed high level accuracy outperforming all the other molecular force fields consideredin the high energy regions. Notably, for the CO2 + CO2 complexes,
MB-nrg displayed a very close to linear relation against the ab initio
calculated data; however, a constant multiplicative factor of roughly
6.3 between the ab initio and MB-nrg calculated results is observed
as shown in the bottom right panel of Fig. 7. This factor is inagreement with the scaling factor of 0.15 for the two-body energies reported to be necessary for accurately simulating liquid andvapor phases.59 All the comparisons considered against other established force field models are constrained to instances wherein themolecules of the dimer complexes are encompassed within theirdefault fragment libraries, and/or those cases where the models do
not fail to run.
While our model demonstrates commendable accuracy in estimating interaction energies for high-energy neutral dimer systems,
its performance diminishes notably for lower energies (2 kJ/mol orless), as illustrated in Fig. S2 through Fig. S5 of the supplementarymaterial. In these instances, our model struggles to predict the formation of any weak bonds among neutral gases. Consequently, forsimulations involving low-energy scenarios, such as liquid-phasesimulations of H2 , N2 , O2 , and/or CO2 , this version of our force fieldmodel still proves to be inadequate and is likely to yield inaccurate
predictions for such systems.
Despite these improvements, our model, as depicted in Fig. S5of the supplementary material, still struggles to predict any covalent bonds, similar to its shortcomings in low energy interaction

The Journalof Chemical Physics

ARTICLEpubs.aip.org/aip/jcp

FIG. 8. Comparison between the total energies of various anion complexes obtained with an order ten ECP variant of our model against its ideal ab initio results. None of the
dimer configurations shown in these comparisons are part of the calibration dataset.

V. CONCLUSIONSWe present an upgraded version of our electron density-based
molecular force field, designed to replicate electrostatic polarization
effects. This improvement hinges on minimizing an energy functional derived from the non-interacting ground state eigendensitiesof all atoms within the system. Notably, this energy functional mirrors the proposed interaction energy, employing distinct optimized
parameters fitted to predict the Mulliken charges of a given system.
To validate the existence of such a functional, we furnish a proof
using a modified version of the Hohenberg–Kohn theorem, establishing a robust theoretical foundation for our polarizable molecularforce field.
In this enhanced version, we have implemented several refinements. These include incorporating electron density contributions

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

devoid of spherical symmetry, acquired from directional derivatives concerning the distance vector between pairs of atoms. In
addition, we have allowed the net atomic contributions to electron
density to vary based on the molecular surroundings. Analogous
to electron–electron and electron–nuclei exchange and correlation
contributions, these non-spherical contributions can be analyticallyresolved. Crucially, by constraining the fitted exchange and correlation coefficients to exact expressions instead of treating them as freeparameters, we have effectively eliminated unphysical oscillations in
predicted interaction energies in our numerical experiments. Thisenhancement is evident across both our full electron and ECP cases,
particularly when extrapolating to exceedingly small interatomic
distances.
Our upgraded model exhibits heightened flexibility and versatility since it can capture polarization effects. The flow of electron
density allows us to model atomic partial charges. Parameters for theelements from the first three rows of the Periodic Table are supplied.
This already allows for application to a wide range of molecularsystems. Furthermore, adapting the model for different ionic statesis straightforward, simply involving the adjustment of the equality
constraint of a Lagrange multiplier.
Notwithstanding the convenience and broad applicability, we
found that the accuracy of the presented model may occasionally deviate from ideal standards. This arises from various factors,
including the inherent complexities of molecular interactions andthe specificity of the calibration dataset.
Despite significant progress, certain aspects warrant further
investigation. For instance, exploring the incorporation of polarization effects to inner core electronic shells could be pivotal, particularly in high-energy molecular dynamics simulations. In addition,
efforts could be directed toward deriving exact or approximate
expressions when employing infinite expansions to mitigate numer-

160, 235101-11

09 August 2024 10:12:59

All the anion dimer complexes used in Fig. 8 are assumed
to be in their doublet states (one unpaired electron), except forthe (O2 + O2 )− complexes that are assumed to be in their quartet
states (three unpaired electrons). This makes sense as one of thetwo interacting O2 monomers will be in its triplet state for verylarge interatomic separations, while the other monomer is, in its
ground state, a doublet. While not presented, doublet state calculations for the (O2 + O2 )− complexes yield similar results, albeit withresults deviating slightly more from a linear regression. The geometric configurations used for the anion complexes are the same as
the ones used for the neutral complexes. For every panel of Fig. 8,
the respective collinear arrangement (θ = 0 for the dimer configurations shown in Fig. S2 through Fig. S5 of the supplementarymaterial) where the interatomic separation is at its maximum (5 Å)is used as the zeroth energy reference point for the interaction energy
calculations.

The Journalof Chemical Physics

ical noise associated with low-energy interactions. Nonetheless, our
ECP calibrated variant demonstrates negligible numerical noisefor up to a tenth-order truncated series expansion across various
dimer systems explored, yielding reasonably accurate predictions forboth interaction energies and partial charges in neutral low-energy
configurations.

SUPPLEMENTARY MATERIAL
In the supplementary material, a table with fitting parameters
for atomic electron densities in the full-electron variant of our model
can be found. The valence electron contribution is represented by
the last three Gaussian contributions for each element. The mathematical expressions for exchange and correlation functionals are also
given, together with the fitted coefficients. In addition, comparisonsof the energies predicted by our model with various small molecules
are provided.

ACKNOWLEDGMENTS

pubs.aip.org/aip/jcp

Probst: Investigation (equal); Validation (equal); Writing – review& editing (equal).
DATA AVAILABILITY
All the data and programming scripts used to make thismanuscript are freely available under a GPL-3.0 license andcan be accessed from the GitHub repository: https://github.com/
JoseRodriguezRomero/PolElectronDensityForceModel.
A software implementation intended for general used in theC++ programming language of this model for calculating point
energies and performing geometry optimization of molecular systems is also freely available under the GPL-3.0 license on GitHub
and Zenodo,76 along with documentation. The software relies on thefollowing C/C++ libraries: Eigen77 for efficient linear algebra operations, libLBFGS78 an implementation of J. Nocedal’s low-memory
BFGS algorithm69–71 in the C programming language, used in all
geometry optimizations, and autodiff79 for forward mode automaticdifferentiation, guaranteeing the calculation of accurate gradients ofthe energy.
APPENDIX A: HOHENBERG–KOHN THEOREM
AND GROUND STATE ELECTRON EIGENDENSITIES
Lemma 1. The eigendensity of any eigenspace associated with a
Hamiltonian operator remains invariant under any rotations applied
to the orthonormal basis being used.
Proof. Let the wavefunctions ψ 1 , . . . , ψ K form an orthonormal
basis of the eigenspace of Ĥ with eigenvalue ε, and let the wavefunctions ψ1′ , . . . , ψK′ be a distinct orthonormal basis of the same
eigenspace. Since both sets of functions serve as a basis for the samefunctional space, there exists an irreducible unitary representation(norm-preserving map) ℛ : Gε → U(K) that maps any element in
the group Gε composed of all wavefunctions that are eigenfunctionsof Ĥ, with eigenvalue ε, to a K × K unitary matrix that operates ona K-dimensional Hilbert space V(K), due to Ĥ being self-adjoint,
such that⎡ ψ ⎤ ⎡ ψ′ ⎤⎢ 1⎥ ⎢ 1⎥⎢ ⎥ ⎢ ⎥⎥ ⎢ ⎥ℛ⎢
⎢ ⋮ ⎥ = ⎢ ⋮ ⎥,⎢ ⎥ ⎢ ′⎥
⎢ψK ⎥ ⎢ψK ⎥⎣ ⎦ ⎣ ⎦(A1)

AUTHOR DECLARATIONSConflict of Interest
The authors have no conflicts to disclose.
Author Contributions
José Romero: Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). Paulo Limão-Vieira:Formal analysis (equal); Supervision (equal); Validation (equal);
Writing – review & editing (equal). Thana Maihom: Formal analysis (equal); Methodology (equal); Validation (equal); Writing –review & editing (equal). Kersti Hermansson: Formal analysis(equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Michael

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

in which taking the norm of both sides of the equation concludes theproof.
◻
Lemma 2. An eigenspace of a Hamiltonian operator cannot bea smaller subset of an eigenspace of another distinct (more than just a
constant) Hamiltonian.
Proof. Let ψ1 , . . . , ψK ∈ L2 (R3N ) be an orthonormal basis of
the eigenspace of a N electron Hamiltonian Ĥ = T̂ + Û EE + Û Ext of
eigenvalue ε and ψ1 , . . . , ψK ′ ∈ L2 (R3N ) an orthonormal basis of theeigenspace of a distinct Hamiltonian Ĥ ′ = T̂ + Û EE + Û ′Ext , in whichK ′ ≤ K. It follows from integration over any Ω ⊆ R3N that it must betrue that

160, 235101-12

09 August 2024 10:12:59

J.R. acknowledges the support from the Portuguese national
funding agency Fundação para a Ciência e a Tecnologia (FCT) for
the PhD scholarship Grant No. PD/BD/142846/2018 and together
with P.L.V. the research Grant No. CEFITEC UIDB/00068/2020(DOI: 10.54499/UIDB/00068/2020). This work was also supportedby the Radiation Biology and Biophysics Doctoral Training Program (RaBBiT, Grant No. PD/00193/2012; UCIBIO Grant No.
UIDB/04378/2020). K.H. acknowledges the support from theSwedish Research Council (Vetenskapsrådet), the Swedish National
Infrastructure for Computing (SNIC/NAISS), and the Swedish
national strategic e-Science program eSSENCE. The results presented here have been achieved, in part, using the LEO High
Performance Computing (HPC) infrastructure of the Universityof Innsbruck. This work was partially carried out within theframework of the EUROfusion Consortium and received funding
from the Euratom Research and Training programme under Grant
Agreement No. 101052200-EUROfusion. The views and opinions
expressed herein do not necessarily reflect those of the European
Commission.

ARTICLEThe Journalof Chemical Physics

ARTICLE

⎧
K′
K′
⎪
⎪
⎪
⎪
ε
⟨ψ
⟩
=
∑ n Ω ∑ ⟨ψn ∣T̂ + Û EE + Û Ext ∣ψn ⟩Ω ,
⎪
⎪
⎪ n=1
n=1
⎨ K′
K′
⎪
⎪
′
⎪
⎪
ε
⟨ψ
⟩
=
⟨ψn ∣T̂ + Û EE + Û ′Ext ∣ψn ⟩Ω ,
∑
∑
n
⎪
Ω
⎪
⎪
n=1
⎩ n=1
K′

⇒

(A2a)

pubs.aip.org/aip/jcp

TABLE II. Comparison of the Vertical Electron Affinity (VEA) values of several few
small molecules obtained from high accuracy ab initio calculations, CCSD(T)/cc-pVTZ
level of theory, and our order ten full electron and ECP variants of our model. All thevalues reported are expressed in eV.

CCSD(T)

ECP variant

Full elec. variant−4.217−3.444−1.433−4.489−3.6762.966−52.189
43.9803.878 × 101
4.772 × 103
6.828 × 104
3.304 × 101

K′′(ε − ε ) ∑ ⟨ψn ⟩Ω = ∑ ⟨ψn ∣Û Ext − Û ′Ext ∣ψn ⟩Ω ,
n=1
n=1
(A2b)

and therefore, Û Ext − Û ′Ext must be simply a constant, implying that
′

H2
N2
O2
CO2

K = K. Furthermore, if we assume the vacuum level to be the zerothenergy level, then Û Ext = Û ′Ext .
◻
Theorem 1. There is a unique functional of the ground state
eigendensity of a Hamiltonian that maps it to its external potential.
Proof. Let Ĥ and Ĥ ′ be two distinct N electron Hamiltonian operators whose respective lowest eigenvalues are ε0 andε′0 . Furthermore, let the degeneracies of their respective groundstate eigenspaces be K and K ′ ; then, their respective ground state
eigendensities are∣ψk (r, r2 , . . . , rN )∣2
dr2 ⋅ ⋅ ⋅ drN ,
K/N
R k=1
K

ρε0 (r) = ∫ 3 ⋅ ⋅ ⋅ ∫ 3 ∑
R

2

K′

(A3b)

and thus, by hypothesis,
ρε0 (r) = ρε′0 (r) = ρ0 (r)∀r ∈ R3.

(A4)

Without loss of generality, we get from Lemma 2 and by integrating
that
⎧
K K′
⎪
⟨ψn ∣T̂ + Û EE + Û Ext ∣ψn ⟩
⎪
⎪
ε
=
,
∑
∑
⎪
0
⎪
⎪
KK ′
⎪
n=1 n′ =1⎨′′
K K ′ ⟨ψ ′ ′ ∣T̂ + Û⎪EE + Û Ext ∣ψn′ ⟩
⎪
′
n
⎪
⎪
,
ε
=
∑
∑
⎪
0
′
⎪
⎪
KK
n=1 n′ =1
⎩
⎧
K K′
⎪
⟨ψn ∣T̂ + Û EE + Û ′Ext ∣ψn ⟩
′
⎪
⎪
,
ε
<
∑ ∑
⎪
0
⎪
⎪
KK ′
⎪
n=1 n′ =1⇒ ⎨′K K ′ ⟨ψ ′ ′ ∣T̂ + Û⎪EE + Û Ext ∣ψn′ ⟩
⎪
n
⎪
⎪
ε
<
,
∑
∑
0
⎪
⎪
⎪
KK ′
n=1 n′ =1
⎩

(A5a)

(A5b)

⎧
K K′
⎪
⟨ψn ∣Û ′Ext − Û Ext ∣ψn ⟩⎪⎪,ε′0 < ε0 + ∑ ∑
⎪
⎪
⎪
KK ′
⎪
n=1 n′ =1
⇒ ⎨
′
′
′
′
K K
⎪
⟨ψn′ ∣Û Ext − Û Ext ∣ψn′ ⟩⎪⎪⎪ε0 < ε′0 + ∑ ∑
,
⎪
⎪
⎪
KK ′
n=1 n′ =1
⎩

(A5c)

⎧
K′
⎪
ρ0 (r)(u′ (r) − u(r))
⎪
′
⎪
⎪
ε
<
ε
dr,
0+ ∑ ∫ 3
0
⎪
⎪
NK ′
⎪
R
n′ =1
⇒ ⎨
′
K
⎪
ρ0 (r)(u(r) − u (r))⎪⎪⎪ε0 < ε′0 + ∑ ∫ 3
dr,
⎪
⎪
⎪
NK
n=1 R
⎩

(A5d)

which is a contradiction identical to the one found in Hohenberg
and Kohn’s original proof.32◻

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

As one would expect, neither the full electron nor ECP orderten variants of our model provide accurate estimations of the vertical electron affinities of the various small molecules as summarizedin Table II. Specifically, the full electron variant shows a much lesser
accuracy, in terms of magnitude, of the two when compared to theequivalent ab initio calculated result with discrepancies with two or
more orders of magnitude higher. These results come to no surprisegiven that our model is unable to model chemical bonds (refer toFig. S6 of the supplementary material). Nevertheless, a detailed
exploration and reporting of the strengths and weaknesses ofour model is essential for those aiming to conduct accurate MD
simulations.

REFERENCES
1

W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein,“Comparison of simple potential functions for simulating liquid water,” J. Chem.
Phys. 79, 926–935 (1983).
2
M. A. González and J. L. F. Abascal, “A flexible model for water based on
TIP4P/2005,” J. Chem. Phys. 135, 224516 (2011).
3
M. M. Ghahremanpour, J. Tirado-Rives, and W. L. Jorgensen, “Refinement of
the optimized potentials for liquid simulations force field for thermodynamics anddynamics of liquid alkanes,” J. Phys. Chem. B 126, 5896–5907 (2022).
4
M. W. Mahoney and W. L. Jorgensen, “A five-site model for liquid water and the
reproduction of the density anomaly by rigid, nonpolarizable potential functions,”
J. Chem. Phys. 112, 8910–8922 (2000).
5
W. L. Jorgensen and J. Tirado-Rives, “Potential energy functions for atomic-level
simulations of water and organic and biomolecular systems,” Proc. Natl. Acad. Sci.
U. S. A. 102, 6665–6670 (2005).
6
L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives, and W. L. Jorgensen, “1.14∗ CM1ALBCC: Localized bond-charge corrected CM1A charges for condensed-phasesimulations,” J. Phys. Chem. B 121, 3864–3870 (2017).
7
L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives, and W. L. Jorgensen, “LigParGen
web server: An automatic OPLS-AA parameter generator for organic ligands,”
Nucleic Acids Res. 45, W331–W336 (2017).
8
D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz,
A. Onufriev, C. Simmerling, B. Wang, and R. J. Woods, “The Amber biomolecularsimulation programs,” J. Comput. Chem. 26, 1668–1688 (2005).
9D. Case, H. Aktulga, K. Belfon, I. Ben-Shalom, J. Berryman, S. Brozell,
D. Cerutti, T. Cheatham III, G. C. Andrés Cisneros, V. Cruzeiro, T. Darden,
R. Duke, G. Giambasu, M. Gilson, H. Gohlke, A. Goetz, R. Harris, S. Izadi, S.
Izmailov, K. Kasavajhala, M. Kaymak, E. King, A. Kovalenko, T. Kurtzman, T. Lee,

160, 235101-13

09 August 2024 10:12:59∣ψk′ (r, r2 , . . . , rN )∣dr2 ⋅ ⋅ ⋅ drN ,K ′ /N
R k=1

ρε′0 (r) = ∫ 3 ⋅ ⋅ ⋅ ∫ 3 ∑
R

(A3a)

APPENDIX B: COMPARING VEA
OF SMALL MOLECULESThe Journalof Chemical Physics

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

pubs.aip.org/aip/jcp

29

S. Naseem-Khan, J.-P. Piquemal, and G. A. Cisneros, “Improvement of theGaussian electrostatic model by separate fitting of Coulomb and exchangerepulsion densities and implementation of a new dispersion term,” J. Chem. Phys.
155, 194103 (2021).
30
S. Vandenbrande, M. Waroquier, V. V. Speybroeck, and T. Verstraelen, “Themonomer electron density force field (MEDFF): A physically inspired model fornoncovalent interactions,” J. Chem. Theory Comput. 13, 161–179 (2017).
31
J. Romero, P. Limão-Vieira, K. Hermansson, and M. Probst, “A simple electrondensity based force field model for high-energy interactions between atoms andmolecules,” J. Phys. Chem. A 128, 1163–1172 (2024).
32
P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136,
B864–B871 (1964).
33
M. Levy and J. P. Perdew, “The constrained search formulation of density functional theory,” in Density Functional Methods in Physics, edited by R. M. Dreizler
and J. da Providência (Springer US, Boston, MA, 1985), pp. 11–30.
34
J. Katriel, F. Zahariev, and K. Burke, “Symmetry and degeneracy in density
functional theory,” Int. J. Quantum Chem. 85, 432–435 (2001).
35
Y. Jung, A. Sodt, P. M. W. Gill, and M. Head-Gordon, “Auxiliary basis expansions for large-scale electronic structure calculations,” Proc. Natl. Acad. Sci. U. S.
A. 102, 6692–6697 (2005).
36
G. A. Cisneros, “Application of Gaussian electrostatic model (GEM) distributed
multipoles in the AMOEBA force field,” J. Chem. Theory Comput. 8, 5072–5080(2012).
37
A. D. Becke, “Density-functional thermochemistry. III. The role of exact
exchange,” J. Chem. Phys. 98, 5648–5652 (1993).
38
S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate
ab initio parametrization of density functional dispersion correction (DFT-D) for
the 94 elements H-Pu,” J. Chem. Phys. 132, 154104 (2010).
39
W. J. Stevens, H. Basch, and M. Krauss, “Compact effective potentials and efficient shared-exponent basis sets for the first- and second-row atoms,” J. Chem.
Phys. 81, 6026–6033 (1984).
40
W. J. Stevens, M. Krauss, H. Basch, and P. G. Jasien, “Relativistic compact effective potentials and efficient, shared-exponent basis sets for the third-fourth-and
fifth-row atoms,” Can. J. Chem. 70, 612–630 (1992).
41
T. R. Cundari and W. J. Stevens, “Effective core potential methods for thelanthanides,” J. Chem. Phys. 98, 5555–5565 (1993).
42
T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations.
I. The atoms boron through neon and hydrogen,” J. Chem. Phys. 90, 1007–1023(1989).
43
D. E. Woon and T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon,” J. Chem. Phys. 98,
1358–1371 (1993).
44
D. E. Woon and T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. V. Core-valence basis sets for boron through neon,” J. Chem.
Phys. 103, 4572–4585 (1995).
45
A. K. Wilson, T. van Mourik, and T. H. Dunning, “Gaussian basis sets for use
in correlated molecular calculations. VI. Sextuple zeta correlation consistent basis
sets for boron through neon,” J. Mol. Struct.: THEOCHEM 388, 339–349 (1996).
46
J. Cížek, in Advances in Chemical Physics, edited by R. LeFebvre and C. Moser(John Wiley & Sons, Inc., 1969).
47
G. D. Purvis and R. J. Bartlett, “A full coupled-cluster singles and doubles model:
The inclusion of disconnected triples,” J. Chem. Phys. 76, 1910–1918 (1982).
48
J. A. Pople, M. Head-Gordon, and K. Raghavachari, “Quadratic configurationinteraction. A general technique for determining electron correlation energies,”
J. Chem. Phys. 87, 5968–5975 (1987).
49
G. E. Scuseria, C. L. Janssen, and H. F. Schaefer, “An efficient reformulation of
the closed-shell coupled cluster single and double excitation (CCSD) equations,”
J. Chem. Phys. 89, 7382–7387 (1988).
50
G. E. Scuseria and H. F. Schaefer, “Is coupled cluster singles and doubles(CCSD) more computationally intensive than quadratic configuration interaction(QCISD)?,” J. Chem. Phys. 90, 3700–3703 (1989).
51
R. A. Kendall, T. H. Dunning, and R. J. Harrison, “Electron affinities of the firstrow atoms revisited. Systematic basis sets and wave functions,” J. Chem. Phys. 96,
6796–6806 (1992).

160, 235101-14

09 August 2024 10:12:59

S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, M. Machado, V. Man, M. Manathunga, K. Merz, Y. Miao, O. Mikhailovskii, G. Monard, H. Nguyen, K. O’Hearn,
A. Onufriev, F. Pan, S. Pantano, R. Qi, A. Rahnamoun, D. Roe, A. Roitberg,
C. Sagui, S. Schott-Verdugo, A. Shajan, J. Shen, C. Simmerling, N. Skrynnikov,
J. Smith, J. Swails, R. Walker, J. Wang, J. Wang, H. Wei, R. Wolf, X. Wu, Y. Xiong,
Y. Xue, D. York, S. Zhao, and P. Kollman, AMBER 2022, 2022.
10
T. A. Halgren, “The representation of van der Waals (vdW) interactions inmolecular mechanics force fields: potential form, combination rules, and vdWparameters,” J. Am. Chem. Soc. 114, 7827–7843 (1992).
11
T. A. Halgren, “Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94,” J. Comput. Chem. 17, 490–519(1996).
12
T. A. Halgren, “Merck molecular force field. II. MMFF94 van der Waals andelectrostatic parameters for intermolecular interactions,” J. Comput. Chem. 17,
520–552 (1996).
13
T. A. Halgren, “Merck molecular force field. III. Molecular geometries andvibrational frequencies for MMFF94,” J. Comput. Chem. 17, 553–586 (1996).
14
T. A. Halgren and R. B. Nachbar, “Merck molecular force field. IV. Conformational energies and geometries for MMFF94,” J. Comput. Chem. 17, 587–615(1996).
15
T. A. Halgren, “Merck molecular force field. V. Extension of MMFF94
using experimental data, additional computational data, and empirical rules,”J. Comput. Chem. 17, 616–641 (1996).
16
F. Maseras and K. Morokuma, “IMOMM: A new integrated ab initio + molecular mechanics geometry optimization scheme of equilibrium structures andtransition states,” J. Comput. Chem. 16, 1170–1179 (1995).
17
S. Dapprich, I. Komáromi, K. Byun, K. Morokuma, and M. J. Frisch,“A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives,” J. Mol.
Struct.: THEOCHEM 461–462, 1–21 (1999).
18
T. Vreven and K. Morokuma, “On the application of the IMOMO (integrated
molecular orbital + molecular orbital) method,” J. Comput. Chem. 21, 1419–1432(2000).
19
M. Svensson, S. Humbel, and K. Morokuma, “Energetics using the single point
IMOMO (integrated molecular orbital + molecular orbital) calculations: Choicesof computational levels and model system,” J. Chem. Phys. 105, 3654–3661 (1996).
20
M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, and
K. Morokuma, “ONIOM: A multilayered integrated MO + MM method forgeometry optimizations and single point energy predictions. A test for DielsAlder reactions and Pt(P(t-Bu)3 )2 + H2 oxidative addition,” J. Phys. Chem. 100,
19357–19363 (1996).
21
T. Matsubara, S. Sieber, and K. Morokuma, “A test of the new ‘integrated
MO + MM’ (IMOMM) method for the conformational energy of ethane andn-butane,” Int. J. Quantum Chem. 60, 1101–1109 (1996).
22
S. Humbel, S. Sieber, and K. Morokuma, “The IMOMO method: Integrationof different levels of molecular orbital approximations for geometry optimizationof large systems: Test for n-butane conformation and SN 2 reaction: RCl+Cl− ,”J. Chem. Phys. 105, 1959–1967 (1996).
23
D. G. Pettifor and I. I. Oleinik, “Analytic bond-order potentials beyond TersoffBrenner. I. Theory,” Phys. Rev. B 59, 8487–8499 (1999).
24
D. G. Pettifor, M. Finnis, D. Nguyen-Manh, D. Murdick, X. Zhou, andH. Wadley, “Analytic bond-order potentials for multicomponent systems,” Mater.
Sci. Eng.: A 365, 2–13 (2004).
25
J. Behler and M. Parrinello, “Generalized neural-network representationof high-dimensional potential-energy surfaces,” Phys. Rev. Lett. 98, 146401(2007).
26
N. Bernstein, G. Csányi, and V. L. Deringer, “De novo exploration and selfguided learning of potential-energy surfaces,” npj Comput. Mater. 5, 99 (2019).
27
G. A. Cisneros, K. T. Wikfeldt, L. Ojamäe, J. Lu, Y. Xu, H. Torabifard, A. P.
Bartók, G. Csányi, V. Molinero, and F. Paesani, “Modeling molecular interactions
in water: From pairwise to many-body potential energy functions,” Chem. Rev.
116, 7501–7528 (2016).
28
H. Gökcan, E. Kratz, T. A. Darden, J.-P. Piquemal, and G. A. Cisneros,“QM/MM simulations with the Gaussian electrostatic model: A density-basedpolarizable potential,” J. Phys. Chem. Lett. 9, 3062–3067 (2018).

ARTICLEThe Journalof Chemical Physics

52

J. Chem. Phys. 160, 235101 (2024); doi: 10.1063/5.0210949© Author(s) 2024

pubs.aip.org/aip/jcp

62

J. A. Nelder and R. Mead, “A simplex method for function minimization,”
Comput. J. 7, 308–313 (1965).
63
J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence properties of the Nelder–Mead simplex method in low dimensions,” SIAM J. Optim. 9,
112–147 (1998).
64
F. Gao and L. Han, “Implementing the Nelder-Mead simplex algorithm withadaptive parameters,” Comput. Optim. Appl. 51, 259–277 (2012).
65
C. G. Broyden, “The convergence of a class of double-rank minimization algorithms 1. General considerations,” IMA J. Appl. Math. 6, 76–90(1970).
66
R. Fletcher, “A new approach to variable metric algorithms,” Comput. J. 13,
317–322 (1970).
67
D. Goldfarb, “A family of variable-metric methods derived by variational
means,” Math. Comput. 24, 23–26 (1970).
68
D. F. Shanno, “Conditioning of quasi-Newton methods for function
minimization,” Mathematics of Computation 24, 647–656 (1970).
69
J. Nocedal, “Updating quasi-Newton matrices with limited storage,” Math.
Comput. 35, 773–782 (1980).
70
D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale
optimization,” Math. Program. 45, 503–528 (1989).
71
J. Nocedal and J. W. Stephen, in Numerical Optimization, 2nd ed. (Springer,
New York, 2006).
72
J. Revels, M. Lubin, and T. Papamarkou, “Forward-mode automatic differentiation in Julia,” arXiv:1607.07892 [cs.MS] (2016).
73
Z.-Z. Yang and C.-S. Wang, “Atom–bond electronegativity equalization
method. 1. Calculation of the charge distribution in large molecules,” J. Phys.
Chem. A 101, 6315–6321 (1997).
74
P. Bultinck, W. Langenaeker, P. Lahorte, F. De Proft, P. Geerlings,
M. Waroquier, and J. P. Tollenaere, “The electronegativity equalization method
I: Parametrization and validation for atomic charge calculations,” J. Phys. Chem.
A 106, 7887–7894 (2002).
75
Maxima, Maxima, a computer algebra system. version 5.47.0, 2023.
76
J. Romero, MolFFSim - A C++ implementation of an OFDFT based molecularforce field model, 2024.
77
G. Guennebaud, B. Jacob et al., Eigen v3, 2010.
78
N. Okazaki, libLBFGS: L-BFGS library written in C, 2014.
79
A. M. M. Leal, autodiff, a modern, fast and expressive C++ library for automaticdifferentiation, 2018.

160, 235101-15

09 August 2024 10:12:59

J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case,“Development and testing of a general amber force field,” J. Comput. Chem. 25,
1157–1174 (2004).
53
N. M. O’Boyle, M. Banck, C. A. James, C. Morley, T. Vandermeersch, and G. R.
Hutchison, “Open Babel: An open chemical toolbox,” J. Cheminf. 3, 33 (2011).
54
M. S. Gordon, M. A. Freitag, P. Bandyopadhyay, J. H. Jensen, V. Kairys, and W.
J. Stevens, “The effective fragment potential method: A QM-based mm approach
to modeling environmental effects in chemistry,” J. Phys. Chem. A 105, 293–307(2001).
55
M. S. Gordon, L. Slipchenko, H. Li, and J. H. Jensen, “Chapter 10. The effective
fragment potential: A general method for predicting intermolecular interactions,”
Annu. Rep. Comput. Chem. 3, 177–193 (2007).
56
L. V. Slipchenko and M. S. Gordon, “Damping functions in the effective
fragment potential method,” Mol. Phys. 107, 999–1016 (2009).
57
I. A. Kaliman and L. V. Slipchenko, “LIBEFP: A new parallel implementation of
the effective fragment potential method as a portable software library,” J. Comput.
Chem. 34, 2284–2292 (2013).
58
M. Riera, E. P. Yeh, and F. Paesani, “Data-driven many-body models for molecular fluids: CO2 /H2 O mixtures as a case study,” J. Chem. Theory Comput. 16,
2246–2257 (2020).
59
S. Yue, M. Riera, R. Ghosh, A. Z. Panagiotopoulos, and F. Paesani,“Transferability of data-driven, many-body models for CO2 simulations in thevapor and liquid phases,” J. Chem. Phys. 156, 104503 (2022).
60
M. Riera, C. Knight, E. F. Bull-Vulpe, X. Zhu, H. Agnew, D. G. A. Smith, A. C.
Simmonett, and F. Paesani, “MBX: A many-body energy and force calculator fordata-driven many-body simulations,” J. Chem. Phys. 159, 054802 (2023).
61
M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R.
Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H.
P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young,
F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D.
Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada,
M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda,
O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F.
Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T.
A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant,
S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi,
J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J.
Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.

ARTICLE