|  Linear Response Properties of Solvated Systems: A Computational Study†
---|---
|  Linda Goletto,a‡ Sara Gómez,b Josefine H. Andersen,c Henrik Koch,a,b∗ and
Tommaso Giovanninib∗
|  We present a computational study of static and dynamic linear
polarizabilities in solution. We use different theoretical approaches to
describe solvent effects, ranging from quantum mechanics/molecular mechanics
(QM/MM) to quantum embedding approaches. In particular, we consider non-
polarizable and polarizable QM/MM methods, the latter based on the fluctuating
charge (FQ) force field. In addition, we use a quantum embedding method
defined in the context of multilevel Hartree-Fock (MLHF), where the system is
divided into active and inactive regions, and combine it with a third layer
described by means of the FQ model. The multiscale approaches are then used as
reference wave functions for equation-of-motion coupled cluster (EOM-CC)
response properties, allowing for the account of electron correlation. The
developed models are applied to the calculation of linear response properties
of two organic moieties—namely, para-nitroaniline and benzonitrile—in non-
aqueous solvents—1,4-dioxane, acetonitrile, and tetrahydrofuran. The computed
polarizabilities are then discussed in terms of the physico-chemical solute-
solvent interactions described by each method (electrostatic, polarization and
Pauli repulsion), and finally compared with the available experimental
references.

##

††footnotetext: a Department of Chemistry, Norwegian University of Science and
Technology (NTNU), NO-7491 Trondheim, Norway††footnotetext: b Scuola Normale
Superiore, Piazza dei Cavalieri 7, IT-56126 Pisa, PI Italy; E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>c DTU Chemistry,
Technical University of Denmark, DK-2800 Kongens Lyngby,
Denmark††footnotetext: ‡ Present address: Scuola Normale Superiore, Piazza dei
Cavalieri 7, I-56126 Pisa, Italy††footnotetext: † Electronic Supplementary
Information (ESI) available: data related to Fig. 3-9. See DOI:
10.1039/cXCP00000x/

##  1 Introduction

The response of a molecular system to an external electric field plays a
fundamental role in a plethora of technological applications.1 In this
context, theoretical chemistry can help understand the underlying physics of
the different phenomena. Among them, linear response properties are the most
basic quantities to be investigated, being the physico-chemical foundation of
many different spectroscopic signals.2, 3, 4, 5, 6 Therefore, the theoretical
modeling can have a pivotal role in gaining insight into how molecules and
complex systems behave in the presence of electromagnetic radiation.2 Of
particular interest are molecular systems embedded in an external environment,
being it a solvent or a biological matrix.7, 8 In fact, in such cases, the
molecular properties of the chromophore, which is usually the target of the
study, can be drastically perturbed by the presence of the environment.9

To face this kind of problems, the most widespread approach is to resort to
the so-called focused models.10, 11 In the specific case of solutions, a high
level of theory is used for the region of interest—the solute—while the
solvent is treated at a less sophisticated level.12 Thus, within the focused
model formalism, the system is usually partitioned into two layers. In
particular, the solute is generally described at the quantum mechanical (QM)
level, while the treatment of the solvent can range from a continuum,11 to
molecular mechanics (MM) approaches,10, 13, 14 to a lower-level QM model.15,
16, 17, 18, 19, 20, 21 Each of these subcategories encompasses many different
methods, where different solute-solvent interactions—e.g. electrostatic
effects, polarization, Pauli repulsion, dispersion—are taken into account, and
different computational efforts are required. In order to correctly describe
strong and specific solute-solvent interactions, the atomistic nature of the
environment usually needs to be retained in the modeling.22 Among the
atomistic models, the most commonly used belong to the family of QM/MM
approaches.23 In their basic formulation—the so-called electrostatic
embedding—the environment electrostatically perturbs the QM density, but not
vice versa.23 To include mutual solute-solvent polarization effects, which
might have a huge influence on the QM properties and spectra, polarizable
QM/MM approaches are exploited.24, 25, 26, 27 In this way, both electrostatic
and polarization contributions are taken into account. However, non-
electrostatic interactions such as Pauli repulsion and dispersion are usually
neglected, although they play a crucial role in many complex systems.28, 29,
30 To recover a theoretically consistent picture of such interactions, which
are intrinsically of quantum nature, quantum embedding approaches can be
used.15, 16, 17, 18, 19, 20, 21, 31 As mentioned above, these models are based
on the description at the QM level—although less sophisticated than the one
used for the solute—of at least a part of the environment. This allows for the
treatment of solute-solvent Pauli repulsion, and in some cases of dispersion
interactions too.32, 33, 34 Due to the quantum description of a larger part of
the system, quantum embedding approaches are generally more computationally
demanding than QM/MM methods. This problem can be solved by three-layer
approaches, where the largest part of the environment, usually the farthest
from the solute, is described by means of classical force fields.35, 21, 36,
37, 38, 39 In this way, within a small portion of the system most interactions
are treated at the QM level, whereas long-range contributions are retained at
the classical level only, providing a physically consistent picture. In this
work, we present a computational investigation of linear response properties
of two organic systems, namely para-nitroaniline (PNA) and benzonitrile
(PhCN), dissolved in dioxane (DIO), acetonitrile (ACN), and tetrahydrofuran
(THF). To quantify the solvent effects on such properties, we present a
hierarchy of solvation approaches, ranging from common QM/MM methods to three-
layer quantum embedding models. As for QM/MM approaches, we consider both non-
polarizable and polarizable frameworks. The latter is based on the fluctuating
charge (FQ) force field,40, 22, 41 which has been recently parametrized for
the selected solvents.42 The three-layer quantum embedding, on the other hand,
is based on the multilevel Hartree-Fock (MLHF) method.16 Within MLHF, the
molecular orbitals (MOs) are partitioned into active and inactive by means of
Cholesky decomposition43, 44, 45, 46 coupled with projected atomic orbitals47,
48 (PAOs) for the virtual space. The computational advantage of such a method
lies in the fact that the active MOs are optimized in the field of the
inactive ones, which are kept frozen, but orthogonal to the active space.
Therefore, electrostatic and Pauli repulsion (and part of the polarization)
active-inactive interactions are automatically taken into account at the HF
level. To refine the picture provided by the basic formulation of MLHF, the
active and inactive orbitals can be localized in their pre-defined spatial
regions by means of an energy-based procedure that we have recently
presented.49 The model obtained is called MLHF-AB. If such a procedure is
applied to an HF optimized wave function, fully accounting for solute-solvent
interactions, the resulting MOs are denoted as fragment localized MOs
(FLMOs).50 To minimize its computational cost, MLHF(-AB) can be coupled to an
external MM layer (MLHF(-AB)/MM).39

To calculate the linear response properties, we use the aforementioned two-
and three-layer wave functions as the reference for a post-HF description of
the solute. In fact, electron correlation has been proven particularly
significant for the accurate modeling of both static and dynamic
(hyper)polarizabilities.51, 52, 53, 54, 55, 56 If the ground state is
dominated by a single-determinant wave function, the coupled cluster (CC)
hierarchy of methods arguably provides one of the most sophisticated
descriptions of electron correlation.57 For this reason, coupled cluster is
often considered the theoretical golden standard for the prediction of linear
response properties, although many other _ab-initio_ methods, ranging from
density functional theory (DFT) to Møller-Plesset (MP) perturbation theory,
have been routinely used for this purpose.58, 59, 60, 61 Note that, when
dealing with excited states and molecular properties, coupled cluster methods
typically follow one of two routes: response theory62, 63, 64 (CCRT) or the
equation-of-motion65, 66 (EOM-CC) formalism, which is exploited here. The two
frameworks result in identical excitation energies, but differ in the
molecular properties, although with a generally small discrepancy.1

The manuscript is organized as follows. In the next section, we detail the
theoretical approaches with a focus on the solvation modeling and the
calculation of linear response properties at the EOM-CC level of theory. Then,
the computational protocol followed in the numerical analysis is presented and
applied to the calculation of static and dynamic polarizabilities of PNA and
PhCN dissolved in DIO, ACN, and THF. A summary and future perspectives end the
manuscript.

##  2 Theory

This section outlines the theoretical basis of the solvation methods employed
to compute the polarizabilities at the EOM-CC2 and EOM-CCSD levels, which are
also briefly described. In particular, we briefly recall the theory of non-
polarizable QM/MM and polarizable QM/FQ, together with that of the three-layer
MLHF-AB/MM method.

###  2.1 Non-polarizable QM/MM and polarizable QM/FQ

As stated above, QM/MM methods rely on the partitioning of the total energy of
the system into a QM ($E_{\text{QM}}$) and an MM ($E_{\text{MM}}$)
contribution23

$E=E_{\text{QM}}+E_{\text{MM}}+E^{\text{int}}_{\text{QM/MM}}\ ,$ (1)

where $E^{\text{int}}_{\text{QM/MM}}$ is the QM/MM interaction energy. In
electrostatic embedding, $E^{\text{int}}_{\text{QM/MM}}$ is limited to the
purely electrostatic interaction, whereas in polarizable embedding, a
polarization contribution is also included. In the former, each MM atom is
endowed with a fixed charge. In the latter, if the FQ force field is
exploited, the charge assigned to each atom can vary as a response to the QM
potential. Since both the non-polarizable QM/MM and polarizable QM/FQ depend
only on charges, the QM/MM interaction energy reads22

$E^{\text{int}}_{\text{QM/MM}}=\sum_{i}q_{i}V_{i}(\mathbf{D})\ ,$ (2)

where $V_{i}(\mathbf{D})$ is the QM potential due to the QM part acting on the
i$th$ MM charge ($q_{i}$). While in electrostatic embedding such charges are
fixed, in QM/FQ their values are obtained by solving the following set of
linear equations22, 41

$\mathbf{M}\mathbf{q}_{\lambda}=-\mathbf{C}_{Q}-\mathbf{V}(\mathbf{D}),$ (3)

where $\mathbf{q}_{\lambda}$ collects the FQ charges $q$ and suitable
Lagrangian multipliers that ensure charge conservation and $\mathbf{M}$ is a
matrix containing charge-charge interactions and Lagrangian blocks. The right
hand side consists instead of $\mathbf{C}_{Q}$, which takes into account
atomic electronegativities and the total charge constraints, and
$\mathbf{V}(\mathbf{D})$—the QM potential. In both QM/MM approaches, the QM
Fock matrix $F_{\mu\nu}$ (in the atomic orbital basis $\\{\chi_{\mu}\\}$) is
modified by the inclusion of the QM/MM interaction22

$F_{\mu\nu}=h_{\mu\nu}+G_{\mu\nu}(\mathbf{D})+\sum_{i}q_{i}V_{i,\mu\nu},$ (4)

where $h_{\mu\nu}$ and $G_{\mu\nu}$ are the one- and two-electron matrix
elements. This additional term is fixed in non-polarizable QM/MM, whereas it
varies at each self consistent field (SCF) step in QM/FQ, because the charges
$q$ depend on the QM density.

###  2.2 Multilevel Hartree-Fock

In the MLHF model,16 the total density matrix ($\mathbf{D}$) of the system is
decomposed into an active ($\mathbf{D}^{A}$) and an inactive
($\mathbf{D}^{B}$) component. Under this assumption, the total energy of a
system described at the HF level can be written as

$\displaystyle E^{\mathrm{TOT}}=$
$\displaystyle\tr\mathbf{h}\mathbf{D}^{A}+\frac{1}{2}\tr\mathbf{D}^{A}\mathbf{G}(\mathbf{D}^{A})+\tr\mathbf{D}^{A}\mathbf{G}(\mathbf{D}^{B})$
(5) $\displaystyle+$
$\displaystyle\tr\mathbf{h}\mathbf{D}^{B}+\frac{1}{2}\tr\mathbf{D}^{B}\mathbf{G}(\mathbf{D}^{B})+h_{\text{nuc}},$

where $\mathbf{h}$ and $\mathbf{G}$ are the one- and two-electron matrices,
and $h_{\text{nuc}}$ is the nuclear repulsion energy. In MLHF, only
$\mathbf{D}^{A}$ is iteratively optimized, whereas $\mathbf{D}^{B}$ is kept
fixed during the SCF optimization. Thus, the last three terms in Eq. 5, i.e.
the inactive energy and the nuclear repulsion, are constant throughout the
procedure. Additionally, the minimization is performed in the MO space of the
active part only, thus intrinsically reducing the computational cost of a full
HF description. As a consequence, the MLHF Fock matrix elements take the
following form

$F_{\mu\nu}=h_{\mu\nu}+G_{\mu\nu}(\mathbf{D}^{A})+G_{\mu\nu}(\mathbf{D}^{B}),$
(6)

where $G_{\mu\nu}(\mathbf{D}^{B})$ describes the interaction between the
active and inactive parts, and is indeed a one-electron term in the Fock
matrix, because $\mathbf{D}^{B}$ is fixed. Within the MLHF framework the
electrostatic, Pauli repulsion, and part of the polarization contributions
between the active and inactive parts are described at the HF level.

To further reduce the computational cost associated with an MLHF description,
MLHF can be coupled with an additional MM layer, yielding the MLHF/MM method
introduced in Ref. 39. Between the MLHF and MM parts, the interaction is
described at the purely electrostatic level, as in Eq. 2. In this case,
$\mathbf{D}$ is obtained as the sum of the active and inactive density
matrices. The MM layer can be equivalently treated at the non-polarizable or
polarizable FQ level. In the latter case, in Eq. 3 $\mathbf{D}$ refers to the
total density matrix ($\mathbf{D}^{A}+\mathbf{D}^{B}$). Finally, the MLHF Fock
matrix is modified by the coupling with the external MM layer as

$F_{\mu\nu}=h_{\mu\nu}+G_{\mu\nu}(\mathbf{D}^{A})+G_{\mu\nu}(\mathbf{D}^{B})+\sum_{i}q_{i}V_{i,\mu\nu}.$
(7)

After converging the MLHF(/MM) wave function, the active and inactive MOs can
be localized in their specific spatial regions by using an energy-based
localization of the MOs. In the resulting MLHF-AB approach,49 the $\mathbf{h}$
contribution is separated in terms of the kinetic operator and the $A$, $B$,
and interaction ($AB$) electron-nuclei potentials. Hence, the total energy can
be rewritten as

$\displaystyle E^{\mathrm{TOT}}$
$\displaystyle=\underbrace{\tr\mathbf{h}^{A}\mathbf{D}^{A}+\frac{1}{2}\tr\mathbf{D}^{A}\mathbf{G}(\mathbf{D}^{A})+h^{A}_{nuc}}_{E_{A}}$
(8)
$\displaystyle+\underbrace{\tr\mathbf{h}^{B}\mathbf{D}^{B}+\frac{1}{2}\tr\mathbf{D}^{B}\mathbf{G}(\mathbf{D}^{B})+h^{B}_{nuc}}_{E_{B}}$
$\displaystyle+\underbrace{\tr\mathbf{V}^{B}\mathbf{D}^{A}+\tr\mathbf{V}^{A}\mathbf{D}^{B}+\tr\mathbf{D}^{A}\mathbf{G}(\mathbf{D}^{B})+h^{AB}_{nuc}}_{E_{AB}}.$

The MOs of the A and B fragments are localized by means of a minimization of
the $E_{A}+E_{B}$ energy in the space spanned by the occupied MOs of the two
fragments. This means that the total density is not changed and that such a
minimization is equivalent to a maximization of the repulsion energy $E_{AB}$.
In this way, the occupied MOs of both fragments are localized in their
specific spatial regions. If a full HF optimization is performed before the
localization procedure, FLMOs are obtained,50 yielding the HFFLMOs approach.
Note that, if an additional MM layer is included in the modeling, it does not
affect the minimization procedure: the total density matrix remains fixed, and
so does the MLHF-MM interaction. The localization procedure outlined here
makes MLHF-AB a promising tool for the calculation of local properties, such
as dipole moments49 and polarizabilities.

###  2.3 EOM-CC2 and -CCSD linear polarizabilities

After the SCF convergence is reached with any of the aforementioned approaches
(QM/MM, MLHF-AB/MM, and HF${}_{\text{FLMOs}}$/MM), the polarizabilities are
computed at the EOM-CC2 and -CCSD level of theory by restraining the coupled
cluster treatment to the QM (in QM/MM) or to the active part (in MLHF-AB/MM
and HF${}_{\text{FLMOs}}$/MM) only.

The coupled cluster wave function is expressed as the exponential
parametrization57

$\ket{\mathrm{CC}}=e^{T}\ket{\mathrm{HF}},$ (9)

where $\ket{\mathrm{HF}}$ is the reference Hartree-Fock wave function and $T$
is the cluster operator—i.e., the sum of all the excitation operators
$\tau_{\nu}$ weighted by their amplitudes $t_{\nu}$

$\displaystyle T=$ $\displaystyle T_{1}+T_{2}+...$ (10) $\displaystyle=$
$\displaystyle\sum_{\nu_{1}}t_{\nu_{1}}\tau_{\nu_{1}}+\sum_{\nu_{2}}t_{\nu_{2}}\tau_{\nu_{2}}+...$

Here, $\nu_{n}$ refers to the $n$-th electronic excitation. In the CC267 and
CCSD68 models, the cluster operator $T$ is truncated after double excitations.
The difference between the two models lies in the fact that in CC2 the double
excitation component, $T_{2}$, is treated perturbatively. In CCSD, the doubles
amplitudes are considered to infinite order, while the CC2 doubles are
considered through first order only. The doubles amplitude equations of the
two models (for a $T_{1}$-transformed Hamiltonian,
$\bar{H}=e^{-T_{1}}He^{T_{1}}$) take the form67

$\displaystyle\mathrm{CCSD:}\quad$
$\displaystyle\bra{\nu_{2}}\bar{H}+[\bar{H},T_{2}]+\frac{1}{2}[[\bar{H},T_{2}],T_{2}]\ket{\mathrm{HF}}=0$
(11) $\displaystyle\mathrm{CC2:}\quad$
$\displaystyle\bra{\nu_{2}}\bar{H}+[F,T_{2}]\ket{\mathrm{HF}}=0$ (12)

while the singles equations remain the same for CC2 as in CCSD:

$\bra{\nu_{1}}\bar{H}+[\bar{H},T_{2}]\ket{\mathrm{HF}}=0$ (13)

where

$\ket{\nu}={\tau}_{\nu}\ket{\mathrm{HF}},\hskip
10.00002pt\bra{\nu}=\bra{\mathrm{HF}}{\tau}^{\dagger}_{\nu}.$ (14)

As a result, CC2 scales as $N^{5}$ compared to CCSD which is a $N^{6}$ model.

The CC2 model has a structure that is compatible with the exact linear
response functions, and the CC2 response properties of a molecule can thus be
computed within the same framework and formalism as CCSD.67 Molecular response
properties arise from the perturbation of an unperturbed system by an external
time-periodic field, and can be expressed in terms of response functions. The
second-order (linear) response to an external electric field gives rise to the
frequency-dependent electronic polarizability.

In the EOM-CC formalism,69, 65 the ground and excited states are explicitly
parametrized as

$\ket{k}=e^{{T}}\sum_{\nu\geq 0}R^{k}_{\nu}\ket{\nu},\hskip
10.00002pt\bra{\bar{k}}=\sum_{\nu\geq 0}\bra{\nu}L^{k}_{\nu}e^{-{T}}.$ (15)

Inserting the EOM states of eq.15 in the exact-state linear response
function,70 and considering the response of the components of the dipole
moment operator $\mu$, the EOM-CC linear electronic polarizability, $\alpha$,
reads66, 64

$\alpha_{XY}(\omega)=\frac{1}{2}P^{XY}\left(\bm{\eta}^{X}\bm{t}^{Y}(\omega)+\bm{\eta}^{X}\bm{t}^{Y}(-\omega)\right),$
(16)

where the permutation operator $P^{XY}$ performs an interchange of the
operators $X$ and $Y$, $\omega$ is the frequency of the external field, $X$
and $Y$ are cartesian components of $\mu$ and, for a generic operator
$\cal{O}$,71, 72

$\displaystyle\eta^{\cal{O}}_{\nu}$
$\displaystyle=\left(\bra{\mathrm{HF}}+\bra{\bar{t}}\right)\bar{\cal{O}}\tau_{\nu}\ket{\mathrm{HF}}-\bar{t}_{\nu}\matrixelement{\mathrm{HF}}{{\bar{\cal{O}}}}{\mathrm{HF}}$
(17)
$\displaystyle=\left(\bra{\mathrm{HF}}+\bra{\bar{t}}\right)[\bar{\cal{O}},\tau_{\nu}]\ket{\mathrm{HF}}+\sum_{\mu>\nu}\bar{t}_{\mu}\matrixelement{\mu}{\tau_{\nu}\bar{\cal{O}}}{\mathrm{HF}}-\bar{t}_{\nu}(\bar{t}\,\xi^{\cal{O}}),$
(18)

where we have introduced the similarity transformed operator
$\bar{\cal{O}}=e^{-T}{\cal{O}}e^{T}$, and
$\bra{\bar{t}}=\sum_{\mu}\bar{t}\bra{\mu}$ with $\bm{\bar{t}}$ being the
ground state Lagrangian multipliers. The amplitude response vectors
$\bm{t}^{\cal{O}}(\omega)$ are obtained by solving the linear equations

$(\mathbf{A}-\omega\mathbf{I})\bm{t}^{\cal{O}}(\omega)=-\bm{\xi}^{\cal{O}},$
(19)

where

$\xi^{\cal{O}}_{\nu}=\bra{\nu}\bar{\cal{O}}\ket{\mathrm{HF}},$ (20)

and $\mathbf{A}$ is the coupled cluster Jacobian matrix with elements

$A_{\nu\mu}=\bra{\nu}[\bar{H},\hat{\tau}_{\mu}]\ket{\mathrm{HF}}.$ (21)

##  3 Computational details

Fig. 1: Molecular structures of the solutes (para-nitroaniline, PNA, and
benzonitrile, PhCN) and solvents (1,4-dioxane, DIO, acetonitrile, ACN, and
tetrahydrofuran, THF).

In this work, we select different organic molecules dissolved in non-aqueous
environments, namely para-nitroaniline dissolved in 1,4-dioxane (PNA-in-DIO)
and benzonitrile dissolved in both acetonitrile (PhCN-in-ACN) and
tetrahydrofuran (PhCN-in-THF)—see Fig. 1. Such systems have been selected
because their measured linear polarizabilities have been previously reported
in the literature.73, 74

In order to correctly sample the solute-solvent phase-space, classical
molecular dynamics (MD) simulations are performed for both PNA and PhCN
dissolved in the different environments. In the case of PNA-in-DIO, the MD
simulation has been performed by following the procedure recently proposed in
Ref. 42. Similarly, for both PhCN-in-ACN and PhCN-in-THF, the General Amber
Force Field (GAFF)75 is used to describe the solute and solvents, for which
charges and parameters are obtained by using the RESP charge-fitting method76
and the Antechamber package,77 respectively. Optimized CAM-B3LYP/aug-cc-pVDZ
geometries are used to generate the force field parameters with the initial
solvent effects incorporated by means of the polarizable continuum model.11
PhCN is kept frozen during the MD runs, similarly to PNA (see Ref. 42). This
choice is justified by their planar and rigid structure, and also avoids any
potentially poor description of dihedral distributions by the classical force
field.78, 79 All simulations are performed using the GROMACS package.80
Following a similar methodology as in Ref. 42, our systems consist of a single
molecule of PhCN surrounded by thousands of solvent molecules and enclosed in
a simulation box of $7\text{\,}\mathrm{nm}$ size. They are minimized for 500
steps, prior to a $2\text{\,}\mathrm{ns}$ equilibration in the isothermal-
isobaric ensemble, keeping the temperature ($300\text{\,}\mathrm{K}$) and the
pressure ($1\text{\,}\mathrm{a}\mathrm{t}\mathrm{m}$) constant by means of a
velocity-rescaling method,81 with a coupling constant of
$0.1\text{\,}\mathrm{ps}$, and the Berendsen barostat,82 with a coupling
constant of $2.0\text{\,}\mathrm{ps}$, respectively. Values of
$9.7\text{\times}{10}^{-5}\text{\,}$ and
$9.6\text{\times}{10}^{-5}\text{\,}{\mathrm{bar}}^{-1}$ are used for the
isothermal compressibilities of THF and ACN, respectively. Afterward, an NVT
production stage of $10\text{\,}\mathrm{ns}$ is performed in order to have a
well-equilibrated system before extracting representative configurations. The
LINCS algorithm83 is used to constrain all bonds of the solute molecule. The
particle-mesh Ewald (PME) algorithm84 is employed to handle long-range
electrostatic interactions. Van der Waals and short-range electrostatic
interactions are truncated with a smoothed $1.4\text{\,}\mathrm{nm}$ spherical
cutoff. The equations of motion are integrated with a $2\text{\,}\mathrm{fs}$
time step.

A set of 20 snapshots is selected from the production stage of each MD
simulation. The time separation between them ($100\text{\,}\mathrm{ps}$ for
PNA-in-DIO and $500\text{\,}\mathrm{ps}$ for PhCN-in-ACN and PhCN-in-THF) is
large enough to ensure that they are uncorrelated.85, 86, 87, 88 A droplet
with a spherical shape of radius $20\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
centered on the solute is cut. Note that the number of selected frames is
enough to guarantee the convergence of the results (see Fig. S1-S3 and S14-S15
in the Electronic Supplementary Information – ESI†). The geometries of all the
frames studied in this work can be found in Ref.89. For each extracted
snapshot, the linear polarizability is then calculated by describing the whole
system at different levels of theory, defined within a hierarchical ladder:
(i) the solute is described at the QM level, whereas the environment is
described by means of electrostatic (QM/EE) or polarizable embedding (by
exploiting the FQ force field – QM/FQ); (ii) The solute and the closest
solvent molecules are included in the MLHF-AB region, while the remaining
solvent molecules are described at the FQ level. The solvent molecules within
a range of $2.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (PNA-in-DIO),
$2.75\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (PhCN-in-THF) and
$3.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (PhCN-in-ACN) from any atom of
the solute are included in the inactive MLHF-AB calculations, whereas the
solute molecule represents the active part. Such an approach is called MLHF-
AB/FQ in what follows. The same solvent molecules and the solute represent the
two regions described at the HF${}_{\text{FLMOs}}$ level in the
HF${}_{\text{FLMOs}}$/FQ approach. The partitioning of the spherical snapshots
at the different levels of theory is graphically depicted in Fig. 2, by taking
PNA-in-DIO as a representative example.

For MLHF-AB/FQ calculations, the protocol outlined in Ref. 39 is followed. A
superposition of molecular densities90 is used as a starting guess. While the
active MO virtual space is constructed at the beginning of the calculation by
means of orthonormalized91 PAOs,47, 48 the active occupied space is firstly
determined by a limited Cholesky decompositon algorithm,43, 44, 45, 46 and
then iteratively adjusted by maximizing the interaction energy $E_{AB}$ in Eq.
8. The number of active MOs is selected to be equal to the correct number of
occupied MOs of the active region.

In order to calculate the linear polarizability of each snapshot, the solute
is described at the EOM-CC2 or EOM-CCSD level with the aug-cc-pVDZ basis set,
by using the HF/EE, HF/FQ, MLHF-AB/FQ, and HF${}_{\text{FLMOs}}$/FQ reference
wave functions. The basis set is selected by following Refs. 92, 93, 94, 95.
For the MLHF-AB/FQ and HF${}_{\text{FLMOs}}$/FQ reference wave functions, the
solvent molecules are described with the cc-pVDZ basis set. In CC/EE, GAFF
atomic charges are used for the EE region.75, 77 For the FQ layer in CC/FQ,
CC-in-MLHF-AB/FQ and CC-in-HF${}_{\text{FLMOs}}$/FQ, the atomic
electronegativity and chemical hardness parameters have been taken from Ref.
42

Fig. 2: Graphical depiction of a snapshot of PNA-in-DIO as partitioned in the
CC/FQ and CC/EE (top) and CC-in-MLHF-AB/FQ and CC-in-HFFLMOs/FQ (bottom)
approaches. The atoms in blue are treated at the CC level, by using HF as
reference wavefunction; the atoms in orange are treated either at the HF
level—in HFFLMOs—or as the inactive MLHF part—in MLHF-AB; the atoms in grey
are treated at the MM level.

For each snapshot, we calculate the static and dynamic isotropic electronic
part of the polarizabilities, which is given by

$\alpha^{\text{iso}}=\frac{1}{3}(\alpha_{xx}+\alpha_{yy}+\alpha_{zz})$ (22)

In the static case, a reorientation term $\alpha_{\mu}$ is added to the purely
electronic term to yield the total static polarizability,
$\alpha_{0}^{\mathrm{tot}}$:

$\displaystyle\alpha_{0}^{\mathrm{tot}}$
$\displaystyle=\alpha^{\text{iso}}_{0}+\alpha_{\mu}.$ (23)
$\displaystyle\alpha_{\mu}$ $\displaystyle=\frac{|\mu|^{2}}{3k_{B}T},$ (24)

where $\mu$ is the molecular dipole moment, $k_{B}$ is the Boltzmann constant
and $T$ is the temperature.96, 7 The final isotropic polarizability is
obtained by averaging the results computed for each snapshot. It is worth
noting that local field effects induced on the active part by the polarization
of solvent under external radiation are not considered in this work, although
they might affect computed linear response properties.97, 92, 98 All the
calculations are performed with a locally modified version of the electronic
structure program $e^{T}$.99 The $e^{T}$ default thresholds are used for the
optimization of the reference and coupled cluster ground state wave functions,
as well as for the dipole moments and EOM-CC polarizabilities. The threshold
for the Cholesky decomposition of the two-electron repulsion integrals is set
to $10^{-4}$.

##  4 Numerical Results

All the aforementioned methods are used to calculate the linear
polarizabilities of PNA-in-DIO, PhCN-in-ACN and PhCN-in-THF. In this section,
the computed results are analyzed in terms of the different physico-chemical
solute-solvent interactions introduced by the different methods. The accuracy
and robustness of the approaches are then tested by comparing the computed
data with the available experimental results.73, 74

###  4.1 PNA-in-DIO

Let us discuss the case of PNA-in-DIO. PNA has been the focus of a large
variety of theoretical100, 101, 102, 92 and experimental103, 104, 105
investigations. It is characterized by a push-pull electronic structure,
presenting an electron-acceptor and an electron-donor functional groups on the
opposite sides of a $\pi$-conjugated aromatic structure (see also Fig. 1).
Such a feature implies that its optical properties are strongly influenced by
solvent effects,7, 106, 107, 108, 109, 110 making PNA a perfect candidate for
studying the performances of the different theoretical approaches.

Fig. 3: PNA-in-DIO snapshot-to-snapshot differences between CC-in-HFFLMOs/FQ,
CC-in-MLHF-AB/FQ, CC/FQ and CC/EE ground state dipole moments results.

In order to highlight the different solute-solvent physico-chemical
interactions taken into account by the different investigated approaches, we
analyze the results computed for each snapshot (see Fig. 3, 4 and 5). In Fig.
3a, the differences between the dipole moments computed at the CC/FQ and CC/EE
levels are graphically depicted. It can be noticed that CC/FQ dipole moments
are larger than CC/EE ones, independently of the solute-solvent
configurations, i.e. the snapshots. This is due to the inclusion of
polarization effects, described by means of the FQ force field, which increase
the magnitude of solute-solvent interactions, and consequently the computed
dipole moments. Remarkably, such an increase highly varies as a function of
the snapshot, ranging from 1 to about 2.2 Debye, and yields an increase of
about $17$-$18\%$ on average.

Moving to the three layers approaches, Fig. 3b depicts the CC-in-
HF${}_{\text{FLMOs}}$/FQ and CC-in-MLHF-AB/FQ $\Delta\mu$ results as a
function of the snapshot. This case provides a different picture. In fact, a
variability both in magnitude and in sign is reported between the two
approaches. However, for most snapshots, the CC-in-HF${}_{\text{FLMOs}}$/FQ
approach predicts larger dipole moments as compared to CC-in-MLHF-AB/FQ. Such
a finding can be explained by considering that within the
HF${}_{\text{FLMOs}}$/FQ reference, all the solute-solvent interactions are
fully accounted for at the HF level. At the MLHF-AB/FQ level, on the contrary,
only part of the solute-solvent polarization effects are taken into account,
and the inactive MOs—those belonging to the solvent—are not fully optimized.
As a consequence, for most snapshots, the full account of polarization of CC-
in-HF${}_{\text{FLMOs}}$/FQ yields large dipole moments. The negative
deviations can be instead related to the optimization of the MOs of the
inactive part, which enhances Pauli repulsion effects.

Finally, we compare the results obtained by exploiting CC-in-
HF${}_{\text{FLMOs}}$/FQ and CC-in-MLHF-AB/FQ to the CC/FQ values, as reported
in Fig. 3c and d, respectively. While CC-in-MLHF-AB/FQ generally yields a
decrease of the dipole moment, a variability in intensity and sign is reported
for CC-in-HF${}_{\text{FLMOs}}$/FQ with respect to CC/FQ. The numerical
results can again be discussed in light of the physico-chemical interactions
included in the different approaches. In the CC/FQ approach, the solute-
solvent interactions are limited to electrostatics and polarization, whereas
in both the three-layer methods Pauli repulsion effects are also taken into
consideration. On the one hand, this explains the average decrease of the
dipole moment reported for the two quantum embedding methods. On the other
hand, the variability depicted in Fig. 3d can suggest that for some snapshots
the CC/FQ approach is not able to fully account for the solute-solvent
polarization effects.

Fig. 4: PNA-in-DIO snapshot-to-snapshot differences between CC-in-HFFLMOs/FQ,
CC-in-MLHF-AB/FQ, CC/FQ and CC/EE results for the electronic component of the
static polarizability. Fig. 5: PNA-in-DIO snapshot-to-snapshot differences
between CC-in-HFFLMOs/FQ, CC-in-MLHF-AB/FQ, CC/FQ and CC/EE electronic dynamic
polarizability results ($589\text{\,}\mathrm{n}\mathrm{m}$).

We now move to the differences in the static and dynamic polarizabilities,
which are graphically depicted as a function of the snapshot in Fig. 4 and 5,
respectively. The above discussion for the dipole moments is generally valid
also for these linear-response properties, but with some noticeable
exceptions. In fact, the differences between CC/FQ–CC/EE, and CC-in-
HF${}_{\text{FLMOs}}$/FQ–CC-in-MLHF-AB/FQ polarizabilities follow the same
trends reported for $\mu$. In particular, the inclusion of polarization in
CC/FQ yields an increase of both static and dynamic polarizabilities, with a
larger effect on the dynamic one. The full account of polarization effects in
HF${}_{\text{FLMOs}}$ provides an overall increase of the computed properties,
albeit with some negative values. However, the effect of using
HF${}_{\text{FLMOs}}$ in place of MLHF-AB is much smaller than the effect of
using FQ in place of EE (on average a $0.2$-$0.4\%$ increase vs. a $4$-$7\%$
increase). The differences between the quantum embedding models and QM/FQ (see
Fig. 4c-d and 5c-d) are negative for all the snapshots, showing that the
inclusion of Pauli repulsion effects provides a general confinement of the
active density. As a consequence, both the static and the dynamic
polarizability values decrease by $2$-$3\%$, on average.

Method | $\mu[\mathrm{D}]$ | $\alpha_{0}^{\mathrm{iso}}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$ | $\alpha^{\mu}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$ | $\alpha_{0}^{\mathrm{tot}}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$ | $\alpha_{0}^{\mathrm{exp}}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$
---|---|---|---|---|---
CC2 in vacuo | $6.7$ | $10.2$ | $221.7$ | $231.9$ |
CCSD in vacuo | $6.9$ | $9.6$ | $230.2$ | $239.9$ |
CC2/EE | $9.1$ | $10.7$ | $403.8$ | $414.5$ | $404\pm 6$
CC2/FQ | $10.8$ | $11.1$ | $565.5$ | $576.6$
CC2-in-MLHF-AB/FQ | $10.5$ | $10.8$ | $536.4$ | $547.2$
CC2-in-HFFLMOs/FQ | $10.8$ | $10.8$ | $573.9$ | $584.8$
CCSD/EE | $9.1$ | $10.0$ | $406.1$ | $416.1$
CCSD/FQ | $10.7$ | $10.3$ | $558.6$ | $568.9$
CCSD-in-MLHF-AB/FQ | $10.4$ | $10.1$ | $532.4$ | $542.5$
CCSD-in-HFFLMOs/FQ | $10.8$ | $10.1$ | $569.4$ | $579.5$
Table 1: Calculated PNA-in-DIO dipole moments ($\mu$), isotropic electronic static polarizabilities ($\alpha_{0}^{\text{iso}}$), reorientation terms ($\alpha^{\mu}$), and total static polarizabilities ($\alpha_{0}^{\text{tot}}$). The experimental reference73 is also provided, along with CC2 and CCSD in vacuo values. Method | $\alpha_{\omega}^{\mathrm{iso}}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$ | $\alpha_{\omega}^{\mathrm{exp}}\big{[}\frac{\mathrm{cm}^{3}}{\mathrm{mol}}\big{]}$
---|---|---
CC2 in vacuo | $11.0$ |
CCSD in vacuo | $10.3$ |
CC2/EE | $11.9$ | $14.1\pm 0.4$
CC2/FQ | $12.6$
CC2-in-MLHF-AB/FQ | $12.2$
CC2-in-HFFLMOs/FQ | $12.3$
CCSD/EE | $10.9$
CCSD/FQ | $11.5$
CCSD-in-MLHF-AB/FQ | $11.2$
CCSD-in-HFFLMOs/FQ | $11.2$
Table 2: Calculated PNA-in-DIO dynamic polarizabilities
($589\text{\,}\mathrm{n}\mathrm{m}$). The experimental reference73 is also
provided, along with CC2 and CCSD in vacuo values.

Finally, let us move to comment on the averaged results, which can be compared
to the available experimental data. In Table 1 the averaged isotropic values
of both the electronic static polarizability ($\alpha_{0}^{\text{iso}}$) and
the dipole moment obtained with the different theoretical methods are
reported, together with the experimental value from Ref. 73. To better
quantify the solvent effects on the computed properties, Table 1 also lists
the vacuo values calculated at the CC2 and CCSD levels. The total static
polarizability, $\alpha_{0}^{\text{tot}}$, is computed using Eq. 23, so it
consists of an orientationally averaged polarizability ($\alpha_{\mu}$) and of
the purely electronic contribution $\alpha^{\text{iso}}_{0}$.

The CC2 and CCSD results mainly differ in two respects: while CC2
$\alpha^{\text{iso}}_{0}$ is larger than the CCSD counterpart, the opposite
holds for the dipole moments in the gas phase. Solvent effects can be
appreciated by comparing the gas-phase results with those obtained by
including a description of the embedding. In particular, an overall increase
of both the static polarizability ($4$-$8\%$) and the dipole moments
($33$-$61\%$) is highlighted. Such a trend is compatible with what has been
reported in the literature for similar systems.111, 108, 112 While the CC/EE
approaches provide the smallest values of $\mu$ ($\sim$
$1.5\text{\,}\mathrm{D}$ away from all other methods) and
$\alpha^{\text{iso}}_{0}$, QM/FQ reports the largest
$\alpha^{\text{iso}}_{0}$, and CC-in-HF${}_{\text{FLMOs}}$/FQ the largest
$\mu$. All the trends between the values computed with different solvation
approaches follow what has already been pointed out for Fig. 3–5. Indeed,
CC/FQ $\mu$ results are larger than CC/EE $\mu$, and a similar trend is
reported for CC-in-HFFLMOs/FQ as compared to CC-in-MLHF-AB/FQ. Such a result
is primarily due to the inclusion (in CC/FQ) and full accounting (in CC-in-
HFFLMOs/FQ) of polarization. When considering $\alpha^{\text{iso}}_{0}$, the
aforementioned trends remain valid. In addition, we note that the use of the
multilevel wave function as a reference decreases $\alpha^{\text{iso}}_{0}$
with respect to the CC/MM data. This is in agreement with Fig. 4 and 5, and is
due to the quantum confinement of the solute density as a result of the
solute-solvent Pauli repulsion contributions introduced by the multilevel
modeling. The differences between CC2 and CCSD dipole moment results are
small, regardless of the solvation model employed.

When comparing to the experimental reference in Ref. 73 (see Table 1), it is
worth remarking that the largest contribution to $\alpha^{tot}_{0}$ is given
by $\alpha^{\mu}\propto\mu^{2}$. Hence, the computed values of the total
static polarizability strongly depend on the numerical values of $\mu$.
Indeed, CC/EE reports the smallest $\alpha^{tot}_{0}$, whereas the largest
$\alpha^{tot}_{0}$ values are given by CC-in-HFFLMOs/FQ, which has a similar
performance as CC-in-MLHB-AB/FQ, and CC/FQ. The differences between CC/EE and
the other solvation methods range between $34$ and $45\%$. This is primarily
due to the differences in the computed dipole moments, as the
$\alpha^{\text{iso}}_{0}$ only slightly affects the final computed property.
Since CC/EE $\alpha^{tot}_{0}$ is closest to the experimental value, the most
important error source for the other methods lies in the overestimation of
$\mu$, which is enhanced when the polarization is included in the modeling,
and lowered by the Pauli repulsion effects introduced in multilevel methods. A
small difference in $\mu$ is reflected in large differences in $\alpha^{\mu}$,
and consequently in $\alpha_{0}^{\text{tot}}$. In fact, as commented above,
$\alpha_{0}^{\text{tot}}$ is calculated as the sum of
$\alpha_{0}^{\text{iso}}$ and $\alpha^{\mu}$, and our results show that
$\alpha^{\mu}$ is approximately 50 times larger than
$\alpha_{0}^{\text{iso}}$.

To remove the dependency of the results on the computed dipole moments, we
move to dynamic polarizabilities, for which $\alpha_{\mu}=0$. Therefore, the
computed $\alpha^{\text{iso}}_{\omega}$ can be directly compared to the
experimental values, which are given in Table 2. Similar trends as reported
for the static electronic polarizability (see Table 1) can be observed also
for the frequency-dependent case. The inclusion of solvent effects in the
modeling increases the computed values independently of the exploited method,
ranging between a $6\%$ (CCSD/EE) and a $15\%$ (CC2/FQ) shift. The gap between
CCSD and CC2 results is larger as compared to the static polarizability, with
the CCSD results being lower than CC2 by up to $9\%$ (CC/FQ). The trends
between the different approaches directly follow those discussed in the static
case: CC/FQ $\alpha^{\text{iso}}_{\omega}$ are larger than the corresponding
CC/EE values, and the same is generally valid between CC-in-
HF${}_{\text{FLMOs}}$/FQ and CC-in-MLHF-AB/FQ, due to the inclusion and the
full accounting of polarization effects, respectively. Similarly to the static
case, CC-in-HF${}_{\text{FLMOs}}$/FQ and CC-in-MLHF-AB provide very similar
results, indicating that the MLHF-AB/FQ method is able to account for most of
the polarization effects. On the contrary, by moving from CC/FQ to the
multilevel methods the value of the computed property decreases ($\sim
2$-$3\%$), because the Pauli repulsion yields a confinement effect in the
reference wave function.

Comparing the results in Table 2 with the experimental counterpart, we first
note that, differently from the static case, all the computed values are
smaller. The CC/EE results present the largest deviations, while the CC/FQ
ones are closest to experiment. As commented above, CC2 polarizabilities are
generally larger than CCSD ones, thus resulting in a better agreement with the
experimental reference. However, it is worth noting that CC2 might reportedly
overestimate linear polarizabilities,113, 114 probably due to an
overestimation of the dispersion interaction.115, 116, 117 Considering the
high level of theory employed in this work, the systematic underestimation of
all methods could be explained by the fact that our model discards the zero-
point correction and the pure vibrational contribution to the linear response
properties. While the latter plays a negligible role in determining dynamic
polarizabilities (due to the unfavourable dependence on the external
frequency), studies at the DFT level have shown that the former can increase
the purely electronic contribution by up to $6\%$.9 Considering that the
discrepancies with the experimental counterpart range between $1.5$ ($\sim
10\%$, with CC2/FQ) and $3.2$ ($\sim 22\%$, with CCSD/EE) cm3/mol, an overall
agreement with the experiment can be reported for almost all methods. Finally,
it is worth noting that our modeling neglects explicit terms arising from
polarization contributions in response equations, which may enhance the
computed linear polarizabilities, as reported in similar contexts.118, 119,
120, 121

###  4.2 PhCN-in-ACN and PhCN-in-THF

We now move to the case of PhCN-in-ACN and PhCN-in-THF polarizabilities, which
have been studied from the experimental point of view in Ref. 74 as a function
of the external frequency (ranging from $1.553\text{\,}$ to
$2.295\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$). We again model solvent
effects by means of CC/EE, CC/FQ and CC-in-MLHF-AB/FQ; the CC-in-HFFLMOs/FQ
methods have not been included in the comparison, as in the previous section
the differences with CC-in-MLHF-AB/FQ in the polarizabilities have been found
to be negligible, but the method is associated with a higher computational
cost.

Fig. 6: PhCN-in-ACN snapshot-to-snapshot differences between CC-in-MLHF-AB/FQ,
CC/FQ and CC/EE electronic dynamic polarizability results
($1.696\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$).

The snapshot-to-snapshot differences CC/FQ–CC/EE and CC-in-MLHF-AB/FQ–CC/FQ
are depicted in Fig. 6 and 7 for the PhCN-in-ACN and PhCN-in-THF systems,
respectively. The differences are presented for a specific frequency equal to
$1.696\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$, which corresponds to the
experimental sodium D line. The plots for all the other frequencies considered
in this work are reported in Fig. S4-S13 of the ESI†, and show similar trends.

For PhCN-in-ACN (see Fig. 6), the differences between CC/FQ and CC/EE are
negligible, ranging between $\pm$ 0.02 cm3/mol. Additionally, the differences
strongly depend on the specific solute-solvent configurations and display a
sign alternation, thus averaging out in the final property. Note that, in
contrast to what has been discussed in the previous section (see Fig. 5), the
sign alternation indicates that the parametrization of the EE modeling
overestimates electrostatics, similarly to other non-polarizable force
fields.122 The mutual polarization between PhCN and ACN appears to play a
minor role in the solute-solvent interaction, as compared to PNA-in-DIO
system. On the other hand, the introduction of the intermediate MLHF-AB layer
between the coupled cluster and FQ regions lowers the polarizabilities results
in all the snapshots, with differences ranging from 0.2 to 0.4 cm3/mol. This
indicates once again the confinement effects provided by the accounting for
solute-solvent Pauli repulsion.

Fig. 7: PhCN-in-THF snapshot-to-snapshot differences between CC-in-MLHF-AB/FQ,
CC/FQ and CC/EE electronic dynamic polarizability results
($1.696\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$).

For PhCN-in-THF, the CC/FQ values are smaller than the CC/EE values in all
snapshots. Such a finding is opposite to the PNA-in-DIO case, and shows that
the parametrization exploited in the EE force field is including
(overestimating) electrostatic effects. However, the CC/FQ–CC/EE differences
only reach the $1\%$ of the total value of the polarizability. Therefore the
THF polarization, albeit numerically more significant than that of ACN (see
Fig. 6), does not play a significant role for such systems. On the contrary,
the Pauli repulsion interactions introduced by the MLHF-AB layer have again a
larger influence on the polarizabilities, lowering the results by
approximately $3\%$ (and numerically by 0.1–0.25 cm3/mol).

Fig. 8: Calculated PhCN-in-ACN polarizabilities as a function of the external
frequency. Standard error bars at the 68% confidence interval are plotted.
PhCN in vacuo polarizabilities are also given, along with the experimental
benchmark from Ref. 74. $\text{Exp}^{a}$ employs the Lorentz local field
correction, while $\text{Exp}^{b}$ employs the Onsager local field correction.
Fig. 9: Calculated PhCN-in-THF polarizabilities as a function of the external
frequency. Standard error bars at the 68% confidence interval are plotted.
PhCN in vacuo polarizabilities are also given, along with the experimental
benchmark from Ref. 74. $\text{Exp}^{a}$ employs the Lorentz local field
correction, while $\text{Exp}^{b}$ employs the Onsager local field correction.

The computed averaged values of the PhCN-in-ACN static/dynamic
polarizabilities are plotted as a function of the external frequency in Fig. 8
(see Table S1 in the ESI${\dagger}$ for the numerical data), together with the
in vacuo data and the experimental results reproduced from Ref. 74. In
particular, we report two different experimental references, which are
obtained by applying the Lorentz ($\text{Exp}^{a}$) and Onsager
($\text{Exp}^{b}$) local field corrections to the measured refractive indexes.
Indeed, it is worth remarking that the reported data are not the measured
quantities—that is, the refractive indexes—, but the electronic part of the
polarizability extracted from them. Therefore, the reference data are
associated with an intrinsic systematic error related to the approach
exploited to extrapolate a microscopic quantity (the polarizability) from a
macroscopic one (the refractive index).

In Fig. 8, the averaged polarizabilities computed by all the different methods
are shown to follow the same trend with respect to the external frequency. In
particular, CC/EE and CC/FQ results are almost identical, whereas the
inclusion of the MLHF-AB layer lowers the polarizabilities by approximately
$3\%$. Indeed, with respect to the in vacuo results, CC/EE and CC/FQ yield an
increase in the polarizabilities, albeit with a negligible deviation ($<1\%$).
On the other hand, the values obtained at the CC-in-MLHF-AB/FQ level are $2$
to $3\%$ lower than the corresponding CC results in vacuo. Thus, depending on
how the environment is treated, solvent effects shift the polarizabilities to
opposite directions. As the main difference between the approaches lies in the
Pauli repulsion between PhCN and ACN being mostly taken into account in CC-in-
MLHF-AB/FQ, this contribution appears to have a significant influence on the
results. CCSD polarizabilities are approximately $4\%$ smaller than the
corresponding CC2 values with all the multiscale methods, following the same
trend observed in vacuo.

We now move to the comparison with the experimental reference.74 Regardless of
the choice of local field correction, being it Lorentz or Onsager, the
experimental results are lower than the computed values. In particular,
$\text{Exp}^{b}$ and $\text{Exp}^{a}$ polarizabilities differ of about $1\%$,
the former being the lowest. The best agreement with the experiment is
obtained by using CCSD-in-MLHF-AB/FQ, which presents $1$-$3\%$ and $2$-$4\%$
deviations from $\text{Exp}^{a}$ and $\text{Exp}^{b}$, respectively. In
particular, CCSD-in-MLHF-AB/FQ reports its largest discrepancy with the
experiment for the static polarizability. Here, the values in Ref. 74 are not
recovered from the experimental permittivity, but obtained by extrapolation
with a Cauchy-type dispersion curve fit, which might introduce further
inaccuracy. The worst agreement with the experiment is given by CC2/EE and
CC2/FQ, which deviate from the experiment by $9\%$ and $12\%$, respectively.
Remarkably, among the considered methods, only the three-layer CC-in-MLHF-
AB/FQ approaches shift the computed polarizabilities from the corresponding in
vacuo results towards the experiment.

Fig. 9 depicts the averaged polarizability values (see Table S2 in the
ESI${\dagger}$ for the numerical data), the in vacuo reference, and
experimental benchmark74 for the PhCN-in-THF system, plotted again as a
function of the external frequencies. Similarly to the previous case, CC/EE
and CC/FQ results are very similar, with CC/FQ providing computed
polarizabilities less than $1\%$ smaller than CC/EE. This is in agreement with
the results reported in Fig. 7. An additional decrease is given by CC-in-MLHF-
AB/FQ, with values $1$-$2\%$ smaller than CC/FQ. The CC2 and CCSD values show
the same trend with respect to a different treatment of the environment, with
CCSD providing polarizabilities $4$-$5\%$ smaller than CC2. The same behaviour
is observed for the in vacuo results. When comparing to the in vacuo
reference, the CC/FQ values are almost identical, with deviations $<0.3\%$.
While being very similar, the CC/EE polarizabilities are slightly larger, even
if the differences fall below $1\%$. On the other hand, the CC-in-MLHF-AB/FQ
methods lower the polarizabilities by $2\%$ with respect to the in vacuo
results.

Additionally, Fig. 9 shows that the experimental references are once again
smaller than the computed values. The comparison between the methods follows
the same trends as in PhCN-in-ACN, with CCSD-in-MLHF-AB/FQ being the closest
(with a $2$-$4\%$ deviation) and CC2/EE being the furthest (with a $9$-$11\%$
deviation) from the experiment. While starting from the in vacuo reference
CC/EE goes in the wrong direction with respect to the experiment, both CC/FQ
and CC-in-MLHF-AB/FQ lower the polarizabilities towards the experiment.

##  5 Summary and conclusions

We have presented a computational investigation of linear polarizabilities of
organic moieties embedded in non-aqueous solvents, employing different
strategies to model solvent effects. We have considered a hierarchy of
solvation approaches, ranging from common QM/MM methods to three-layer
approaches based on a multilevel partitioning of the reference wave function.
In particular, we have considered both non-polarizable and polarizable QM/MM
approaches, the latter based on the FQ force field, suitably parameterized for
the selected solvents. The three-layer approaches are instead based on a
partitioning of the system into three portions: an active region (the solute),
an inactive region (the solvent molecules closest to the solute), and an MM
region (the rest of the solvent), treated by means of the FQ force field. In
this way, the external MM layer accounts for long range interactions with an
electrostatic and polarization nature. The MLHF-AB reference wave function
introduces the electron repulsion effects between solute and solvent, as well
as an approximated HF treatment of the polarization.

The approaches have been applied to the calculation of static and dynamic
linear polarizabilities (at the CC2 and CCSD levels) of the PNA-in-DIO, PhCN-
in-ACN, and PhCN-in-THF systems. To sample the solute-solvent phase-space, the
calculations have been run on various snapshots extracted from classical MD
simulations. Overall, this protocol gives a consistent physical description of
the properties and interactions at work in a complex environment. The computed
results have been rationalized in terms of the different solute-solvent
physico-chemical interactions modeled by each solvation approach and compared
with the available experimental data. In all cases, we have obtained an
overall good agreement with the reference measurements, in particular when
Pauli repulsion effects, which are introduced in the three-layer approaches,
are taken into account.

To further increase the agreement with experimental results, the three-layer
model could be further improved to include dispersion effects between the
solute and the solvent. Also, the quality of long-range electrostatics and
polarization effects can be increased by including an additional source of
polarization in the FQ force field in terms of fluctuating dipoles (FQF$\mu$)
to account for anisotropic interactions.123, 124, 125 Finally, the protocol
can be extended to the treatment of higher-order properties, such as first-
hyperpolarizabilities ($\beta$).

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We thank Eirik F. Kjønstad, Sarai D. Folkestad and Alexander C. Paul for their
contributions to the $e^{T}$ code. J.H.A. acknowledges Sonia Coriani (DTU) for
discussions. We acknowledge funding from the Marie Sklodowska-Curie
Interational Training Network “COSINE - COmputational Spectroscopy In Natural
sciences and Engineering”, Grant Agreement No. 765739, and from the Research
Council of Norway through the FRINATEK project 275506, TheoLight. We
acknowledge computing resources through UNINETT Sigma2—the National
Infrastructure for High Performance Computing and Data Storage in Norway,
through project number NN2962k. We also acknowledge Chiara Cappelli (SNS) for
computing resources, and the Center for High Performance Computing (CHPC) at
SNS for providing the computational infrastructure.

## Notes and references

  * Helgaker _et al._ 2012 T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen and K. Ruud, _Chem. Rev._ , 2012, 112, 543–631.
  * Norman _et al._ 2018 P. Norman, K. Ruud and T. Saue, _Principles and practices of molecular properties: Theory, modeling, and simulations_ , John Wiley & Sons, 2018\.
  * Amos 1986 R. Amos, _Chem. Phys. Lett._ , 1986, 124, 376–381.
  * Jensen and Jensen 2009 L. L. Jensen and L. Jensen, _J. Phys. Chem. C_ , 2009, 113, 15182–15190.
  * Norman _et al._ 1997 P. Norman, Y. Luo, D. Jonsson and H. Ågren, _J. Chem. Phys._ , 1997, 106, 8788–8791.
  * Norman _et al._ 2003 P. Norman, A. Jiemchooroj and B. E. Sernelius, _J. Chem. Phys._ , 2003, 118, 9167–9174.
  * Cammi _et al._ 2000 R. Cammi, B. Mennucci and J. Tomasi, _J. Phys. Chem. A_ , 2000, 104, 4690–4698.
  * Pedersen _et al._ 2014 M. N. Pedersen, E. D. Hedegard, J. M. H. Olsen, J. Kauczor, P. Norman and J. Kongsted, _J. Chem. Theor. Comput._ , 2014, 10, 1164–1171.
  * Egidi _et al._ 2014 F. Egidi, T. Giovannini, M. Piccardo, J. Bloino, C. Cappelli and V. Barone, _J. Chem. Theory Comput._ , 2014, 10, 2456–2464.
  * Warshel and Levitt 1976 A. Warshel and M. Levitt, _J. Mol. Biol._ , 1976, 103, 227–249.
  * Tomasi _et al._ 2005 J. Tomasi, B. Mennucci and R. Cammi, _Chem. Rev._ , 2005, 105, 2999–3094.
  * Mennucci and Corni 2019 B. Mennucci and S. Corni, _Nat. Rev. Chem._ , 2019, 3, 315–330.
  * Lin and Truhlar 2007 H. Lin and D. G. Truhlar, _Theor. Chem. Acc._ , 2007, 117, 185–199.
  * Senn and Thiel 2009 H. M. Senn and W. Thiel, _Angew. Chem. Int. Ed._ , 2009, 48, 1198–1229.
  * Marrazzini _et al._ 2021 G. Marrazzini, T. Giovannini, M. Scavino, F. Egidi, C. Cappelli and H. Koch, _J. Chem. Theory Comput._ , 2021, 17, 791–803.
  * Sæther _et al._ 2017 S. Sæther, T. Kjærgaard, H. Koch and I.-M. Høyvik, _J. Chem. Theor. Comput._ , 2017, 13, 5282–5290.
  * Svensson _et al._ 1996 M. Svensson, S. Humbel, R. D. Froese, T. Matsubara, S. Sieber and K. Morokuma, _J. Phys. Chem._ , 1996, 100, 19357–19363.
  * Govind _et al._ 1999 N. Govind, Y. A. Wang and E. A. Carter, _J. Chem. Phys._ , 1999, 110, 7677–7688.
  * Wesolowski _et al._ 2015 T. A. Wesolowski, S. Shedge and X. Zhou, _Chem. Rev._ , 2015, 115, 5891–5928.
  * Sun and Chan 2016 Q. Sun and G. K.-L. Chan, _Acc. Chem. Res._ , 2016, 49, 2705–2712.
  * Olsen _et al._ 2015 J. M. H. Olsen, C. Steinmann, K. Ruud and J. Kongsted, _J. Phys. Chem. A_ , 2015, 119, 5344–5355.
  * Giovannini _et al._ 2020 T. Giovannini, F. Egidi and C. Cappelli, _Chem. Soc. Rev._ , 2020, 49, 5664–5677.
  * Senn and Thiel 2009 H. M. Senn and W. Thiel, _Angew. Chem. Int. Ed._ , 2009, 48, 1198–1229.
  * Curutchet _et al._ 2009 C. Curutchet, A. Muñoz-Losa, S. Monti, J. Kongsted, G. D. Scholes and B. Mennucci, _J. Chem. Theory Comput._ , 2009, 5, 1838–1848.
  * Bondanza _et al._ 2020 M. Bondanza, M. Nottoli, L. Cupellini, F. Lipparini and B. Mennucci, _Phys. Chem. Chem. Phys._ , 2020, 22, 14433–14448.
  * Olsen and Kongsted 2011 J. M. H. Olsen and J. Kongsted, _Adv. Quantum Chem._ , 2011, 61, 107–143.
  * Olsen _et al._ 2010 J. M. Olsen, K. Aidas and J. Kongsted, _J. Chem. Theory Comput._ , 2010, 6, 3721–3734.
  * Giovannini _et al._ 2017 T. Giovannini, P. Lafiosca and C. Cappelli, _J. Chem. Theory Comput._ , 2017, 13, 4854–4870.
  * Slipchenko 2016 L. V. Slipchenko, in _Many-Body Effects and Electrostatics in Biomolecules_ , Pan Stanford, 2016, pp. 147–187.
  * Giovannini _et al._ 2019 T. Giovannini, M. Ambrosetti and C. Cappelli, _J. Phys. Chem. Lett._ , 2019, 10, 5823–5829.
  * Reinholdt _et al._ 2017 P. Reinholdt, J. Kongsted and J. M. H. Olsen, _J. Phys. Chem. Lett._ , 2017, 8, 5949–5958.
  * Folkestad and Koch 2019 S. D. Folkestad and H. Koch, _J. Chem. Theory Comput._ , 2019, 16, 179–189.
  * Myhre _et al._ 2014 R. H. Myhre, A. M. Sánchez de Merás and H. Koch, _J. Chem. Phys._ , 2014, 141, 224105.
  * Myhre and Koch 2016 R. H. Myhre and H. Koch, _J. Chem. Phys._ , 2016, 145, 044111\.
  * Wanko _et al._ 2008 M. Wanko, M. Hoffmann, T. Frauenheim and M. Elstner, _J. Phys. Chem. B_ , 2008, 112, 11462–11467.
  * Bennie _et al._ 2016 S. J. Bennie, M. W. van der Kamp, R. C. Pennifold, M. Stella, F. R. Manby and A. J. Mulholland, _J. Chem. Theory Comput._ , 2016, 12, 2689–2697.
  * Nogueira _et al._ 2018 J. J. Nogueira, S. Roßbach, C. Ochsenfeld and L. Gonzalez, _J. Chem. Theor. Comput._ , 2018, 14, 4298–4308.
  * Macetti and Genoni 2021 G. Macetti and A. Genoni, _J. Phys. Chem. A_ , 2021, 125, 6013–6027.
  * Goletto _et al._ 2021 L. Goletto, T. Giovannini, S. D. Folkestad and H. Koch, _Phys. Chem. Chem. Phys._ , 2021, 23, 4413–4425.
  * Cappelli 2016 C. Cappelli, _Int. J. Quantum Chem._ , 2016, 116, 1532–1542.
  * Giovannini _et al._ 2020 T. Giovannini, F. Egidi and C. Cappelli, _Phys. Chem. Chem. Phys._ , 2020, 22, 22864–22879.
  * Ambrosetti _et al._ 2021 M. Ambrosetti, S. Skoko, T. Giovannini and C. Cappelli, _J. Chem. Theory Comput._ , 2021, 17, 7146–7156.
  * Beebe and Linderberg 1977 N. H. Beebe and J. Linderberg, _Int. J. Quantum Chem._ , 1977, 12, 683–705.
  * Sánchez de Merás _et al._ 2010 A. M. Sánchez de Merás, H. Koch, I. G. Cuesta and L. Boman, _J. Chem. Phys._ , 2010, 132, 204105.
  * Aquilante _et al._ 2006 F. Aquilante, T. Bondo Pedersen, A. Sánchez de Merás and H. Koch, _J. Chem. Phys._ , 2006, 125, 174101.
  * Folkestad _et al._ 2019 S. D. Folkestad, E. F. Kjønstad and H. Koch, _J. Chem. Phys._ , 2019, 150, 194112.
  * Pulay 1983 P. Pulay, _Chem. Phys. Lett._ , 1983, 100, 151–154.
  * Saebo and Pulay 1993 S. Saebo and P. Pulay, _Annu. Rev. Phys. Chem._ , 1993, 44, 213–236.
  * Giovannini and Koch 2020 T. Giovannini and H. Koch, _J. Chem. Theory Comput._ , 2020, 17, 139–150.
  * Giovannini and Koch 2022 T. Giovannini and H. Koch, _J. Chem. Theory Comput._ , 2022, 18, 4806–4813.
  * Champagne _et al._ 2005 B. Champagne, E. Botek, M. Nakano, T. Nitta and K. Yamaguchi, _J. Chem. Phys._ , 2005, 122, 114315.
  * Pecul 2005 M. Pecul, _Chem. Phys. Lett._ , 2005, 404, 217–226.
  * Christiansen and Mikkelsen 1999 O. Christiansen and K. V. Mikkelsen, _J. Chem. Phys._ , 1999, 110, 8348–8360.
  * Kongsted _et al._ 2002 J. Kongsted, A. Osted, K. V. Mikkelsen and O. Christiansen, _Mol. Phys._ , 2002, 100, 1813–1828.
  * Kongsted _et al._ 2003 J. Kongsted, A. Osted, K. V. Mikkelsen and O. Christiansen, _J. Chem. Phys._ , 2003, 118, 1620–1633.
  * Hrsak _et al._ 2018 D. Hrsak, J. M. H. Olsen and J. Kongsted, _J. Chem. Theory Comput._ , 2018, 14, 1351–1360.
  * Helgaker _et al._ 2014 T. Helgaker, P. Jørgensen and J. Olsen, _Molecular electronic-structure theory_ , John Wiley & Sons, 2014.
  * Limacher _et al._ 2009 P. A. Limacher, K. V. Mikkelsen and H. P. Lüthi, _J. Chem. Phys._ , 2009, 130, 194114.
  * Baranowska-Łaczkowska _et al._ 2013 A. Baranowska-Łaczkowska, W. Bartkowiak, R. W. Góra, F. Pawłowski and R. Zaleśny, _J. Comput. Chem._ , 2013, 34, 819–826.
  * Wormer and Rijks 1986 P. E. Wormer and W. Rijks, _Phys. Rev. A_ , 1986, 33, 2928.
  * Olsen and Jørgensen 1985 J. Olsen and P. Jørgensen, _J. Chem. Phys._ , 1985, 82, 3235–3264.
  * Monkhorst 1977 H. J. Monkhorst, _Int. J. Quantum Chem._ , 1977, 12, 421–432.
  * Koch and Jørgensen 1990 H. Koch and P. Jørgensen, _J. Chem. Phys._ , 1990, 93, 3333–3344.
  * Pedersen and Koch 1997 T. B. Pedersen and H. Koch, _J. Chem. Phys._ , 1997, 106, 8059–8072.
  * Stanton and Bartlett 1993 J. F. Stanton and R. J. Bartlett, _J. Chem. Phys._ , 1993, 98, 7029–7039.
  * Kobayashi _et al._ 1994 R. Kobayashi, H. Koch and P. Jørgensen, _Chem. Phys. Lett._ , 1994, 219, 30–35.
  * Christiansen _et al._ 1995 O. Christiansen, H. Koch and P. Jørgensen, _Chem. Phys. Lett._ , 1995, 243, 409–418.
  * Purvis III and Bartlett 1982 G. D. Purvis III and R. J. Bartlett, _J. Chem. Phys._ , 1982, 76, 1910–1918.
  * Stanton and Bartlett 1993 J. F. Stanton and R. J. Bartlett, _J. Chem. Phys._ , 1993, 99, 5178–5183.
  * Linderberg and Öhrn 2004 J. Linderberg and Y. Öhrn, _Propagators in quantum chemistry_ , John Wiley & Sons, 2004.
  * Pawłowski _et al._ 2015 F. Pawłowski, J. Olsen and P. Jørgensen, _J. Chem. Phys._ , 2015, 142, 114109.
  * Faber and Coriani 2019 R. Faber and S. Coriani, _J. Chem. Theory Comput._ , 2019, 15, 520–528.
  * Wortmann _et al._ 1993 R. Wortmann, P. Krämer, C. Glania, S. Lebus and N. Detzer, _Chem. Phys._ , 1993, 173, 99–108.
  * Alvarado _et al._ 2003 Y. Alvarado, P. Labarca, N. Cubillán and A. Karam, _Z. Naturforsch. A_ , 2003, 58, 68–74.
  * Wang _et al._ 2004 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, _J. Comput. Chem._ , 2004, 25, 1157–1174.
  * Bayly _et al._ 1993 C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, _J. Phys. Chem._ , 1993, 97, 10269–10280.
  * Wang _et al._ 2006 J. Wang, W. Wang, P. A. Kollman and D. A. Case, _J. Mol. Graph. Model._ , 2006, 25, 247–260.
  * Kjellgren _et al._ 2018 E. R. Kjellgren, J. M. Haugaard Olsen and J. Kongsted, _J. Chem. Theor. Comput._ , 2018, 14, 4309–4319.
  * Giovannini _et al._ 2018 T. Giovannini, G. Del Frate, P. Lafiosca and C. Cappelli, _Phys. Chem. Chem. Phys._ , 2018, 20, 9181–9197.
  * Pronk _et al._ 2013 S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, _Bioinformatics_ , 2013, 29, 845–854.
  * Bussi _et al._ 2007 G. Bussi, D. Donadio and M. Parrinello, _J. Chem. Phys._ , 2007, 126, 014101.
  * Berendsen and Van Gunsteren 1986 H. J. Berendsen and W. F. Van Gunsteren, _Molecular-dynamics simulation of statistical-mechanical systems_ , 1986, 43–65.
  * Hess _et al._ 1997 B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, _J. Comput. Chem._ , 1997, 18, 1463–1472.
  * Darden _et al._ 1993 T. Darden, D. York and L. Pedersen, _J. Chem. Phys._ , 1993, 98, 10089–10092.
  * Reinholdt _et al._ 2018 P. Reinholdt, M. S. Nørby and J. Kongsted, _J. Chem. Theor. Comput._ , 2018, 14, 6391–6404.
  * Harczuk _et al._ 2015 I. Harczuk, O. Vahtras and H. Ågren, _Phys. Chem. Chem. Phys._ , 2015, 17, 7800–7812.
  * Skoko _et al._ 2020 S. Skoko, M. Ambrosetti, T. Giovannini and C. Cappelli, _Molecules_ , 2020, 25, 5853.
  * Puglisi _et al._ 2019 A. Puglisi, T. Giovannini, L. Antonov and C. Cappelli, _Phys. Chem. Chem. Phys._ , 2019, 21, 15504–15514.
  * Goletto _et al._ 2022 L. Goletto, S. Gómez, J. H. Andersen, H. Koch and T. Giovannini, _Geometries for ‘Linear Response Properties of Solvated Systems: A Computational Study’_ , 2022, https://doi.org/10.5281/zenodo.7230203.
  * Neugebauer _et al._ 2005 J. Neugebauer, M. J. Louwerse, E. J. Baerends and T. A. Wesolowski, _J. Chem. Phys._ , 2005, 122, 094115.
  * Löwdin 1970 P.-O. Löwdin, in _Adv. Quantum Chem._ , Elsevier, 1970, vol. 5, pp. 185–199.
  * Egidi _et al._ 2014 F. Egidi, M. Segado, H. Koch, C. Cappelli and V. Barone, _J. Chem. Phys._ , 2014, 141, 224114.
  * Cuesta _et al._ 2004 I. G. Cuesta, T. B. Pedersen, H. Koch and A. S. de Merás, _Chem. Phys. Lett._ , 2004, 390, 170–175.
  * Alparone 2011 A. Alparone, _Comput. Theor. Chem._ , 2011, 976, 188–190.
  * Alparone 2013 A. Alparone, _Chem. Phys._ , 2013, 410, 90–98.
  * Liptay _et al._ 1982 W. Liptay, J. Becker, D. Wehning, W. Lang and O. Burkhard, _Z. Naturforsch. A_ , 1982, 37, 1396–1408.
  * List _et al._ 2016 N. H. List, H. J. A. Jensen and J. Kongsted, _Phys. Chem. Chem. Phys._ , 2016, 18, 10070–10080.
  * Tomasi _et al._ 2002 J. Tomasi, R. Cammi, B. Mennucci, C. Cappelli and S. Corni, _Phys. Chem. Chem. Phys._ , 2002, 4, 5697–5712.
  * Folkestad _et al._ 2020 S. D. Folkestad, E. F. Kjønstad, R. H. Myhre, J. H. Andersen, A. Balbi, S. Coriani, T. Giovannini, L. Goletto, T. S. Haugland, A. Hutcheson, I.-M. Høyvik, T. Moitra, A. C. Paul, M. Scavino, A. S. Skeidsvoll, Åsmund H. Tveten and H. Koch, _J. Chem. Phys._ , 2020, 152, 184103\.
  * Daniel and Dupuis 1990 C. Daniel and M. Dupuis, _Chem. Phys. Lett._ , 1990, 171, 209–216.
  * Karna _et al._ 1991 S. P. Karna, P. N. Prasad and M. Dupuis, _J. Chem. Phys._ , 1991, 94, 1171–1181.
  * Champagne 1996 B. Champagne, _Chem. Phys. Lett._ , 1996, 261, 57–65.
  * Millefiori _et al._ 1977 S. Millefiori, G. Favini, A. Millefiori and D. Grasso, _Spectrochim. Acta A_ , 1977, 33, 21–27.
  * Carsey _et al._ 1979 T. Carsey, G. Findley and S. McGlynn, _J. Am. Chem. Soc._ , 1979, 101, 4502–4510.
  * Woodford _et al._ 1997 J. Woodford, M. Pauley and C. Wang, _J. Phys. Chem. A_ , 1997, 101, 1989–1992.
  * Cammi _et al._ 1998 R. Cammi, B. Mennucci and J. Tomasi, _J. Am. Chem. Soc._ , 1998, 120, 8834–8847.
  * Painelli and Terenziani 2001 A. Painelli and F. Terenziani, _Synth. Met._ , 2001, 124, 171–173.
  * Cammi _et al._ 2003 R. Cammi, L. Frediani, B. Mennucci and K. Ruud, _J. Chem. Phys._ , 2003, 119, 5818–5827.
  * Giovannini _et al._ 2019 T. Giovannini, R. R. Riso, M. Ambrosetti, A. Puglisi and C. Cappelli, _J. Chem. Phys._ , 2019, 151, 174104.
  * Kosenkov and Slipchenko 2010 D. Kosenkov and L. V. Slipchenko, _J. Phys. Chem. A_ , 2010, 115, 392–401.
  * Mikkelsen _et al._ 1994 K. V. Mikkelsen, Y. Luo, H. Ågren and P. Jørgensen, _J. Chem. Phys._ , 1994, 100, 8240–8250.
  * Targema _et al._ 2013 M. Targema, N. O. Obi-Egbedi and M. D. Adeoye, _Comput. Theor. Chem._ , 2013, 1012, 47–53.
  * Christiansen _et al._ 1999 O. Christiansen, J. Gauss and J. F. Stanton, _Chem. Phys. Lett._ , 1999, 305, 147–155.
  * Salek _et al._ 2005 P. Salek, T. Helgaker, O. Vahtras, H. Ågren, D. Jonsson and J. Gauss, _Mol. Phys._ , 2005, 103, 439–450.
  * Rocha-Rinza and Christiansen 2009 T. Rocha-Rinza and O. Christiansen, _Chem. Phys. Lett._ , 2009, 482, 44–49.
  * Schmies _et al._ 2011 M. Schmies, A. Patzer, M. Fujii and O. Dopfer, _Phys. Chem. Chem. Phys._ , 2011, 13, 13926–13941.
  * Kołaski _et al._ 2013 M. Kołaski, C. Arunkumar and K. S. Kim, _J. Chem. Theory Comput._ , 2013, 9, 847–856.
  * Cammi 2012 R. Cammi, _Int. J. Quantum Chem._ , 2012, 112, 2547–2560.
  * Caricato 2018 M. Caricato, _Int. J. Quantum Chem._ , 2018, e25674.
  * Caricato 2019 M. Caricato, _Int. J. Quantum Chem._ , 2019, 119, e25710.
  * Caricato 2020 M. Caricato, _Wiley Interdiscip. Rev. Comput. Mol. Sci._ , 2020, 10, e1463.
  * Mobley _et al._ 2007 D. L. Mobley, É. Dumont, J. D. Chodera and K. A. Dill, _J. Phys. Chem. B_ , 2007, 111, 2242–2254.
  * Giovannini _et al._ 2019 T. Giovannini, A. Puglisi, M. Ambrosetti and C. Cappelli, _J. Chem. Theory Comput._ , 2019, 15, 2233–2245.
  * Giovannini _et al._ 2019 T. Giovannini, L. Grazioli, M. Ambrosetti and C. Cappelli, _J. Chem. Theory Comput._ , 2019, 15, 5495–5507.
  * Marrazzini _et al._ 2020 G. Marrazzini, T. Giovannini, F. Egidi and C. Cappelli, _J. Chem. Theory Comput._ , 2020, 16, 6993–7004.


0.24 0.24 0.24 0.24

Figure :  Figure : Qt Figure : Android Figure : OpenStack Figure :  Comparison
of perdition performance (AUC) between the specific models and the general
models. The model with a prefix “G” means using the trained general model to
predict changes with a particular intent. “b” represents ‘Bug Fix’, “re”
represents ‘Resource’, “f” represents ‘Feature’, “t” represents ‘Test’, “r”
represents ‘Refactor’, “m” represents ‘Merge’, “d” represents ‘Deprecate’, and
“o” represents ‘Others’. For example G-b means using the trained general model
to predict changes with the ‘Bug Fix’ intent. S-b is the specific model
trained and evaluated on changes with the ‘Bug Fix’ intent.