††thanks: Contributed equally††thanks: Contributed equally

# Non-covalent interactions between molecular dimers (S66) in electric fields

Max Schwilk  Faculty of Physics, University of Vienna, Kolingasse 14-16, 1090
Vienna, Austria  Institute of Physical Chemistry and National Center for
Computational Design and Discovery of Novel Materials (MARVEL), Department of
Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel,
Switzerland  Pál D. Mezei  Institute of Physical Chemistry and National Center
for Computational Design and Discovery of Novel Materials (MARVEL), Department
of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel,
Switzerland  Diana N. Tahchieva  Institute of Physical Chemistry and National
Center for Computational Design and Discovery of Novel Materials (MARVEL),
Department of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056
Basel, Switzerland  O. Anatole von Lilienfeld
<EMAIL_ADDRESS>Faculty of Physics, University of Vienna,
Kolingasse 14-16, 1090 Vienna, Austria  Institute of Physical Chemistry and
National Center for Computational Design and Discovery of Novel Materials
(MARVEL), Department of Chemistry, University of Basel, Klingelbergstrasse 80,
CH-4056 Basel, Switzerland

###### Abstract

Fine tuning and microscopic control of van der Waals interactions through
oriented external electric fields (OEEF) mandates an accurate and systematic
understanding of intermolecular response properties. Having taken exploratory
steps into this direction, we present a systematic study of interaction
induced dipole electric properties of all molecular dimers in the S66 set,
relying on CCSD(T)-F12b/aug-cc-pVDZ-F12 as reference level of theory. For
field strengths up to $\approx$5 GV m-1 the interaction induced electric
response beyond second order is found to be insignificant. Large interaction
dipole moments (i.e. dipole moment changes due to van der Waals binding) are
observed in the case of hydrogen bonding oriented along the intermolecular
axis, and mostly small interaction dipole moments are found in dimers bonded
by $\pi$-stacking or London dispersion. The interaction polarizabilities (i.e.
polarizability changes due to van der Waals binding) were generally found to
be small but always with a positive-valued principal component approximately
aligned with the intermolecular axis, and two other negative-valued
components. Energy decompositions according to symmetry adapted perturbation
theory (SAPT0/jun-cc-pVDZ) suggest that electrostatics dominates the
interaction dipole moment, with exchange and induction contributing on a
smaller scale, and with dispersion having the smallest effect. First-order
SAPT0 decomposition into monomer-resolved contributions enables us to
establish a quantitative link between electric properties of monomers and
dimers, which is found to be in qualitative agreement with the coupled cluster
reference method. Using the aug-cc-pVQZ basis and non-empirical PBE semilocal
exchange-correlation kernels, we also assess how density functional
approximations in the nonlocal exchange and correlation parts affect the
predictive accuracy: While dRPA@PBE0 based predictions are in excellent
overall agreement with coupled cluster results, the computationally more
affordable LC-$\omega$PBE0-D3 level of theory also yields reliable results
with relative errors below 5%. PBE alone, even when dispersion corrected,
produces larger errors in interaction dipole moments ($\approx$ 10%) and
polarizabilities ($\approx$ 20%). We also resolve the mutual impact of the
three dimensions of the OEEF, and we present a discussion of the
intermolecular distance dependence of the perturbations.

††preprint: AIP/123-QED

##  I Introduction

Modulating chemical interactions of molecules with an external electric field
is a rather intuitive concept that nevertheless has not yet found its way to
widespread synthetic applications. An external electric field can be present
as a result of a crystal lattice or a solvent environment (with field
strengths of up to 8 GV m-1),spackman2007use ; fried2017electric  and it also
occurs between the tip and the sample in scanning tunneling
microscopy.hong1998electron  It was only in the 1970s that the capability of
external electric fields to catalyze chemical reactions was demonstrated by
exploring the “intramolecular salt effect”.pocker1970electro ; SMITH19912617
The experimental work also triggered computational efforts on probing the
chemical modulation capabilities of the oriented external electric field
(OEEF), such as accelerating proton transfer reactions,arabi2011effects
catalyzing the Diels-Alder reaction,meir2010oriented  or controlling reaction
selectivity.shaik2004external  Experimental progress included the activation
of enzymesthuren1987triggering  and the reversible isomerization by the OEEF
at the tip of a scanning tunneling microscope.alemani2006electric  In 2016,
Aragonès _et al._aragones2016  succeeded in using the OEEF at the tip of the
scanning tunneling microscope to catalyze the Diels-Alders reaction. This
breakthrough led to a recent surge in efforts to unveil the full potential of
OEEFs as a chemical tool. Impressive progress has been made due to synergies
of experimental and theoretical efforts.shaik2016oriented ; Nejad2014 ;
shaik2020 ; kaestner2020

Strong electric fields can also occur as an intramolecular effect. An internal
electric field can arise in the cavity of frustrated Lewis
pairsgrimme2010mechanism . An internal electric field of mainly dipolar
character is generated at an active site of enzymes if two oppositely charged
residues are in its close vicinity.lehle2005probing ; suydam2006electric  In
proteins, induced internal electric fields can be measured by the vibrational
Stark effect.fried2015measuring  An internal electric field can also be
adjusted by point mutations in proteins (such as replacing residues in the
amino acid sequence) for example to accelerate photosynthetic
reactions,gopher1985effect ; popovic1986electric ; franzen1990electric  or to
alter the enzyme activity.moroney1984effect

Noncovalent interactions play an important role for understanding supra-
molecular processes, developing new synthesis pathways,knowles2010attractive
designing molecular recognition processes,metrangolo2008halogen , or
engineering materials.zhang2003fabrication  For example, within a recent
study, molecular dynamics simulations revealed that a novel paracetamol
polymorph can be created in the presence of a strong electric field (1.5 GV
m-1).parks2017molecular  Or, Kleshchonok and Tkatchenko studied intermolecular
interactions in the presence of point charges in biologically relevant
noncovalent dimers,kleshchonok2018tailoring  obtaining field induced
dispersion effects amounting to up to 35% of the total interaction energy.
Even more surprisingly, they found that positive and negative point charges in
the vicinity of noncovalent dimers may weaken or strengthen the dispersion
interaction. Similarly, non-trivial static effects can also emerge when it
comes to the screening of dispersion in layered systems ambrosetti2018hidden .

Kozłowska _et al._kozlowska  computed interaction induced dipole moments of
endohedral small polar molecules in carbon nanotubes of varying diameter with
the PNO-LCCSD-F12Schwilk2017JCTC ; Ma2017JCTC_1  local coupled cluster method
and compared various density functional approximations (DFAs), addressing the
impact of exchange-correlation kernels, orbital relaxation, and dispersion
treatments. Their results underline the sensitivity of DFAs with respect to
the dispersion treatment for electric properties. In a follow-up study, DFAs
were benchmarked for interaction-induced electric properties of small
hydrogen-bonded dimers confined by a simple quadratic potential.kozlowska2  It
should be noted that, in general, interaction induced electric properties may
also contain important contributions from nuclear relaxation and molecular
vibration.dft_interaction_induced3

In this work, we computationally quantify and examine the effect of a dipolar
OEEF on noncovalent interaction energies for the S66 set with 66 noncovalent
dimers of small and medium size molecules.s66 ; s66x8  The focus of our study
lies on probing the capability of a dipolar OEEF to modulate the
intermolecular interaction, on probing computationally more efficient
methodologies to describe the effect, as well as on obtaining a qualitative
understanding of the most relevant physical interactions.

##  II Methodology and computational details

The equilibrium S66 noncovalent dimer geometries are taken from the
literature,s66x8  and no further geometry optimization is applied. For all S66
dimers, we calculate the field $\bm{F}$ dependent total potential energies of
interaction $\Delta E$ according to eq 1.

$\Delta E(\bm{F})=E^{AB}(\bm{F})-E^{A}(\bm{F})-E^{B}(\bm{F})$ (1)

In order to describe the change in the interaction energy induced by the
applied OEEFs, we compute the interaction dipole moments ($\Delta\bm{\mu}$)
and interaction polarizabilities ($\Delta\bm{\alpha}$) components
$i,j\in\\{x,y,z\\}$ according to eqs 2 and 3 by finite differentiation. The
accuracy of this second order expansion is verified.

$\Delta\mu_{i}=-\frac{\partial\Delta E(\bm{F})}{\partial
F_{i}}\Bigr{|}_{\bm{F}=\bm{0}}$ (2)

$\Delta\alpha_{ij}=-\frac{\partial^{2}\Delta E(\bm{F})}{\partial F_{i}\partial
F_{j}}\Bigr{|}_{\bm{F}=\bm{0}}$ (3)

An OEEF is applied along the intermolecular (z) axis (going through the two
monomer centers of mass and directed such that the z-component of the sum of
the two monomer dipole moments is positive) with values of 0.000 a.u.,
$\pm$0.001 a.u., $\pm$0.005 a.u., and $\pm$0.01 a.u. (0.01 a.u. $\approx$ 5.14
GV m-1), covering a range of frequently encountered field strengths. We also
apply a three-dimensional grid of OEEF strengths $\\{F_{x},F_{y},F_{z}\\}$
with values of 0.00 a.u. and $\pm$0.01 a.u. in each dimension.

Furthermore, we analyze orientation and strength of the interaction dipole
moment $\Delta\bm{\mu}$ as well as the interaction polarizability principal
components $\bm{\nu}^{r}$, $r\in\\{a,b,c\\}$ as given in eq 4.

$\displaystyle\left[\begin{matrix}v^{a}_{x}&v^{b}_{x}&v^{c}_{x}\\\
v^{a}_{y}&v^{b}_{y}&v^{c}_{y}\\\
v^{a}_{z}&v^{b}_{z}&v^{c}_{z}\end{matrix}\right]^{T}\left[\begin{matrix}\Delta\alpha_{xx}&\Delta\alpha_{xy}&\Delta\alpha_{xz}\\\
\Delta\alpha_{yx}&\Delta\alpha_{yy}&\Delta\alpha_{yz}\\\
\Delta\alpha_{zx}&\Delta\alpha_{zy}&\Delta\alpha_{zz}\end{matrix}\right]\left[\begin{matrix}v^{a}_{x}&v^{b}_{x}&v^{c}_{x}\\\
v^{a}_{y}&v^{b}_{y}&v^{c}_{y}\\\
v^{a}_{z}&v^{b}_{z}&v^{c}_{z}\end{matrix}\right]$ (4)
$\displaystyle=\left[\begin{matrix}\Delta\alpha_{aa}&0&0\\\
0&\Delta\alpha_{bb}&0\\\ 0&0&\Delta\alpha_{cc}\end{matrix}\right]$

For the sake of completeness, we also present two simpler but more frequently
used quantities in the supporting information (SI): the isotropic and
anisotropic components of the interaction polarizability in the $(x,y,z)$
frame. It has been shown that fourth order perturbation theory effects are
important in the accurate description of field induced electric
properties.dft_interaction_induced  We therefore choose the coupled cluster
singles, doubles, and perturbative triples [CCSD(T)]mhg1989  method for
reference. The coupled cluster reference calculations are performed with the
Molpro quantum chemistry softwareMOLPRO_brief ). Since the basis set
convergence of CCSD(T) is fairly slow and comes with high computational
overhead even for relatively small systems, we use the explicitly correlated
F12b methodknizia:2009  [3*A(fix,nox) variant as implemented in Molpro] with
the aug-cc-pVDZ-F12Dunning1989 ; peterson2008systematically ; nitai2017  basis
set and the aug-cc-pVTZ-F12 auxiliary basis set.nitai2017 ; Yousaf2008 ;
Yousaf2009  This level of theory has shown to yield a reasonably small
interaction energy error (0.048 kcal mol-1 on average) at zero field strength
for the benchmark S66 energies.s66x8

We use the coupled cluster reference for a systematic assessment of different
density functional theory (DFT) approximations based on the semilocal
nonempirical PBE functional.pbe  Due to the rigorous roots in fundamental
physics of the PBE functional we expect it to cover short range exchange and
correlation at least qualitatively correct in the context of electric response
properties.dft_dipole_moment ; dft_polarizability  Approximations tested
include the PBE-D3,pbe ; dftd3 ; dftd3bj , PBE-VV10,pbe ; vv10  PBE0-D3,pbe0 ;
dftd3 ; dftd3bj , PBE0-VV10,pbe0 ; vv10  HSE06-D3,hse03 ; hse06 ; dftd3 ;
dftd3bj  LC-$\omega$PBE-D3,lcwpbe ; dftd3 ; dftd3bj  PBE0-2,pbe02 ,
dRPA75,drpa75 , dRPA@PBE, and dRPA@PBE0drpa_drccd ; drpa_equivalence ;
drpa_implementation  functionals. Hence, the numerical results can be expected
to be fairly transparent regarding effects due to dispersion corrections,
exact exchange mixing, and nonlocal correlation components, and encode the
prospect that conclusions are mostly transferable to other density functional
families. All DFT computations are performed with the aug-cc-pVQZ basis
setDunning1989  in the MRCC quantum chemistry program suite.MRCC

To gain further insights, we also employ symmetry-adapted perturbation theory
(SAPT)jeziorski1994perturbation  to decompose the interaction energy, dipole
moment, and polarizability into physically meaningful contributions (such as
electrostatic, exchange, induction, and dispersion). In SAPT, the perturbed
Hamiltonian as given in eq 5 depends on two coupling constants $\zeta$ and
$\lambda$, where $\hat{V}$ is the intermonomer interaction operator, $\hat{F}$
is the unperturbed operator defined as the sum of monomer Fock operators, and
$\hat{W}$ is the intramonomer correlation operator defined as the difference
between the sum of monomer Hamiltonians and the sum of monomer Fock operators.

$\displaystyle\hat{H}(\zeta,\lambda)=\hat{F}+\zeta\hat{V}+\lambda\hat{W}$ (5)

The interaction energy is obtained by a twofold perturbation expansion given
in eq 6 where v and w define the orders of perturbation in $\zeta$ and
$\lambda$. Hence, at zeroth order, the SAPT solution consists of two non-
interacting monomers, each of which is described by a separate Hartree-Fock
(HF) wave function.

$\displaystyle\Delta
E^{\mathrm{SAPT}}=\sum_{v=1}^{\infty}\sum_{w=0}^{\infty}{\Delta E^{(vw)}}$ (6)

In our interaction energy decomposition, we apply the SAPT0sapt0_ref1 ;
sapt0_ref2  method that contains only zeroth-order terms in $\lambda$ as given
in eq 7. The jun-cc-pVDZ basis setseasonal  was used and all SAPT computations
were performed with the Psi4 quantum chemistry software.psi4  The accuracy of
SAPT0 as a qualitative computational method is also verified.

$\displaystyle\Delta E^{\mathrm{SAPT0}}=\Delta
E^{\mathrm{SAPT0}}_{\mathrm{elst}}+\Delta
E^{\mathrm{SAPT0}}_{\mathrm{exch}}+\Delta
E^{\mathrm{SAPT0}}_{\mathrm{ind}}+\Delta E^{\mathrm{SAPT0}}_{\mathrm{disp}}$
(7)

We introduce the field-dependent one-particle charge density of a monomer $M$
$\rho^{M}_{\mathrm{tot}}(\mathbf{r},\mathbf{r}^{\prime},\bm{F})$ as given in
eq 8 with $\rho^{M}$ the monomer electron density and the fixed nuclear charge
density expressed in terms of the nuclear charge $Z_{I}$ of atom I and
$\delta(\mathbf{r})$ the three-dimensional Dirac delta function.

$\displaystyle\rho^{M}_{\mathrm{tot}}(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime},\bm{F})=-\rho^{M}(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime},\bm{F})+\sum_{I\in
M}Z_{I}\delta(\mathbf{r}_{1}^{\prime}-\mathbf{r}_{I})$ (8)

The first-order (in $\zeta$) contribution to the field-dependent SAPT0
interaction energy $\Delta E^{(10)}$ contains the mean field electrostatic and
exchange interaction of the two unperturbed monomer one-particle charge
densities. In the closed shell case, we can express this in a concise form as
given in eq 9, introducing the Coulomb interaction operator
$\hat{\tilde{V}}_{12}=|{\mathbf{r}_{1}-\mathbf{r}_{2}}|^{-1}\big{(}1-\frac{1}{2}\hat{P}^{\prime}_{12}\big{)}$
for the 1-2 particle index pair with the $\hat{P}^{\prime}_{12}$ permutation
operator that interchanges the two electron indices, if both present, but
evaluates to zero otherwise.

$\displaystyle\Delta E^{(10)}(\bm{F})=\Delta
E^{\mathrm{SAPT0}}_{\mathrm{elst}}(\bm{F})+\Delta
E^{\mathrm{SAPT0}}_{\mathrm{exch}}(\bm{F})=$ (9)
$\displaystyle\int_{\mathbf{r}^{\prime}_{1}=\mathbf{r}_{1}}\int_{\mathbf{r}^{\prime}_{2}=\mathbf{r}_{2}}\hat{\tilde{V}}_{12}\rho^{A}_{\mathrm{tot}}(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime},\bm{F})\rho^{B}_{\mathrm{tot}}(\mathbf{r}_{2},\mathbf{r}_{2}^{\prime},\bm{F})\mathrm{d^{3}}\mathbf{r}_{1}\mathrm{d^{3}}\mathbf{r}_{2}$

Note that the operator $\hat{\tilde{V}}_{12}$ does not act on functions of
primed variables, and the prime is removed only before integration. The
exchange interactions arise via the contributions of $\hat{P}^{\prime}_{12}$.
The one-particle electron density can be conveniently represented as matrix
$(D_{rs})$ in the atomic orbital basis set. Eq 10 introduces the the monomer
electron density matrix $(D_{rs})$.

$\rho^{M}(\mathbf{r},\mathbf{r}^{\prime},\bm{F})=\sum_{rs}D^{M}_{rs}(\bm{F})\phi^{*}_{r}(\mathbf{r}^{\prime})\phi_{s}(\mathbf{r})$
(10)

The second and higher order contributions (in $\zeta$) to the SAPT0
interaction energy contain terms from the polarization of
$\rho^{M}(\mathbf{r},\mathbf{r}^{\prime})$ due to the intermolecular
interaction. Pure induction contributions $E^{\mathrm{SAPT0}}_{\mathrm{ind}}$
(i. e. pure static polarization interaction) up to the formally infinite order
in $\zeta$ can be computed as the remaining HF contributions to the
interaction energy as given in eq 11.

$\displaystyle E^{\mathrm{SAPT0}}_{\mathrm{ind}}=\Delta E^{\mathrm{HF}}-\Delta
E^{(10)}\quad.$ (11)

The lowest order dispersion contributions (second order in $\zeta$), arise
from two-body dynamic polarization between the monomers’ charge densities and
are also taken into account in SAPT0. They contain a Coulomb interaction and
an exchange interaction part and are given in eq 12.

$\displaystyle\Delta E^{\mathrm{SAPT0}}_{\mathrm{disp}}=\Delta
E^{(20)}_{\mathrm{disp}}+\Delta E^{(20)}_{\mathrm{exch-disp}}\quad.$ (12)

The components of the SAPT0 interaction dipole moment $\Delta\mu^{(10)}$ and
interaction polarizability $\Delta\alpha_{ij}^{(10)}$ can be expressed in
terms of the monomer charge densities
$\rho^{M}_{\mathrm{tot}}(\mathbf{r},\mathbf{r}^{\prime},\bm{F})$ by inserting
eq 9 into eqs 2 and 3. Since the Coulomb interaction operator
$\hat{\tilde{V}}_{12}$ does not depend on the external field, the field
dependence of $\Delta E^{(10)}$ can be expressed in the set of atomic orbital
density matrix elements $\\{D^{A}_{rs}(\bm{F}),D^{B}_{rs}(\bm{F})\\}$, given
for the electron-electron interaction part of $\Delta E^{(10)}$ in eq 13.

$\displaystyle\Delta$ $\displaystyle E^{(10)}_{\mathrm{el-
el}}(\bm{F})=\sum_{rstu}D^{A}_{rs}(\bm{F})\ D^{B}_{tu}(\bm{F})$ (13)
$\displaystyle\int_{\mathbf{r}^{\prime}_{1}=\mathbf{r}_{1}}\int_{\mathbf{r}^{\prime}_{2}=\mathbf{r}_{2}}\phi^{*}_{r}(\mathbf{r}_{1}^{\prime})\phi_{s}(\mathbf{r}_{1})\hat{\tilde{V}}_{12}\phi^{*}_{t}(\mathbf{r}_{2}^{\prime})\phi_{u}(\mathbf{r}_{2})\mathrm{d}^{3}\mathbf{r}_{1}\mathrm{d}^{3}\mathbf{r}_{2}$

Using the chain rule and introducing the short-hand notations
$\partial_{(Mrs)}\Delta E^{(10)}\equiv\frac{\partial\Delta
E^{(10)}(\bm{F})}{\partial D^{M}_{rs}(\bm{F})}\bigr{|}_{\bm{F}=\bm{0}}$,
$\partial_{(Mrs)(Ntu)}\Delta E^{(10)}\equiv\frac{\partial^{2}\Delta
E^{(10)}(\bm{F})}{\partial D^{M}_{rs}(\bm{F})\partial
D^{N}_{tu}(\bm{F})}\bigr{|}_{\bm{F}=\bm{0}}$,
$\partial_{i}D^{M}_{rs}\equiv\frac{\partial D^{M}_{rs}(\bm{F})}{\partial
F_{i}}\bigr{|}_{\bm{F}=\bm{0}}$, and
$\partial_{ij}D^{M}_{rs}\equiv\frac{\partial^{2}D^{M}_{rs}(\bm{F})}{\partial
F_{i}\partial F_{i}}\bigr{|}_{\bm{F}=\bm{0}}$ one can express
$\Delta\mu^{(10)}_{i}$ and $\Delta\alpha^{(10)}_{ij}$ in terms of monomer-
resolved density polarizations in the electric field as given in eqs 14 and
15.

$\displaystyle\Delta\mu^{(10)}_{i}=-\sum_{rs}\big{[}\partial_{(Ars)}\Delta
E^{(10)}\partial_{i}D^{A}_{rs}+\partial_{(Brs)}\Delta
E^{(10)}\partial_{i}D^{B}_{rs}\big{]}$ (14)

$\displaystyle\Delta\alpha^{(10)}_{ij}=-\sum_{rs}\big{[}$
$\displaystyle\partial_{(Ars)}\Delta
E^{(10)}\partial_{ij}D^{A}_{rs}+\partial_{(Brs)}\Delta
E^{(10)}\partial_{ij}D^{B}_{rs}\big{]}$ (15)
$\displaystyle-\sum_{rstu}\Big{[}$ $\displaystyle\partial_{(Ars)(Btu)}\Delta
E^{(10)}\big{[}\partial_{i}D^{A}_{rs}\partial_{j}D^{B}_{tu}+\partial_{j}D^{A}_{rs}\partial_{i}D^{B}_{tu}\big{]}\
\Big{]}$

Hence, there are two separate contributions to $\Delta\mu^{(10)}_{i}$, arising
from the linear polarization of monomer A and B, respectively. There are
contributions to $\Delta\alpha^{(10)}_{ij}$ originating from quadratic
polarization on only one monomer, as well as from “mixed” contributions that
stem from the combined linear polarization of monomers A and B. Electrostatic
and exchange monomer contributions to the SAPT0 interaction-induced electric
properties can be obtained separately (cf. eq 9). The monomer decomposition
analysis of the first-order terms enables a better qualitative understanding
of the link between monomer electric properties and interaction induced
electric properties. It is in this context interesting to note that there is a
“twofold” physical interpretation of the interaction induced electric
properties,heijmen1996  e. g. for the interaction dipole moment: First, one
can see it as an expectation value of the dimer due to the mutual charge
polarization among the monomers; second, one can see it as the first order
change in the interaction energy due to the monomer charge density
polarization by the external field, as expressed in eq 14.

##  III Result and discussion

###  III.1 CCSD(T) interaction-induced properties along the z-axis

The interaction energies of the S66 set of non-covalent dimers are computed as
a function of the OEEF categorized as "H-bonding" (dimers 1-23),
"$\pi$-stacking" (dimers 24-33), "London dispersion" (dimers 34-46), and
"mixed interactions" (dimers 47-66). All dimers are depicted as 3D structures
in Figure 1 (hydrogen bonded dimers), Figure 2 ($\pi$-stacking and London
dispersion), and Figure 3 (mixed interaction). The dimers are in these Figures
arranged in matrices consisting of monomer interaction types in its rows and
columns, so that the relation between the dimers in terms of monomer
composition and conformation is visualized.

Figure 1:  Interaction energy of hydrogen bonded dimers in the S66 set as a
function of $F_{z}$ with the dimer number given in the plot background. The
plots are arranged in a matrix with rows and columns expanding the monomer
conformations in the respective combination of a dimer. Brackets denote “H
dnr” for H-bond donor, “H acpt” for H-bond acceptor, “H d & a” for H-bond
donor and acceptor. Dimer 12 not shown since it would overlap with dimer 3
within our matrix arrangement. Corresponding 3D structural depictions, linked
through black arrows, visualize the interaction dipole moment vector (red
arrow) and the interaction polarizability principal components (blue ellipsoid
with black and white rods for positive and negative eigenvalues, respectively)
computed at the dRPA@PBE0/aug-cc-pVQZ level of theory. The 3D structures were
created using the PyMol software package.PyMOL  Figure 2:  Interaction energy
of $\pi$-stacking and London dispersion bonded dimers in the S66 set as a
function of $F_{z}$ with the dimer number given in the plot background. The
plots are arranged in a matrix with rows and columns expanding the monomer
conformations in the respective combination of a dimer. Brackets denote “side
on” and “end on” for parallel and orthogonal relative placements,
respectively. Corresponding 3D structural depictions, linked through black
arrows, visualize the interaction dipole moment vector (red arrow) and the
interaction polarizability principal components (blue ellipsoid with black and
white rods for positive and negative eigenvalues, respectively) computed at
the dRPA@PBE0/aug-cc-pVQZ level of theory.  Figure 3:  Interaction energy of
mixed interaction bonded dimers in the S66 set as a function of $F_{z}$ with
the dimer number given in the plot background. The plots are arranged in a
matrix with rows and columns expanding the monomer conformations in the
respective combination of a dimer. Brackets denote “H dnr” for H-bond donor,
“H acpt” for H-bond acceptor, “H d & a” for H-bond donor and acceptor.
Corresponding 3D structural depictions, linked through black arrows, visualize
the interaction dipole moment vector (red arrow) and the interaction
polarizability principal components (blue ellipsoid with black and white rods
for positive and negative eigenvalues, respectively) computed at the
dRPA@PBE0/aug-cc-pVQZ level of theory.

According to our CCSD(T) results, the magnitude of the first interaction
hyperpolarizability along the $z$-axis is usually small (0.165 Å5 on average
obtained by a polynomial fit up to third order which means a 0.002 kcal mol-1
effect in the interaction energy at our largest field strength). Hence, the
interaction energy as a function of the OEEF is fully determined by the zero
field interaction energy, interaction dipole moment, and interaction
polarizability. For further details on the numerical fitting procedure we
refer to the SI.

The reference CCSD(T) interaction energies, dipole moments, and
polarizabilities along the $z$-axis are shown in Figure 4. Among the
interaction types, the magnitude of the interaction energy is large in the
case of hydrogen bonding ($\approx$9 kcal mol-1 on average), while moderately
large in the other cases ($\approx$3-4 kcal mol-1 on average).

We observe positive interaction dipole moments in most of the cases, which
means the interaction-induced polarization enhances the monomer polarity.
However, negative interaction dipole moments occur for the hydrogen bonded
AcOH - uracil dimer, and for some of the dimers corresponding to the London
dispersion and mixed interaction classes. Furthermore, (near) zero interaction
dipole moments occur for some homodimers due to symmetry (or near-symmetry).
The largest magnitudes of interaction dipole moments occur for hydrogen bonded
systems ($\approx$0.7 D on average), dispersion bonded cases show rather small
magnitudes ($\approx$0.1 D on average), and mostly moderate magnitudes are
found in the mixed influenced cases ($\approx$0.4 D on average). The large
interaction dipole moment for asymmetric hydrogen bonded systems (e. g. dimers
1, 15, and 17) or hydrogen-$\pi$-interacting systems (e. g. dimers 47 and 57)
suggests that the OEEF can considerably affect the interaction energy in ice,
proteins, or double stranded DNA.
Parallel displaced and T-shaped $\pi$-stacking homodimers (e. g. dimers 24 and
47) have similar interaction energies. The former however have a near zero
interaction dipole moment along the $z$-axis while the latter can bear a
significant interaction dipole moment. This may be utilized for selective
growth of molecular crystals with fishbone lattice structures from melts in
strong OEEFs.
In most of the cases, the magnitude of the interaction polarizability is small
($\approx$0.6-0.7 Å3 on average), and it has only a small effect on the shape
of the interaction energy function with respect to the field strength. Its
influence is somewhat more important in the $\pi$-$\pi$-interaction and London
dispersion dominated dimers ($\approx$0.9 Å3 on average), where the
interaction dipole moment is small. The largest interaction polarizabilities
occur in dimers of aromatic heterocycles with the ring planes oriented along
the $z$ axis (dimers 17 and 58) with an increase in interaction energy of up
to 0.87 kcal mol-1 at the largest considered field strength.

Figure 4: CCSD(T)-F12b/aug-cc-pVDZ interaction energies, dipole moments, and
polarizabilities for the S66 noncovalent dimers along the z-axis. Vertical
black dashed lines separate hydrogen bonding, $\pi$-stacking, London
dispersion, and mixed interaction classes.

###  III.2 SAPT0 decomposition of the interaction-induced properties along the
z-axis

For a qualitative understanding of the most relevant physical interactions, we
compute the SAPT0 decomposition of the interaction induced electric
properties. Interaction induced properties of Watson-Crick amino acid pairs
have also been studied with higher orders of SAPTczyznikowska2013  and also
with taking into account nuclear relaxation effects for small hydrogen- and
halogen-bonded dimers in the order up to interaction first
hyperpolarizabilities.medved2020 ; zalensky2018  The SAPT0 mean absolute
errors for the interaction energy, dipole moment, and polarizability are,
respectively, 0.50 kcal mol-1, 0.039 D, and 0.15 Å3. Hence, SAPT0 is in
qualitative agreement with the coupled cluster reference.
Figure 5 depicts the SAPT0 electrostatic, exchange, induction, and dispersion
contributions for $\Delta E$, $\Delta\mu_{z}$, and $\Delta\alpha_{z}$. In the
case of hydrogen bonding, all four interaction contributions to the
interaction energy are significant. The sum of the (first-order) electrostatic
and exchange parts is usually attractive; however, both the induction and
dispersion parts are larger in magnitude than the first-order interaction
energy. In case of dispersion bonded dimers the first-order interaction energy
is repulsive and the dimers are bound by the dominating attractive dispersion
interaction.
Interestingly, the contributions to the interaction dipole moment strongly
differ from those to the interaction energy. In the hydrogen bonded dimers,
the interaction dipole moment is always positive and dominated by the
electrostatic part. Exchange contributions are rather small and tend to
slightly moderate the electrostatic contribution. Higher-order static
polarization effects described by the induction contributions further enhance
the interaction dipole moment. Intermolecular dispersion at the SAPT0 level of
theory gives very small contributions. In the $\pi$-stacking and London
dispersion bonded dimers, the electrostatic part of the interaction dipole
moment is quenched by the exchange part to some extent, and the direction of
the interaction dipole moment does not seem to correlate with the direction of
the sum of the monomer dipole moments. The induction and intermolecular
dispersion contributions are small for SAPT0. It should be pointed out here
that SAPT0 does not cover higher-order terms that include many-body and
intramonomer dispersion. For the interaction polarizability of cytosine dimers
in stacked alignments, it was indeed found that electrostatic effects from
intramonomer correlation are even larger than the lowest order intermolecular
dispersion contributions.czyznikowska2013
Regardless the interaction class, the interaction polarizability contains
important static polarization contributions in the lowest order
(electrostatics and exchange) as well as at higher orders (induction). The
SAPT0 dispersion contribution is found to be only significant in the case of
$\pi$-$\pi$-interactions. Czyżnikowska _et al._czyznikowska2013  found for
Watson-Crick amino acid pairs very strong exchange and induction contributions
to the interaction polarizability which underlines the exceptionally strong
interaction in these dimer systems. Our findings are in agreement with an
earlier computation for the HCN dimer where the interaction polarizability was
found to be greatly influenced by electrostatic
effects.dft_interaction_induced3

Figure 5: SAPT0 decomposition of the interaction energy, dipole moment, and
polarizability for the S66 noncovalent dimers along the z-axis. Vertical black
dashes lines separate hydrogen bonding, $\pi$-stacking, London dispersion, and
mixed interaction classes.

###  III.3 SAPT0 monomer contributions to the first-order interaction-induced
properties along the z-axis

A link between monomer electric properties and interaction-induced electric
properties is provided for the SAPT0 electrostatic and exchange energy
contributions by eqs 14 and 15. The monomers are assigned the labels “A” and
“B” as indicated in the 3D structures of the dimers in Figure 1 (hydrogen
bonded dimers), Figure 2 ($\pi$-stacking and London dispersion), and Figure 3
(mixed interaction). Figures 6 and 7 show the monomer contributions to the
interaction dipole moment and the interaction polarizability, respectively. In
singly hydrogen bonded dimers (dimers 1-16 and 18-19), the hydrogen bond donor
(monomer A) and acceptor (monomer B) contributions to the electrostatic
interaction dipole moment have the same sign. The donor contribution to the
electrostatic interaction dipole moment is relatively less pronounced in the
case of OH groups (31% on average with the water and methanol H-bond donors),
while somewhat more pronounced in the case of NH groups (53% on average with
the peptide and methyl amine H-bond donors). In (nearly or fully) symmetric
doubly hydrogen bonded dimers (dimers 17 and 20–23) or $\pi$-stacking/London
dispersion dominated homodimers (dimers 24-26, 34, 36, and 38), the oppositely
signed monomer contributions to the electrostatic interaction dipole moment
mostly cancel each other. The directionality of hydrogen bonds is thus a
strong vector for designing the first-order interaction electric response. The
monomer contributions to the exchange interaction dipole moment always have
the opposite sign and approximately (or exactly, in case of symmetric
homodimers) the same magnitude. This is a manifestation of the increase in
Pauli repulsion when the monomer electron density polarizes towards the other
monomer. As this results generally in only a small net exchange contribution
to the interaction dipole moment, Pauli repulsion appears to act in a rather
“symmetric” way upon monomer polarization. This observation holds for all
dimer systems with the exception of some $\pi$-stacking dimers (dimers 28–33).
It is apparent that the interaction-induced static polarization of the
hydrogen bond donor and acceptor is the main contribution to large interaction
dipole moments in the hydrogen-bonded cases. The hydrogen donor-acceptor
conformation in these dimers leads to a strong electrostatic contribution to
the interaction energy (cf. Figure 5) and rather short dimer distances where
the charge distributions of the monomers arrange optimally one to each other.
This plays in favour of a strong monomer electrostatic polarization. For the
systems with $\pi$-$\pi$ interactions and London dispersion interaction, the
picture is more heterogeneous regarding the electrostatic monomer
contributions; they may be of the same or opposite sign and are quite large in
magnitude for some cases.
The electrostatic interaction polarizability is composed only of positive
contributions for both monomer second-order polarization and combined monomer
first-order polarization. For the monomer second-order polarizations this is
an interesting finding, as it suggests that the monomer static polarizability
always increases along the axis of interaction, i. e. the monomer polarization
due to the interaction always leads to a more polarizable monomer electron
density. These monomer second-order polarizations are especially strong in
hydrogen-bonded and $\pi$-stacking dimers. It is remarkable that monomer first
hyperpolarizabilites apparently play an important role, as they are at the
origin of the monomer second order polarization, while interaction first
hyperpolarizabilities are observed to be insignificant. However, in the case
of exceptionally “field-sensitive” systems, such as electron-rich, polar,
strongly hydrogen-bonded dimers that are aligned along the $z$-axis,
significant interaction first hyperpolarizabilities have been
predicted.gora2013
In London dispersion and mixed interaction dimers the combined monomer
polarization dominates mostly and is large for all dimers. This is to be
expected, as the “cooperative” polarization of both monomers under the
electric field $F_{z}$ increases the charge density mobility across the dimer,
as compared to the isolated monomer charge density mobility. Remarkably, the
monomer coupling term in the electrostatic interaction polarizability
$\Delta\alpha_{zz}$ correlates well with the product of the corresponding
monomer polarizabilities ($R^{2}\approx 0.92$, see the SI for numerical
details). This proportionality may be a suitable simple design principle for
the second-order behavior of the interaction energy in OEEFs in direction of
the interaction axis $z$.
The individual monomer contributions to the exchange interaction
polarizability dominate by far and are negative. This can be understood by
considering that the repulsive intermolecular exchange interaction diminishes
the polarization capacity of the individual monomer charge densities. The
monomer coupling contributions to the exchange interaction polarizability give
small positive contributions, indicating that slight cooperative effects in
monomer polarization capacity may also be mediated by the exchange
interactions.

Figure 6: Monomer first-order polarization contributions to the SAPT0
electrostatic (top panel) and exchange interaction (bottom panel) dipole
moments for the S66 noncovalent dimers along the z-axis. Vertical black dashed
lines separate hydrogen bonding, $\pi$-stacking, London dispersion, and mixed
interaction classes. Figure 7: Monomer second-order polarization and combined
monomer first-order polarization contributions to the SAPT0 electrostatic and
exchange interaction polarizabilities for the S66 noncovalent dimers along the
z-axis. Vertical black dashed lines separate hydrogen bonding, $\pi$-stacking,
London dispersion, and mixed interaction classes.

###  III.4 Performance of density functional approximations for interaction
induced electric properties

The shortcomings (and corrections) of many of the original density functional
approximations for van der Waals energies are well established. Their
performance for interaction dipole moments and polarizabilities, however, is
not that known. While the accuracy of calculated electrostatic dipole (and
multiple) moments of monomers is typically quite decent ChemistsGuidetoDFT ,
the predictive power of density functionals can vary dramatically when it
comes to polarizability van1996improved ; champagne1998assessment ;
willow2016mp2 ; hait2018accurate ; akter2021well . In order to facilitate
future studies dealing with improved density functional approximations, we
therefore have also included a selection of common DFT approximations in this
study, and we now discuss the corresponding numerical estimates of van der
Waals interaction energies, dipole moments, and polarizabilities.

We have placed our main focus on results related to the PBE semilocal
exchange-correlation kernel, augmented by various nonlocal exchange and
correlation ingredients. In terms of methodology aiming at the quantitative
predictive power needed to account properly for long-range dispersion, we
decided to focus our DFT benchmark study on this latter part. Figure 8 shows
the mean absolute percentage error (MAPE) of the considered density
functionals, including purely semilocal PBE and PBE0 for comparison. Semilocal
functionals like PBE cannot capture long-range dispersion interactions due the
nonlocal nature of the latter. The D3 dispersion correction can give a
satisfactory description of long-range dispersion in some use cases, but is
insensitive to the electric field. The VV10 dispersion correction is based on
the electron density, but its effect on the interaction dipole moments is
quite small, and the interaction polarizabilities can become even more
erroneous if the calculated electron density is too diffuse. The exact
exchange mixing in hybrid functionals significantly lowers the interaction
dipole moments and polarizabilities errors. This may indicate that the quality
of the calculated electron density is better compared to the electron
densities obtained by purely semilocal functionals (the latter being too
delocalized and diffuse due to the large self-interaction error). However, an
additional dispersion correction is clearly needed for reliable predictions.
The interaction-induced properties obtained with the screened hybrid HSE06
functional are very similar to the ones with the global hybrid PBE0 with the
same fraction of exact exchange at short-range, due to the weak screening by
the very small range-separation parameter. The long-range correction in the
LC-$\omega$PBE functional considerably improves the interaction dipole moments
and polarizabilities, but such functionals still require a dispersion
correction. Our observations up to this point are in agreement with the
literature, where the long-range corrected LC-BLYP methodb88 ; lyp ; lc_scheme
was reported to yield more accurate interaction-induced properties than the
semilocal BLYP methodb88 ; lyp  for hydrogen bonded
dimers,dft_interaction_induced3  and global hybrid density functionals were
found to be qualitatively correct for interaction dipole moments and mean
interaction polarizabilities, but not for the anisotropy of the interaction
polarizabilities for noble gas dimers.dft_interaction_induced4  Our study
suggests that long-range corrected functionals are generally more accurate
than global hybrid ones for interaction-induced properties even though global
hybrids can be occasionally better than long-range corrected hybrids depending
on the parent semilocal functional.dft_interaction_induced2  The double hybrid
PBE0-2 functional with partial second-order perturbative correlation and a
large fraction of exact exchange delivers somewhat less accurate interaction
energies but more accurate interaction dipole moments and polarizabilities.
The dual-hybrid dRPA75 functional with full dRPA correlation and a similarly
large fraction of exact exchange yields similarly accurate results. Finally,
the dRPA method with full exact exchange and dRPA correlation (dRPA@PBE0)
leads to very accurate results for all three quantities as compared to the
coupled cluster reference (MAPE$\leq$3% for all three quantities).
Nevertheless it should be noted that the dRPA interaction dipole moments and
polarizabilities are quite sensitive to the reference orbitals: Omitting the
exact exchange (dRPA@PBE) introduces considerable errors to the interaction
polarizabilities. The sensitivity of density functional approximations with
respect to orbital relaxation and higher-order long-range dispersion effects
was also found for the computation of electric properties of endohedral small
polar molecules in carbon nanotubes.kozlowska

In summary, the LC-$\omega$PBE functional is a surprisingly reliable method
given its low computational cost (MAPE$\leq$5% for all three quantities), but
other computationally efficient methods (PBE0-D3, PBE0-VV10, HSE06-D3) also
perform fairly well (MAPE$\leq$10% for all three quantities).

Figure 8:  Mean absolute percentage error (MAPE) of the interaction energies,
dipole moments, and polarizabilities obtained with various density functional
approximations for the S66 noncovalent dimers along the z-axis relative to the
average magnitude of the same properties calculated with the CCSD(T)-F12b/aug-
cc-pVDZ reference method.

###  III.5 Interaction-induced properties in three dimensions

Given the good performance of dRPA@PBE0/aug-cc-pVQZ with respect to the
coupled cluster reference, this level of theory was chosen to compute all the
elements of the interaction dipole moment vector and interaction
polarizability matrix with a three-dimensional grid (as described in the
methods section) that allows for fitting of a quadratic function. The
interaction dipole moments are depicted in vector form (direction and
magnitude) for the S66 dimers in Figure 1 (hydrogen bonded dimers), Figure 2
($\pi$-stacking and London dispersion), and Figure 3 (mixed interaction). The
interaction polarizabilities with directions and magnitude of the three
principal components depicted as ellipsoids are given in the same Figures. We
tested the quadratic dependency of the interaction energy with respect to the
OEEF for the CCSD(T) reference method also in the xy-plane on a denser grid
for some selected dimers. The magnitude of the first interaction
hyperpolarizability components is in this case 0.046 Å5 on average, obtained
by a polynomial fit up to third order and resulting in a 0.0007 kcal mol-1
effect on the interaction energy at our largest field strength (further
details are given in the SI). Hence, the first interaction hyperpolarizability
appears also to be negligible in the $xy$-plane.

The average value of the angle of the interaction dipole moment vector with
respect to the z-axis is 28.0° for the whole database, while it is 12.6° for
hydrogen bonding, 38.4° for $\pi$-$\pi$-interaction, and 48.4° for London
dispersion dominated dimers. The individual angles for the S66 dimers are
given in Figure 9. The dominant interaction dipole moment component along the
$z$-axis for H-bonded systems is a result of the hydrogen bond’s
directionality. In the case of $\pi$-$\pi$\- and London dispersion
interactions the orientation of the interaction dipole moment mainly
originates from the monomers’ orientation, as the lowest-order contribution to
the $\Delta\mu_{i}$ vector component depends on the monomer polarizations
$\partial_{i}D^{M}_{rs}$ (see eq 14).

Figure 9:  Interaction dipole moments for the S66 noncovalent dimers
calculated with dRPA. Magnitude of the interaction dipole moment vector (top),
interaction dipole moment component perpendicular to the z-axis (mid), and
angle of the interaction dipole moment vectors with respect to the z-axis
(bottom). Vertical black dashed lines separate hydrogen bonding,
$\pi$-stacking, London dispersion, and mixed interaction classes.

The principal axis corresponding to the first principal component includes a
small angle (0-22.8°) with the z-axis and is shown in Figure 10. This angle is
particularly small for London dispersion bonded dimers (3.3° on average),
while somewhat larger for hydrogen bonded (8.8° on average) and
$\pi$-$\pi$-interaction bonded (7.1° on average) dimers. The first principal
component ($\bm{\nu}^{c}$) of the interaction polarizability has a positive
eigenvalue. The second ($\bm{\nu}^{b}$) and third ($\bm{\nu}^{a}$) principal
components of the interaction polarizability have negative eigenvalues.
$\bm{\nu}^{b}$ and $\bm{\nu}^{a}$ are smaller than $\bm{\nu}^{c}$ in the
hydrogen bonded cases, indicating the dominance of the directionality imposed
by the hydrogen bonds. In the case of London dispersion and mixed
interactions, the relative magnitude of the components depends on the relative
orientation of the monomers. This is in line with our observations on how the
monomers’ linear and quadratic polarization components strongly determine the
interaction polarizability (compare for example dimers 17 and 26, dimers 30
and 40, or dimers 34 and 36 in Figures 7 and 10).

The sign of the interaction polarizability principal components is mostly in
line with the intuitive assumption that dimer polarization parallel to the
$z$-axis aligns dipolar polarization at opposite poles and increases the
intermolecular bonding while dimer polarization in the $xy$-plane aligns
dipolar polarization at identical poles and decreases the intermolecular
bonding. Even though this simple and intuitive design principle for the
modulation of the intermolecular interaction via a dipolar field is here
confirmed for the S66 set, very strong fields would be required for its
successful application.

Figure 10: The angle of the interaction polarizability first principal axis
with respect to the z-axis as well as the three principal components of the
interaction polarizability for the S66 noncovalent dimers calculated with
dRPA@PBE0 (“c” first principal axis, “b” second principal axis, “a” third
principal axis). Vertical black dashed lines separate hydrogen bonding,
$\pi$-stacking, London dispersion, and mixed interaction classes.

Figures 1 (hydrogen bonded dimers), Figure 2 ($\pi$-stacking and London
dispersion), and Figure 3 (mixed interaction) depict the quadratic dependency
of the interaction energy with respect to our employed field strength of up to
$\approx$5 GV m-1 along the $z$-axis. It visualizes the range of modulation of
the interaction energies up to a moderately strong OEEF. The range of
modulation of the interaction energies by an OEEF in the $xy$-plane can be
obtained by inspection of the interaction dipole moment arrows and interaction
polarizability ellipsoids in the Figures and also by comparing with Figures 9
and 10.
Monomers with more extended $\pi-$conjugations than those present in the S66
set (for example fused rings) generally show a larger charge density
polarizability and cooperative polarization coupling via many-body
effectskozlowska ; zawada2011  and are therefore capable of larger modulation
effects of interaction energies by an OEEF.

###  III.6 Intermolecular distance dependence of the interaction induced
properties

The effect of the OEEF on the noncovalent interaction potential energy curves
of the S66 dimers is briefly discussed in this section. Table 1 depicts values
of the averaged interaction energy, dipole moment ($z$-component), and
polarizability ($z$-component) discretized for four intermolecular distances
$r/r_{\mathrm{e}}$ ($r_{\mathrm{e}}$ being the equilibrium distance) along the
$z$-axis. The interaction energy is given for zero field strengths, as well as
for our strongest employed field strengths ($\pm\ \approx$ 5.14 109 V m-1). On
average, and in line with the aforementioned concavity with the field, the
interaction strengthens upon application of the field in either direction.
This observation also holds for the London dispersion dominated dimers (see
also Fig. 2). For the hydrogen bonded dimers (shown in Fig. 1), interaction
energies are substantially larger in magnitude, and their interaction strength
can be either decreased or weakened, depending on the field’s direction. As
one would expect (due to the cancellation of the nuclear repulsion terms), the
magnitude of the interaction induced electric properties decreases
monotonically with the intermolecular distance. Not surprisingly, the
interaction dipole moment in hydrogen bonded dimers is largest in magnitude
and still quite sizeable at twice the equilibrium distance. By contrast,
interaction dipole moments are negligible for the London dispersion dominated
dimers. At such a small absolute scale, we think it unlikely that their
changing sign for some distances carries much meaning due to all numbers in
Table 1 corresponding to averages over dozens of chemistries. For interaction
polarizabilities, an interesting trend is observed: London dispersion
dominated dimers exhibit significantly larger values than hydrogen bonded
dimers for short distances, at twice the equilibrium distance, however, the
interaction polarizability of the hydrogen bonded dimers is larger. Presumably
this is due to London dispersion decaying with a faster power law than
electrostatic interactions.

Table 1:  Interaction energies $\Delta E$ [kcal/mol], dipole moments $\Delta\mu$ [D], and polarizabilities $\Delta\alpha$ [Å3] at CCSD(T)-F12b/aug-cc-pVDZ level for four intermolecular distances, averaged over the entire S66 dimer sets, as well as for dimer sub sets HB and $\pi$ corresponding to being bonded predominantly through hydrogen bonding (shown in Fig. 1) or London dispersion (shown in Fig. 2), respectively. Averaged interaction energies $\Delta E^{n}$, $\Delta E^{z}$, and $\Delta E^{p}$ correspond to results obtained in electric fields with negative, zero, and positive strengths $F_{z}=$ -0.01, 0.00, and +0.01 a.u., respectively.  $r/r_{\mathrm{e}}$ | 9/10 | 1 | 5/4 | 2
---|---|---|---|---
$\Delta E^{n}$ | -5.059 | -5.927 | -4.107 | -0.915
$\Delta E^{z}$ | -4.556 | -5.487 | -3.786 | -0.790
$\Delta E^{p}$ | -4.679 | -5.630 | -3.925 | -0.855
$\Delta E^{n}_{\rm HB}$ | -10.093 | -10.451 | -7.641 | -2.003
$\Delta E^{z}_{\rm HB}$ | -8.322 | -8.986 | -6.684 | -1.677
$\Delta E^{p}_{\rm HB}$ | -7.073 | -8.024 | -6.158 | -1.555
$\Delta E^{n}_{\pi}$ | -2.438 | -3.796 | -2.310 | -0.322
$\Delta E^{z}_{\pi}$ | -2.266 | -3.609 | -2.142 | -0.252
$\Delta E^{p}_{\pi}$ | -2.820 | -4.088 | -2.477 | -0.375
$\Delta\mu$ | 0.356 | 0.292 | 0.189 | 0.067
$\Delta\mu_{\rm HB}$ | 0.823 | 0.671 | 0.427 | 0.147
$\Delta\mu_{\pi}$ | 0.004 | -0.001 | -0.002 | 0.001
$\Delta\alpha$ | 0.740 | 0.688 | 0.541 | 0.225
$\Delta\alpha_{\rm HB}$ | 0.615 | 0.595 | 0.510 | 0.241
$\Delta\alpha_{\pi}$ | 0.858 | 0.786 | 0.595 | 0.228

##  IV Conclusion

Interaction induced electric properties along the binding axis for the diverse
dimers of the S66 set have been presented using high-level coupled cluster
computations [CCSD(T)-F12b/aug-cc-pVDZ-F12] as a reference. Our results
suggest that the intermolecular energy’s response to variations in the
electric field along the intermolecular axis is largely dominated by first
(dipole) and second (polarizability) order effects, while higher order
contributions (hyperpolarizabilities) play only a negligible role. Based on
symmetry adapted perturbation theory (SAPT0), a qualitative assessment has
been performed for the most relevant physical interactions that govern the
observed interaction induced properties. Analysis of first order SAPT
contributions computed on a monomer-resolved level reveals some essential
links between composing monomer and dimer properties. The reference results
have served as a benchmark for testing the performance of frequently employed
density functional approximations (DFAs) of non-local exchange and correlation
in combination with the semi-local PBE exchange-correlation kernel and the
aug-cc-pVQZ basis set. We find that exact exchange contributions play an
important role for the performance of the DFAs, raw PBE showing an unsigned
relative mean error above 10%. Among the computationally more efficient DFAs,
the LC-$\omega$PBE0-D3 level of theory stands out providing reliable values
for interaction dipole moments and dipole polarizabilities throughout
(relative error below 5%). A very small relative error of at most a few
percent is obtained at the computationally more demanding dRPA@PBE0 level of
theory. The dRPA@PBE0 level of theory has been used to compute the full
interaction dipole moment vector, as well as the principal components of the
interaction polarizability tensors. Our tensorial property analysis indicates
that the vector of directionality given by hydrogen bonds for manipulating
dominates the first order interaction response. It also reveals that the
principal components of the interaction polarizability stay well aligned with
the interaction axis in all cases (at most $\approx$20∘ deviation), and
consistently offer one vector for increasing the interaction energy
(approximately along the binding axis) and two perpendicular vectors for
decreasing the interaction energy. Results for intermolecular distance
dependent averages of interaction energies, dipole moments, and polarizability
also suggest interesting trends for hydrogen bonded vs. London dispersion
bonded systems. In particular, we find that fields can only strengthen
interaction energies for London dispersion bonded systems, while they can
strengthen and weaken the hydrogen bonded systems (depending on direction).
Interaction dipole moments decay monotonically with distance for hydrogen
bonded systems, and are negligible for London dispersion bound systems.
Interaction polarizability trends depend on distance: In the equilibrium
distance they are larger for London dispersion bound systems, for larger
intermolecular distances hydrogen bonded systems exhibit larger interaction
polarizabilities. We believe that our study of interaction induced electric
properties, focusing on accurate reference data, DFA non-local exchange-
correlation performance testing, diverse interaction types, and qualitative
interpretation of the most relevant physical interactions, bears essential
exploratory insights for future endeavours in intermolecular manipulation and
design via oriented external electric fields.

##  V Supporting information

A spreadsheet containing numerical results for all properties discussed in
this work is included.

##  VI Acknowledgment

We acknowledge support from the European Research Council (ERC-CoG grant QML).
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under Grant Agreement #772834. This work was
partly supported by the NCCR MARVEL, funded by the Swiss National Science
Foundation. We also acknowledge support by the Swiss National Science
foundation (No. PP00P2_138932, 407540_167186 NFP 75 Big Data, 200021_175747,
NCCR MARVEL). Some calculations were performed at sciCORE
(http://scicore.unibas.ch/) scientific computing core facility at University
of Basel.

## References

  * (1) M. A. Spackman, P. Munshi, and D. Jayatilaka, “The use of dipole lattice sums to estimate electric fields and dipole moment enhancement in molecular crystals,” Chem. Phys. Lett., vol. 443, no. 1-3, pp. 87–91, 2007.
  * (2) S. D. Fried and S. G. Boxer, “Electric fields and enzyme catalysis,” Annu. Rev. Biochem., vol. 86, pp. 387–415, 2017.
  * (3) Y. Hong, J. Hahn, and H. Kang, “Electron transfer through interfacial water layer studied by scanning tunneling microscopy,” J. Chem. Phys., vol. 108, no. 11, pp. 4367–4370, 1998.
  * (4) Y. Pocker and R. F. Buchholz, “Electrostatic catalysis of ionic aggregates. I. Ionization and dissociation of trityl chloride and hydrogen chloride in lithium perchlorate-diethyl ether solutions,” J. Am. Chem. Soc., vol. 92, no. 7, pp. 2075–2084, 1970.
  * (5) P. J. Smith and C. S. Wilcox, “The chemistry of functional group arrays. Electrostatic catalysis and the “intramolecular salt effect”.,” Tetrahedron, vol. 47, no. 14, pp. 2617 – 2628, 1991.
  * (6) A. A. Arabi and C. F. Matta, “Effects of external electric fields on double proton transfer kinetics in the formic acid dimer,” Phys. Chem. Chem. Phys., vol. 13, no. 30, pp. 13738–13748, 2011.
  * (7) R. Meir, H. Chen, W. Lai, and S. Shaik, “Oriented electric fields accelerate Diels–Alder reactions and control the endo/exo selectivity,” Chem. Phys. Chem., vol. 11, no. 1, pp. 301–310, 2010.
  * (8) S. Shaik, S. P. de Visser, and D. Kumar, “External electric field will control the selectivity of enzymatic-like bond activations,” J. Am. Chem. Soc., vol. 126, no. 37, pp. 11746–11749, 2004.
  * (9) T. Thuren, A. P. Tulkki, J. A. Virtanen, and P. K. Kinnunen, “Triggering of the activity of phospholipase A2 by an electric field,” Biochem., vol. 26, no. 16, pp. 4907–4910, 1987.
  * (10) M. Alemani, M. V. Peters, S. Hecht, K.-H. Rieder, F. Moresco, and L. Grill, “Electric field-induced isomerization of azobenzene by STM,” J. Am. Chem. Soc., vol. 128, no. 45, pp. 14446–14447, 2006.
  * (11) A. C. Aragonès, N. L. Haworth, N. Darwish, S. Ciampi, N. J. Bloomfield, G. G. Wallace, I. Diez-Perez, and M. L. Coote, “Electrostatic catalysis of a Diels–Alder reaction,” Nature, vol. 531, no. 7592, pp. 88–91, 2016\.
  * (12) S. Shaik, D. Mandal, and R. Ramanan, “Oriented electric fields as future smart reagents in chemistry,” Nat. Chem., vol. 8, no. 12, p. 1091, 2016\.
  * (13) C. Foroutan-Nejad and R. Marek, “Potential energy surface and binding energy in the presence of an external electric field: modulation of anion–$\pi$ interactions for graphene-based receptors,” Phys. Chem. Chem. Phys., vol. 16, pp. 2508–2514, 2014.
  * (14) S. Shaik, D. Danovich, J. Joy, Z. Wang, and T. Stuyver, “Electric-Field Mediated Chemistry: Uncovering and Exploiting the Potential of (Oriented) Electric Fields to Exert Chemical Catalysis and Reaction Control,” J. Am. Chem. Soc., vol. 142, no. 29, pp. 12551–12562, 2020.
  * (15) C. A. Dietrich, R. Schuldt, D. Born, H. Solodenko, G. Schmitz, and J. Kästner, “Evaporation and Fragmentation of Organic Molecules in Strong Electric Fields Simulated with DFT,” J. Phys. Chem. A, vol. 124, no. 41, pp. 8633–8642, 2020.
  * (16) S. Grimme, H. Kruse, L. Goerigk, and G. Erker, “The mechanism of dihydrogen activation by frustrated Lewis pairs revisited,” Angew. Chem. Int. Ed., vol. 49, no. 8, pp. 1402–1405, 2010.
  * (17) H. Lehle, J. M. Kriegl, K. Nienhaus, P. Deng, S. Fengler, and G. U. Nienhaus, “Probing electric fields in protein cavities by using the vibrational Stark effect of carbon monoxide,” Biophys. J., vol. 88, no. 3, pp. 1978–1990, 2005.
  * (18) I. T. Suydam, C. D. Snow, V. S. Pande, and S. G. Boxer, “Electric fields at the active site of an enzyme: Direct comparison of experiment with theory,” Science, vol. 313, no. 5784, pp. 200–204, 2006.
  * (19) S. D. Fried and S. G. Boxer, “Measuring electric fields and noncovalent interactions using the vibrational Stark effect,” Acc. Chem. Res., vol. 48, no. 4, pp. 998–1006, 2015.
  * (20) A. Gopher, Y. Blatt, M. Schönfeld, M. Okamura, G. Feher, and M. Montal, “The effect of an applied electric field on the charge recombination kinetics in reaction centers reconstituted in planar lipid bilayers,” Biophys. J., vol. 48, no. 2, pp. 311–320, 1985.
  * (21) Z. Popovic, G. Kovacs, P. Vincett, G. Alegria, and P. Dutton, “Electric field dependence of recombination kinetics in reaction centers of photosynthetic bacteria,” Chem. Phys., vol. 110, no. 2-3, pp. 227–237, 1986.
  * (22) S. Franzen, R. F. Goldstein, and S. G. Boxer, “Electric field modulation of electron transfer reaction rates in isotropic systems: long distance charge recombination in photosynthetic reaction centers,” J. Phys. Chem., vol. 94, no. 12, pp. 5135–5149, 1990.
  * (23) P. M. Moroney, T. A. Scholes, and P. C. Hinkle, “Effect of membrane potential and pH gradient on electron transfer in cytochrome oxidase,” Biochem., vol. 23, no. 21, pp. 4991–4997, 1984.
  * (24) R. R. Knowles and E. N. Jacobsen, “Attractive noncovalent interactions in asymmetric catalysis: Links between enzymes and small molecule catalysts,” Proc. Natl. Acad. Sci., vol. 107, no. 48, pp. 20678–20685, 2010.
  * (25) P. Metrangolo and G. Resnati, “Halogen versus hydrogen,” Science, vol. 321, no. 5891, pp. 918–919, 2008.
  * (26) S. Zhang, “Fabrication of novel biomaterials through molecular self-assembly,” Nat. Biotechnol., vol. 21, no. 10, p. 1171, 2003.
  * (27) C. Parks, A. Koswara, H.-H. Tung, N. Nere, S. Bordawekar, Z. K. Nagy, and D. Ramkrishna, “Molecular dynamics electric field crystallization simulations of paracetamol produce a new polymorph,” Cryst. Growth Des., vol. 17, no. 7, pp. 3751–3765, 2017.
  * (28) A. Kleshchonok and A. Tkatchenko, “Tailoring van der Waals dispersion interactions with external electric charges,” Nat. Commun., vol. 9, no. 1, p. 3017, 2018.
  * (29) A. Ambrosetti and P. L. Silvestrelli, “Hidden by graphene–towards effective screening of interface van der waals interactions via monolayer coating,” Carbon, vol. 139, pp. 486–491, 2018.
  * (30) J. Kozłowska, M. Schwilk, A. Roztoczyńska, and W. Bartkowiak, “Assessment of DFT for endohedral complexes’ dipole moment: PNO-LCCSD-F12 as a reference method,” Phys. Chem. Chem. Phys., vol. 20, no. 46, pp. 29374–29388, 2018.
  * (31) M. Schwilk, Q. Ma, C. Köppl, and H.-J. Werner, “Scalable electron correlation methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNO-LCCSD),” J. Chem. Theory Comput., vol. 13, pp. 3650–3675, 2017.
  * (32) Q. Ma, M. Schwilk, C. Köppl, and H.-J. Werner, “Scalable electron correlation methods. 4. Parallel explicitly correlated local coupled cluster with pair natural orbitals (PNO-LCCSD-F12),” J. Chem. Theory Comput., vol. 13, pp. 4871–4896, 2017.
  * (33) J. Kozłowska, P. Lipkowski, A. Roztoczyńska, and W. Bartkowiak, “DFT and spatial confinement: a benchmark study on the structural and electrical properties of hydrogen bonded complexes,” Phys. Chem. Chem. Phys., vol. 21, pp. 17253–17273, 2019.
  * (34) R. Zaleśny, M. Garcia-Borràs, R. W. Góra, J. M. Luis, et al., “On the physical origins of interaction-induced vibrational (hyper) polarizabilities,” Phys. Chem. Chem. Phys., vol. 18, no. 32, pp. 22467–22477, 2016.
  * (35) J. Rezác, K. E. Riley, and P. Hobza, “S66: A well-balanced database of benchmark interaction energies relevant to biomolecular structures,” J. Chem. Theory Comput., vol. 7, no. 8, pp. 2427–2438, 2011.
  * (36) B. Brauer, M. K. Kesharwani, S. Kozuch, and J. M. Martin, “The S66x8 benchmark for noncovalent interactions revisited: Explicitly correlated ab initio methods and density functional theory,” Phys. Chem. Chem. Phys., vol. 18, no. 31, pp. 20905–20925, 2016.
  * (37) A. Chantzis and G. Maroulis, “Interaction-induced electric properties in Kr–Ne from ab initio and DFT calculations. Is there a discrepancy between theory and experiment for the dipole moment?,” Chem. Phys. Lett., vol. 507, no. 1-3, pp. 42–47, 2011.
  * (38) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, “A fifth-order perturbation comparison of electron correlation theories,” Chem. Phys. Lett., vol. 157, no. 6, pp. 479–483, 1989.
  * (39) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al., “MOLPRO, version 2019.2, a package of ab initio programs,” 2019.  see https://www.molpro.net.
  * (40) G. Knizia, T. B. Adler, and H.-J. Werner, “Simplified CCSD(T)-F12 methods: Theory and benchmarks,” J. Chem. Phys., vol. 130, no. 5, p. 054104, 2009\.
  * (41) T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,” J. Chem. Phys., vol. 90, no. 2, pp. 1007–1023, 1989.
  * (42) K. A. Peterson, T. B. Adler, and H.-J. Werner, “Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, B–Ne, and Al–Ar,” J. Chem. Phys., vol. 128, no. 8, p. 084102, 2008.
  * (43) N. Sylvetsky, M. K. Kesharwani, and J. M. L. Martin, “The aug-cc-pVnZ-F12 basis set family: Correlation consistent basis sets for explicitly correlated benchmark calculations on anions and noncovalent complexes,” J. Chem. Phys., vol. 147, no. 13, p. 134106, 2017.
  * (44) K. E. Yousaf and K. A. Peterson, “Optimized auxiliary basis sets for explicitly correlated methods,” J. Chem. Phys., vol. 129, no. 18, p. 184108, 2008.
  * (45) K. E. Yousaf and K. A. Peterson, “Optimized complementary auxiliary basis sets for explicitly correlated methods: aug-cc-pVnZ orbital basis sets,” Chem. Phys. Lett., vol. 476, no. 4, pp. 303–307, 2009.
  * (46) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett., vol. 77, no. 18, p. 3865, 1996.
  * (47) D. Hait and M. Head-Gordon, “How accurate is density functional theory at predicting dipole moments? an assessment using a new database of 200 benchmark values,” J. Chem. Theory Comput., vol. 14, no. 4, pp. 1969–1981, 2018.
  * (48) D. Hait and M. Head-Gordon, “How accurate are static polarizability predictions from density functional theory? An assessment over 132 species at equilibrium geometry,” Phys. Chem. Chem. Phys., vol. 20, no. 30, pp. 19800–19810, 2018.
  * (49) 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., vol. 132, no. 15, p. 154104, 2010.
  * (50) S. Grimme, S. Ehrlich, and L. Goerigk, “Effect of the damping function in dispersion corrected density functional theory,” J. Comput. Chem., vol. 32, no. 7, pp. 1456–1465, 2011.
  * (51) O. A. Vydrov and T. Van Voorhis, “Nonlocal van der Waals density functional: The simpler the better,” J. Chem. Phys., vol. 133, no. 24, p. 244103, 2010\.
  * (52) C. Adamo and V. Barone, “Toward reliable density functional methods without adjustable parameters: The PBE0 model,” J. Chem. Phys., vol. 110, no. 13, pp. 6158–6170, 1999.
  * (53) J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hybrid functionals based on a screened Coulomb potential,” J. Chem. Phys., vol. 118, no. 18, pp. 8207–8215, 2003.
  * (54) J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)],” J. Chem. Phys., vol. 124, no. 21, p. 219906, 2006.
  * (55) O. A. Vydrov and G. E. Scuseria, “Assessment of a long-range corrected hybrid functional,” J. Chem. Phys., vol. 125, no. 23, p. 234109, 2006.
  * (56) J.-D. Chai and S.-P. Mao, “Seeking for reliable double-hybrid density functionals without fitting parameters: the PBE0-2 functional,” Chem. Phys. Lett., vol. 538, pp. 121–125, 2012.
  * (57) P. D. Mezei, G. I. Csonka, A. Ruzsinszky, and M. Kállay, “Construction and application of a new dual-hybrid random phase approximation,” J. Chem. Theory Comput., vol. 11, no. 10, pp. 4615–4626, 2015.
  * (58) G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, “The ground state correlation energy of the random phase approximation from a ring coupled cluster doubles approach,” J. Chem. Phys., vol. 129, no. 23, p. 231101, 2008.
  * (59) G. Jansen, R.-F. Liu, and J. G. Ángyán, “On the equivalence of ring-coupled cluster and adiabatic connection fluctuation-dissipation theorem random phase approximation correlation energy expressions,” J. Chem. Phys., vol. 133, no. 15, p. 154106, 2010.
  * (60) A. Heßelmann, “Random-phase-approximation correlation method including exchange interactions,” Phys. Rev. A, vol. 85, no. 1, p. 012517, 2012\.
  * (61) M. Kállay, Z. Rolik, J. Csontos, P. Nagy, G. Samu, D. Mester, J. Csóka, B. Szabó, I. Ladjánszki, L. Szegedy, B. Ladóczki, K. Petrov, M. Farkas, P. D. Mezei, and B. Hégely, “Mrcc, a quantum chemical program suite.”  See also Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, J. Chem. Phys. 139, 094105 (2013), as well as: www.mrcc.hu.
  * (62) B. Jeziorski, R. Moszynski, and K. Szalewicz, “Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes,” Chem. Rev., vol. 94, no. 7, pp. 1887–1930, 1994.
  * (63) E. G. Hohenstein, R. M. Parrish, C. D. Sherrill, J. M. Turney, and H. F. Schaefer III, “Large-scale symmetry-adapted perturbation theory computations via density fitting and Laplace transformation techniques: Investigating the fundamental forces of DNA-intercalator interactions,” J. Chem. Phys., vol. 135, no. 17, p. 174107, 2011.
  * (64) E. G. Hohenstein and C. D. Sherrill, “Density fitting and Cholesky decomposition approximations in symmetry-adapted perturbation theory: Implementation and application to probe the nature of $\pi$-$\pi$ interactions in linear acenes,” J. Chem. Phys., vol. 132, no. 18, p. 184111, 2010.
  * (65) E. Papajak, H. R. Leverentz, J. Zheng, and D. G. Truhlar, “Efficient Diffuse Basis Sets: cc-pV x Z+ and maug-cc-pV x Z,” J. Chem. Theory Comput., vol. 5, no. 5, pp. 1197–1202, 2009.
  * (66) J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, et al., “Psi4: an open-source ab initio electronic structure program,” Wiley Interdiscip. Rev. Comput. Mol. Sci., vol. 2, no. 4, pp. 556–565, 2012.
  * (67) T. G. A. Heijmen, R. Moszynski, P. E. S. Wormer, and A. V. D. Avoird, “Symmetry-adapted perturbation theory applied to interaction-induced properties of collisional complexes,” Mol. Phys., vol. 89, no. 1, pp. 81–110, 1996.
  * (68) Schrödinger, LLC, “The PyMOL molecular graphics system, version 1.8.” 2015\.
  * (69) Ż. Czyżnikowska, R. W. Góra, R. Zaleśny, W. Bartkowiak, A. Baranowska-Ła̧zkowska, and J. Leszczynski, “The effect of intermolecular interactions on the electric dipole polarizabilities of nucleic acid base complexes,” Chem. Phys. Lett., vol. 555, pp. 230–234, 2013.
  * (70) M. Medved’, A. Iglesias-Reguant, H. Reis, R. W. Góra, J. M. Luis, and R. Zaleśny, “Partitioning of interaction-induced nonlinear optical properties of molecular complexes. II. Halogen-bonded systems,” Phys. Chem. Chem. Phys., vol. 22, pp. 4225–4234, 2020.
  * (71) R. Zaleśny, M. Medved’, R. W. Góra, H. Reis, and J. M. Luis, “Partitioning of interaction-induced nonlinear optical properties of molecular complexes. I. Hydrogen-bonded systems,” Phys. Chem. Chem. Phys., vol. 20, pp. 19841–19849, 2018.
  * (72) R. W. Góra and B. Błasiak, “On the Origins of Large Interaction-Induced First Hyperpolarizabilities in Hydrogen-Bonded $\pi$-Electronic Complexes,” J. Phys. Chem. A, vol. 117, no. 31, pp. 6859–6866, 2013.
  * (73) W. Koch and M. C. Holthausen, A Chemist’s Guide to Density Functional Theory.  Wiley-VCH, 2002.
  * (74) S. Van Gisbergen, V. Osinga, O. Gritsenko, R. Van Leeuwen, J. Snijders, and E. Baerends, “Improved density functional theory results for frequency-dependent polarizabilities, by the use of an exchange-correlation potential with correct asymptotic behavior,” The Journal of chemical physics, vol. 105, no. 8, pp. 3142–3151, 1996.
  * (75) B. Champagne, E. A. Perpète, S. J. Van Gisbergen, E.-J. Baerends, J. G. Snijders, C. Soubra-Ghaoui, K. A. Robins, and B. Kirtman, “Assessment of conventional density functional schemes for computing the polarizabilities and hyperpolarizabilities of conjugated oligomers: An ab initio investigation of polyacetylene chains,” The Journal of chemical physics, vol. 109, no. 23, pp. 10489–10498, 1998.
  * (76) S. Y. Willow, X. C. Zeng, S. S. Xantheas, K. S. Kim, and S. Hirata, “Why is mp2-water “cooler” and “denser” than dft-water?,” The journal of physical chemistry letters, vol. 7, no. 4, pp. 680–684, 2016.
  * (77) D. Hait and M. Head-Gordon, “How accurate are static polarizability predictions from density functional theory? an assessment over 132 species at equilibrium geometry,” Physical Chemistry Chemical Physics, vol. 20, no. 30, pp. 19800–19810, 2018.
  * (78) S. Akter, J. A. Vargas, K. Sharkas, J. E. Peralta, K. A. Jackson, T. Baruah, and R. R. Zope, “How well do self-interaction corrections repair the overestimation of static polarizabilities in density functional calculations?,” Physical Chemistry Chemical Physics, vol. 23, no. 34, pp. 18678–18685, 2021.
  * (79) A. D. Becke, “Density-functional exchange-energy approximation with correct asymptotic behavior,” Phys. Rev. A, vol. 38, no. 6, p. 3098, 1988.
  * (80) C. Lee, W. Yang, and R. G. Parr, “Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density,” Phys. Rev. B, vol. 37, no. 2, p. 785, 1988.
  * (81) H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, “A long-range correction scheme for generalized-gradient-approximation exchange functionals,” J. Chem. Phys., vol. 115, no. 8, pp. 3540–3544, 2001.
  * (82) D. Xenides, A. Hantzis, and G. Maroulis, “Comparison of high-level post-Hartree–Fock and DFT methods on the calculation of interaction-induced electric properties of Kr–He,” Chem. Phys., vol. 382, no. 1-3, pp. 80–87, 2011.
  * (83) A. Zawada, A. Kaczmarek-Kędziera, and W. Bartkowiak, “Assessment of DFT functionals for the calculation of interaction-induced electric properties of molecular complexes,” Chem. Phys. Lett., vol. 503, no. 1-3, pp. 39–44, 2011.
  * (84) A. Zawada and W. Bartkowiak, “Many-body interactions and the electric response of hydrogen-bonded molecular chains,” Comput. Theor. Chem., vol. 967, no. 1, pp. 120–128, 2011.