# Assessing the Quality of QM/MM Approaches to Describe Vacuo-to-water
Solvatochromic Shifts

Luca Nicoli  Tommaso Giovannini<EMAIL_ADDRESS>Chiara Cappelli
<EMAIL_ADDRESS>Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126
Pisa, Italy.

###### Abstract

The performance of different Quantum Mechanics/Molecular Mechanics embedding
models to compute vacuo-to-water solvatochromic shifts are investigated. In
particular, both non-polarizable and polarizable approaches are analyzed and
computed results as compared to reference experimental data. We show that none
of the approaches outperforms the others and that errors strongly depend on
the nature of the molecular transition. Thus, we prove that the best choice of
embedding model highly depends on the molecular system, and that the use a
specific approach as a black-box can lead to significant errors and sometimes
totally wrong predictions.

††preprint: AIP/123-QED

##  I Introduction

Focused models have a long standing tradition in computational chemistry for
the simulation of spectral properties of complex systems. Warshel and Levitt
(1976); Senn and Thiel (2009); Giovannini, Egidi, and Cappelli (2020a);
Cappelli (2016) Among them, quantum mechanics/molecular mechanics (QM/MM)
approaches have become very popular, Warshel and Levitt (1976); Gao (1997);
Lin and Truhlar (2007); Senn and Thiel (2009); Mennucci and Corni (2019) due
to their strengths in dealing with many diverse external environments, ranging
from strongly interacting solventsGiovannini, Egidi, and Cappelli (2020a) to
biomolecular environments.Bondanza  _et al._ (2020); Curutchet and Mennucci
(2017); Cupellini  _et al._ (2019) Indeed, the increasing popularity of QM/MM
is linked to their ability to describe target/environment interactions with an
atomistic detail.Senn and Thiel (2009); Boulanger and Harvey (2018)

When applied to solvated systems, the most common QM/MM partition consists of
treating the solute at the QM level, and the solvent in terms of classical MM
force fields. For a given QM level, the quality of QM/MM results strongly
depends on the physics lying behind the specific approach which is exploited
to model the interaction between the QM and MM layers.Mennucci and Corni
(2019) The latter is generally limited to electrostatic terms, being non-
electrostatic contributions only rarely taken into account.Giovannini,
Lafiosca, and Cappelli (2017); Giovannini  _et al._ (2019a); Giovannini,
Ambrosetti, and Cappelli (2019); Olsen  _et al._ (2015)

The MM layer can be modeled in terms of a set of fixed multipoles placed at
atomic sites, thus yielding the so-called Electrostatic Embedding (EE)
approach. Senn and Thiel (2009) As a consequence, the MM layer polarizes the
QM density, but not viceversa. However, a correct physical description of an
interacting solute-solvent systems requires mutual solute-solvent polarization
effects to be considered.Jensen, Van Duijnen, and Snijders (2003a); Boulanger
and Thiel (2012, 2014) Thus, many different polarizable embedding have been
proposed and amply tested. Jensen, Van Duijnen, and Snijders (2003a); Loco
_et al._ (2016); Olsen and Kongsted (2011); Boulanger and Thiel (2012, 2014);
Giovannini  _et al._ (2019b); Bondanza  _et al._ (2020); Wu  _et al._ (2017);
Giovannini, Egidi, and Cappelli (2020a); Loco  _et al._ (2021); Bondanza  _et
al._ (2020)

In polarizable QM/MM approaches, MM fragments are endowed with polarizable
multipolar charge distributions which are modified as a result of the
interaction with the QM density, and viceversa.Mennucci and Corni (2019) The
physically consistent description which is then obtained, permits to compute
remarkably accurate values of many spectroscopic signals, especially when
polarizable QM/MM approaches are coupled to accurate procedures to sample the
configurational phase-space. Jensen and Van Duijnen (2005); Caricato  _et al._
(2013); Giovannini  _et al._ (2019c); Segatta  _et al._ (2019); Nottoli  _et
al._ (2021) The various QM/MM approaches differ from the specific way the
electrostatic and polarization terms are modeled. The latter not only modifies
the solute’s ground state density, but also its response properties.

Despite the increasing interest in exploiting QM/MM approaches to describe
spectral properties, the performance of the different QM/MM approaches has
only rarely been investigated.Reinholdt, Nørby, and Kongsted (2018); Prioli
and Kongsted (2021) Therefore, the ideal model to be employed for a given
application has not been clearly defined yet.

In this work, we present extensive comparison of the results obtained by
applying a selection of QM/MM embedding models to the calculation of vacuo-to-
water solvatochromic shifts. The approaches are chosen because they
conceptually span diverse classes of models that are employed in the
literature. In particular, we employ the EE (as specified by means of the
TIP3P parametrization), Jorgensen  _et al._ (1983) where MM atoms are
described in terms of fixed charges, the polarizable Fluctuating Charges (FQ),
Rick, Stuart, and Berne (1994); Chen and Martínez (2009); Chen, Hundertmark,
and Martínez (2008); Giovannini, Egidi, and Cappelli (2020a) where
polarization effects are described in terms of a set of charges that vary as a
response to the external electric potential.Cappelli (2016) Discrete Reaction
Field (DRF) Jensen, Van Duijnen, and Snijders (2003a, b) is an example of
amply used approaches which model polarization effects in terms of a set of
induced dipoles assigned to MM atoms.Jensen, Van Duijnen, and Snijders (2003a,
b); Olsen and Kongsted (2011); Curutchet  _et al._ (2009); List, Olsen, and
Kongsted (2016) More sophisticated models are used to refine DRF electrostatic
description in terms of fixed multipolar expansions. Loco  _et al._ (2016);
Ren and Ponder (2002, 2003) The last approach is the Fluctuating Charges and
Fluctuating Dipoles (FQF$\mu$) model,Giovannini  _et al._ (2019b) where each
MM atom is assigned a charge and dipole which can vary as a result of
polarization effects. While EE and DRF directly follow from an electrostatic
multipolar expansion of the energy,Stone (2013) FQ is grounded in conceptual
DFT,Geerlings, De Proft, and Langenaeker (2003) and FQF$\mu$ can be seen as a
pragmatical extension of FQ.Giovannini  _et al._ (2019b)

Each embedding approach models QM/MM interactions according to the order of
the multipolar expansion of the MM variables (fixed and/or polarizable). From
the numerical point of view, such an interaction also depends on the
parameters defining the specific model: fixed atomic charge $q$ (for EE and
DRF), atomic electronegativity $\chi$ and chemical hardness $\eta$ (for FQ and
FQF$\mu$), and atomic polarizability $\alpha$ (for DRF and FQF$\mu$). The
numerical values of such parameters clearly determine the QM/MM interaction,
and in turn computed spectroscopic signals.Ambrosetti  _et al._ (2021) Thus,
in this paper a total of eight different parameter sets, which are
specifically developed for the aqueous environment, are compared. Jorgensen
_et al._ (1983); Rick, Stuart, and Berne (1994); Carnimeo, Cappelli, and
Barone (2015); Giovannini  _et al._ (2019a); Jensen, Van Duijnen, and Snijders
(2003a); Giovannini  _et al._ (2019b)

The manuscript is organized as follows. In the next section, we briefly recap
the theoretical foundations of the QM/MM embedding approaches which are
exploited in this work. Then, their performance are tested to describe vacuo-
to-water solvatochromic shift of a set of 11 medium-to-large molecules, for
which experimental data are available in the literature. Results are also
discussed in terms of the physico-chemical description of the QM/MM
interaction and the nature of the solute’s transition.

##  II Theoretical Modelling

The total energy of a QM/MM system reads: Giovannini, Egidi, and Cappelli
(2020a)

$\displaystyle E^{tot}$ $\displaystyle=E_{QM}+E_{MM}+E_{QM/MM}^{int}$ (1)

where $E_{QM}$ and $E_{MM}$ are the energies of the QM and MM portions,
respectively. By neglecting non-electrostatic (dispersion/repulsion)
interactions, the QM-MM interaction energy $E^{int}_{QM/MM}$ can be expressed
as:

$\displaystyle E_{QM/MM}^{int}$
$\displaystyle=E_{QM/MM}^{ele}+E_{QM/MM}^{pol}$ (2)

where the electrostatic $E_{QM/MM}^{ele}$ and possibly polarization
$E_{QM/MM}^{pol}$ energy terms are highlighted. In a generic definition of a
force field, MM atoms can be endowed with a fixed multipolar distribution
$\mathbf{M}$ (charges, dipoles, quadrupoles, …) and additional quantities
$\mathbf{D}$, accounting for polarization effects. By this, the various
polarizable or non-polarizable QM/MM approaches differ in the way they define
$\mathbf{M}$ and $\mathbf{D}$, and because they account or neglect
polarization terms (i.e. $\mathbf{D}$). By assuming a classical electrostatic
interaction between the QM and MM portions, the total energy in Eq. 1 can be
rewritten as:

$\displaystyle\mathcal{E}^{tot}[\rho,\mathbf{D}]=$ $\displaystyle
E_{QM}[\rho(\mathbf{r})]+\mathbf{M}^{\dagger}\int\mathbf{T}_{M}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}$
(3)
$\displaystyle+\frac{1}{2}\mathbf{D}^{\dagger}\mathbf{A}\mathbf{D}+\mathbf{D}^{\dagger}\int\mathbf{T}_{D}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}+\mathbf{D}^{\dagger}\mathbf{T}\mathbf{M}$

where the $\mathbf{A}$ matrix describes the self interaction of the
polarization sources; $\mathbf{T}$ is a block matrix, which takes into account
the interaction between the fixed and polarizable MM distributions.
$\mathbf{T}_{\xi}(\mathbf{r})$ ($\xi=M,D$) collects QM/MM electrostatic
interaction kernels Stone (2013) (see Sec. S1.1 in the Supplementary Material
– SM for more details).

Within a Kohn–Sham (KS) density functional theory (DFT) formulation, by
differentiating Eq. 3 with respect to $\rho$, the QM/MM Fock Matrix
$\tilde{F}$ is recovered. By minimizing Eq. 3 with respect to $\mathbf{D}$,
the equations which describe the polarization of the MM portion are obtained.
This allows us to define the coupled QM/MM equations:

$\displaystyle\frac{\delta\mathcal{E}^{tot}[\rho,\mathbf{D}]}{\delta\rho(\mathbf{r})}=h_{KS}^{0}[\rho(\mathbf{r})]+\hat{v}^{emb}(\mathbf{r})=\tilde{F}$
(4)
$\displaystyle\frac{\delta\mathcal{E}^{tot}[{\rho},\mathbf{D}]}{\delta\mathbf{D}}=\bm{\Theta}[{\rho},\mathbf{D}]=0$
(5)

where $h^{0}_{KS}$ is the common KS operator, given by:

$\displaystyle h_{\mathrm{KS}}^{0}$
$\displaystyle=-\frac{1}{2}\nabla^{2}-\sum_{m}\frac{Z_{m}}{\left|\mathbf{r}-\mathbf{R}_{m}\right|}+\int\frac{\rho(\mathbf{r}^{\prime})}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}d\mathbf{r}^{\prime}+\frac{\delta
E_{\mathrm{XC}}}{\delta\rho(\mathbf{r})}$ (6)

where $E_{XC}$ is the exchange-correlation energy functional. In Eqs. 4 and 5,
$\hat{v}^{emb}(\mathbf{r})$ and $\bm{\Theta}[\rho,\mathbf{D}]$ are defined as:

$\displaystyle\hat{v}^{emb}(\mathbf{r})$
$\displaystyle=\mathbf{M}^{\dagger}\mathbf{T}_{M}(\mathbf{r})+\mathbf{D}^{\dagger}\mathbf{T}_{D}(\mathbf{r})$
(7) $\displaystyle\bm{\Theta}[\rho,\mathbf{D}]$
$\displaystyle=\mathbf{A}\mathbf{D}+\int\mathbf{T}_{D}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}+\mathbf{T}\mathbf{M}$
(8)

The solutions of Eqs. 4 and 5 define the ground state (GS) QM density and the
polarization vector $\mathbf{D}$.

Vertical excitation energies can be computed by resorting to the linear
response (LR) formulation of the time-dependent DFT (TDDFT) formalism (see
Sec. S1 in the SM for more details). Olsen and Kongsted (2011); Giovannini,
Egidi, and Cappelli (2020a); Jensen, Van Duijnen, and Snijders (2003b);
Giovannini  _et al._ (2019d). In the case of QM/MM approaches, LR-TDDFT
equations Casida  _et al._ (1998) are modified to account for the presence of
the MM layer.Giovannini  _et al._ (2019d) In particular, the MM environment
modifies excitation energies through two mechanisms: (i) modification of
energy and spatial distribution of GS molecular orbitals (MOs), usually
referred to as indirect effect; (ii) inclusion of additional terms in LR-TDDFT
equations, which account for the mutual interaction between the MM layer and
the transition QM density. The latter contribution is usually called “direct
effect”, and is only present in case of polarizable embedding approaches. Note
that state-specific formulations of polarizable embedding models have also
been proposed. They specifically account for the relaxation of the solute
density in the excited state of interest, while discarding the dynamical
aspects of solute-solvent interactions, which are instead considered in the LR
formalism.Giovannini  _et al._ (2019d); Caricato  _et al._ (2006) It is
finally worth noting that local field effects induced on the QM moiety due to
the polarization of the MM portion to the external radiation field, are not
taken into account in this work, although they may affect computed oscillator
strengths.List, Jensen, and Kongsted (2016)

###  II.1 Embedding Models

The equations reported in the previous section are general enough to
constitute a unified framework, which can be specified for the various
embedding approaches that are exploited in the present work. The latter differ
in the way $\mathbf{D}$, $\mathbf{M}$ are defined.

  1. 1.

Electrostatic embedding (EE): each atom in the MM region is endowed with a
fixed charge i.e. $\mathbf{M}=[\mathbf{q}_{M}]$ and $\mathbf{D}=[\mathbf{0}]$.
Therefore, the MM layer polarizes the QM density but not viceversa, thus it
only indirectly affects the QM solute’s response properties.

  2. 2.

Fluctuating Charges (FQ) approach: each MM atom is endowed with a charge,
whose value is not fixed, but varies as a result polarization effects. Rick
_et al._ (1995); Cappelli (2016); Giovannini, Egidi, and Cappelli (2020a, b)
Thus, $\mathbf{M}=[\mathbf{0}]$ and $\mathbf{D}=\left[\mathbf{q}\right]$, with
$\mathbf{q}$ being the polarizable charges. The parameters entering the FQ
models, thus determining the $\mathbf{q}$ charges, are the atomic
electronegativity $\chi$ and chemical hardness $\eta$, which are theoretically
defined in conceptual DFT.Geerlings, De Proft, and Langenaeker (2003)
Polarization follows from the electronegativity equalization principle,
Mortier, Van Genechten, and Gasteiger (1985); Sanderson (1951) which allows to
define atomic partial charges in terms of the constrained minimum of a
suitable energy functional. Giovannini, Egidi, and Cappelli (2020a) More
details on the FQ model can be found in section S2 in the SM.

  3. 3.

Discrete Reaction Field (DRF): each MM atom is endowed with a fixed charge $q$
and a polarizable dipole $\bm{\mu}$, Jensen, Van Duijnen, and Snijders (2003a,
b) This approach to model polarization effects is exploited also by other
polarizable QM/MM approaches.Olsen and Kongsted (2011); Boulanger and Thiel
(2012); Mennucci and Corni (2019); Ren and Ponder (2003) Thus, in this case,
$\mathbf{M}=\mathbf{q}_{M}$ and $\mathbf{D}=[\mathbf{0},\bm{\mu}]$. Additional
details about DRF can be found in section S3 in the SM.

  4. 4.

Fluctuating Charges and Fluctuating Dipoles (FQF$\mu$): each MM atom is
endowed both with a polarizable charge $q$ and a polarizable dipole
$\bm{\mu}$. Giovannini  _et al._ (2019b, d, e); Giovannini, Ambrosetti, and
Cappelli (2019); Marrazzini  _et al._ (2020) FQF$\mu$ is a pragmatic extension
of the FQ model, where $\mathbf{D}=\left[\mathbf{q},\bm{\mu}\right]$. The
parameters that need to be set are the atomic electronegativity $\chi$,
chemical hardness $\eta$ and atomic polarizability $\alpha$. Additional
information on FQF$\mu$ is reported in section S4 in the SM.

To better understand analogies and differences between the aforementioned
approaches, they are schematically specified for the water molecule in Fig. 1.

Figure 1: Graphical explanation of the variables associated to EE, FQ, DRF and
FQF$\mu$ FFs for the water molecule.

##  III Computational details

We apply the aforementioned QM/MM approaches to the calculation of vacuo-to-
water solvatochromic shifts. To this end, we select eleven molecules (see Fig.
2) for which experimental UV-Vis absorption spectra in aqueous solution are
available in the literature.Aidas  _et al._ (2008); Millefiori  _et al._
(1977); Kovalenko  _et al._ (2000); Novaki and El Seoud (1996); Soujanya,
Krishna, and Samanta (1992); Politzer  _et al._ (1989); Abdel-Halim and Awad
(2009); Zakerhamidi  _et al._ (2014); Reichardt and Welton (2011); Arroyo  _et
al._ (2009); Abdel-Halim and Awad (2005); Jacques (1986); Martínez-Fernández
_et al._ (2016) The variety of molecular size, together with the different
sign of experimentally measured solvatochromic shifts, makes this set an ideal
test-bed for embedding models.

Figure 2: Molecules studied in this work. I acrolein, II para-nitroaniline,
III 1-methyl-8-oxyquinolinium betaine, IV 4-aminophtalimide, V
syn-$CH_{3}$-$CH_{3}$ bimane, VI
4-2-[4-(dimethylamino)phenyl]ethenyl-1-methylpyridinium, VII doxorubicin, VIII
Reichardt’s betaine, IX 4,4-difluoro-4-bora-3a,4a-diaza-s-indacene (BODIPY), X
1-methyl-4-[(oxocyclohexadienylidene)ethylidene]-1,4-dihydropyridine or
Brooker’s merocyanine (MB), XI 5-methylcytidine.

In order to sample the solute-solvent phase-space, molecular dynamics (MD)
simulations are performed. In particular, we run MD simulations of I, II, IV,
V, VI, VII, IX, and XI in aqueous solution without imposing any constraints on
the solute’s geometry. On the other hand, the solute geometry is kept frozen
at the PCM Tomasi, Mennucci, and Cammi (2005) optimized structure during MD
runs of III, VIII, and X, because only minor geometrical distortions are
expected, due to their limited flexibility. All MDs are performed according to
the protocols previously reported by some of the present authors (see Refs.
14, 26, 42, 44). For each molecule, a series of uncorrelated snapshots are
extracted from MD runs, and for each snapshot, a spherical droplet with a
variable radius depending on the solute intrinsic size is cut. The droplet
radius ranges from 15 Å (molecule I) to 25 Å (molecule VIII), and it is set to
retain all solute-solvent interactions (see Tab. 1 for the average number of
water molecules included for each system). For each snapshot, the absorption
spectrum is calculated and convoluted with a Gaussian function of FWHM of 0.3
eV. The final UV/Vis absorption spectrum, is then obtained as the average over
the set of uncorrelated snapshots extracted from MD runs. Note that the
convergence of the computed spectra as a function of the number of snapshots
(see Tab. 1) has been checked. Similarly to previous studies,Aidas  _et al._
(2008); Giovannini, Ambrosetti, and Cappelli (2019); Giovannini  _et al._
(2019d) for each investigated molecule we perform an extra set of calculations
in which all water molecules that are placed at a distance lower than 3.5 Å
from each solute atom are described at the QM level, whereas the remaining
ones are treated by using the FQF$\mu$ force field. The resulting approach is
called QM/QMw/FQF$\mu$ (QMw). Note that, within this approach, a proper QM
description of hydrogen bonding interactions is introduced. QM/QMw/FQF$\mu$
results are obtained as an average on the minimal number of geometries (70
structures for I, 60 structures for II, and 20 structures for the remaining
molecules) which guarantee the convergence of the spectra (see Section S5.1 in
the SM).

For each system, the solvatochromic shift ($\Delta E$) is calculated as:

$\Delta E=E_{VAC}-E^{max}_{WAT}$ (9)

where $E_{VAC}$ is the excitation energy in gas-phase and $E^{max}_{WAT}$ is
the energy of the maximum of the absorption band computed in aqueous solution.
Moreover, in order to classify and to deeply investigate the nature of the
electronic transitions, we calculate the $\Delta r$ charge transfer (CT) index
by Guido et al. Guido  _et al._ (2013); Guido, Cortona, and Adamo (2014). This
index is defined in terms of the MOs involved in the electronic transition.

| IGiovannini, Ambrosetti, and Cappelli (2019) | IIGiovannini  _et al._ (2019d) | IIIAmbrosetti  _et al._ (2021) | IVGiovannini  _et al._ (2019c) | VGiovannini  _et al._ (2019c) | VIGiovannini  _et al._ (2019c) | VIIGiovannini  _et al._ (2019c) | VIIIAmbrosetti  _et al._ (2021) | IXGiovannini  _et al._ (2019c) | XAmbrosetti  _et al._ (2021) | XIGiovannini  _et al._ (2019c)
---|---|---|---|---|---|---|---|---|---|---|---
Nframes | 200 | 100 | 100 | 200 | 200 | 200 | 200 | 100 | 200 | 100 | 200
N${}_{\mathrm{H_{2}O}}$ | 492 | 491 | 868 | 725 | 729 | 1179 | 1208 | 2249 | 727 | 1534 | 730

Table 1: Number of exploited frames (Nframes) and average number of solvent
molecules (N${}_{\mathrm{H_{2}O}}$) included in the droplets used in QM/MM
calculations of each considered systems.

All QM/MM calculations are performed by using a locally modified version of
ADF adf ; Te Velde  _et al._ (2001) engine within the Amsterdam Modeling Suite
(AMS). Ams  The QM part is treated by exploiting the CAMY-B3LYP functional
Akinaga and Ten-no (2008); Seth and Ziegler (2012) combined with the TZ2P
basis set. Van Lenthe and Baerends (2003) Solvent molecules within the MM
region are described by means of the aforementioned classical force fields. In
QM electrostatic embedding calculations, the water molecules are described by
the TIP3P force field. Mark and Nilsson (2001) Three parameterizations are
employed in QM/FQ calculations: FQ1 (from Ref. 32), FQ2 (from Ref. 43) and FQ3
(from Ref. 13). Three parameterizations are also employed for QM/DRF
calculations: DRF1, DRF2 and DRF3 (all from Ref. 16). The parameterization
presented in Ref. 21 is exploited in QM/FQF$\mu$ and QMw calculations.

##  IV Results

In this section we examine vacuo-to-water solvatochromic shifts. First, we
discuss reference computed data in vacuo. Second, we compare computed
solvatochromic shifts obtained with each embedding model with their
experimental counterparts, by resorting to selected statistical estimators to
quantitatively analyze the performance of the various models. Finally, we
rationalize our findings in light of the different physico-chemical
description provided by the employed approaches. Remarkably, QM/MM results are
commonly directly compared with experimental findings, and for this reason we
primarily assess the quality of the embedding models by taking the
experimental solvatochromic shift as a reference. However, it is worth noting
that experimental shifts might be affected by a variety of physical effects
(temperature, Franck-Condon broadening, repulsion, dispersion, …) which are
only partially included (or even absent) in our QM/MM modeling. Then, QM/MM
results are also compared to QMw data, which may be seen as complementary to
experimental measurements. However, it is worth remarking that such a
reference is not optimal, because QMw provides a QM description of solute-
solvent interactions, thus including Pauli repulsion and charge transfer
effects, which are not taken into account by a purely QM/classical approach.

###  IV.1 Excitation energies in vacuo

Experimental vacuum excitation energies are only available in the literature
for I and II. In order to obtain reference "experimental" values for the other
compounds, we resort to an extrapolation procedure, through a linear fit of
the experimentally available excitation energies measured in different
solvents, as a function of the solvent polarity indicator $E^{N}_{T}$. The
value of $E^{N}_{T}$ in gas phase is set to -0.111 according to Reichardt et
al.. Reichardt (1994); Reichardt and Welton (2011) The resulting values are
reported in Tab. 2, together with the fitting $R^{2}$ coefficients (see
Sections S5.2 – S5.12 in the SM for more details).

| VAC | WTR | Shift | R2
---|---|---|---|---
molecule | Calc | Exp | Exp | Exp |
I${}_{n\to\pi^{*}}$ | 3.78 | 3.69 Aidas  _et al._ (2008) | 3.94Aidas  _et al._ (2008) | -0.25 | –
I${}_{\pi\to\pi^{*}}$ | 6.46 | 6.41Aidas  _et al._ (2008) | 5.89Aidas  _et al._ (2008) | 0.52 | –
II | 4.34 | 4.24 Millefiori  _et al._ (1977) | 3.26 Kovalenko  _et al._ (2000) | 0.98 | –
III | 2.06 | 1.98∗ | 2.80 Novaki and El Seoud (1996) | -0.82 | 0.998
IV | 4.00 | 3.50∗ | 3.35 Soujanya, Krishna, and Samanta (1992) | 0.15 | 0.537
V | 3.99 | 3.54∗ | 3.20Politzer  _et al._ (1989) | 0.34 | 0.999
VI | 2.76 | 2.64∗ | 2.78Abdel-Halim and Awad (2009) | -0.14 | 1.000
VII | 2.97 | 2.50∗ | 2.49 Zakerhamidi  _et al._ (2014) | 0.01 | 0.196
VIII | 1.61 | 1.18 Reichardt and Welton (2011) | 2.74 Reichardt and Welton (2011) | -1.56 | 1.000
IX | 3.15 | 2.46∗ | 2.49Arroyo  _et al._ (2009) | -0.03 | 0.437
X | 2.94 | 1.56∗ | 2.80 Abdel-Halim and Awad (2005); Jacques (1986) | -1.24 | 0.960
XI | 4.79 | 4.30∗ | 4.46Martínez-Fernández  _et al._ (2016) | -0.16 | 1.000
Table 2: Calculated (Calc) and experimental (Exp) excitation energies (eV) in
vacuo (VAC) and aqueous solution (WAT). The experimental values reported with
$*$ are extrapolated by a linear fit (R2 coefficients reported in the last
column). Experimental solvatochromic shifts (eV) are also given.

Acrolein (see I in Fig. 2) is experimentally characterized by a dark
$n\to\pi^{*}$ and a bright $\pi\to\pi^{*}$ transition, which are placed at
3.69 and 6.41 eV, respectively.Aidas  _et al._ (2008) Our calculations in gas-
phase are in very good agreement with such findings, reporting a first dark
transition at 3.78 eV and a second bright one at 6.46 eV (see Tab. 2). In the
following analysis, such transitions are named I${}_{n\to\pi^{*}}$ and
I${}_{\pi\to\pi^{*}}$, respectively. For para-nitroaniline (pNA, see II in
Fig. 2), $\pi\to\pi^{*}$ transition is analyzed,Sok  _et al._ (2011);
Giovannini  _et al._ (2019d) which is experimentally measured at 4.24
eV,Millefiori  _et al._ (1977) and is well reproduced by our calculations
(4.34 eV). The same occurs for 1-methyl-8-oxyquinolinium betaine (see III in
Fig. 2) for which our computed value is 2.06 eV, in good agreement with the
extrapolated experimental value of 1.98 eV (see Sec. S5.4 in the SM).

The experimental excitation energy of 4-aminophtalimide (see IV in Fig. 2) is
extrapolated at 3.50 eV (see Sec. S5.5 in the SM). Our computed value is
placed at about 4.00 eV, differing of almost 0.5 eV from the extrapolated
experimental value. Such a discrepancy is probably due to the choice of
functional/basis set which may be not ideal for such system. We remark however
that when studying solvatochromic shifts, which are obtained as energy
differences, such systematic errors should not affect much the results.Loco
_et al._ (2016)

The first transition of syn-CH3-CH3 (bimane, see V in Fig. 2), is computed at
about 3.99 eV and is experimentally extrapolated at 3.54 eV (see Sec. S5.6 in
the SM).

4-2-[4-(dimethylamino)phenyl]ethenyl-1-methylpyridinium (see VI in Fig. 2) is
characterized by a computed first transition in vacuo at 2.76 eV, in good
agreement with the corresponding experimental extrapolated excitation energy
at 2.64 eV (see Sec. S5.7 in the SM).

The first transition of doxorubicin (see VII in Fig. 2) in gas-phase falls at
3.00 eV , which differs of almost 0.5 eV from the extrapolated experimental
excitation energy (see Sec. S5.8 in the SM). Reichardt’s dye (see VIII in Fig.
2) first transition is calculated at about 1.61 eV in good agreement with its
experimental counterpart (1.175 eV, see Sec. S5.9 in the SM). By moving to
BODIPY (see IX in Fig. 2), we report its first transition in vacuo at 3.15 eV.
Its experimental excitation energy in vacuo is extrapolated at 2.46 eV (see
Sec. S5.10 in the SM). For Brooker’s merocyanine (see X in Fig. 2), the first
transition has a $\pi\to\pi^{*}$ character and its maximum is placed at 2.94
eV, which largely differs from the extrapolated experimental excitation energy
(1.56 eV, see Sec. S5.11 in the SM). Note however that the computed value is
in agreement with previously reported values (2.92 eV from Ref. 42).

Finally, we studied the first transition of 5-methylcytidine (see XI in Fig.
2) which is computed at 4.79 eV. Its experimental excitation energy in vacuo
is extrapolated at 4.30 eV (see Sec. S5.12 in the SM).

###  IV.2 Solvatochromic Shifts

Figure 3: Computed absorption spectra of molecules I–XI in gas-phase (black)
and in aqueous solution, as predicted by different embedding models.

In Fig. 3, computed spectra of I–XI (see Fig. 2) in gas-phase and aqueous
solution are reported. Raw data and pictures of the MOs involved in the
transitions are shown in Sections S5.2 – S5.12 in the SM. The data depicted in
Fig. 3 clearly show that vertical excitation energies, oscillator strengths
and band-widths, vary as a function of the specific embedding approach (and
parametrization) which is exploited. To rationalize such findings, we first
investigate vacuo-to-water solvatochromic shifts ($\Delta E$) for I – VIII, as
computed by exploiting the different embedding approaches discussed above. The
results are graphically depicted in panels I – VIII in Fig. 4 (see Tabs.
S4-S18 in the SM for the raw data).

Figure 4: Computed solvatochromic shifts ($\Delta E$, in eV) with different
QM/MM methods (see key). Experimental reference values are also plotted. Red
asterisk denotes that the reference experimental gas-phase excitation energy
has been extrapolated by a linear fit (see text).

All QM/MM approaches are able to qualitatively grasp experimental
solvatochromic shifts, independently of the nature of the shift, i.e.
batochromic (I${}_{\text{b}}$, II, IV, V, VII) or hypsochromic
(I${}_{\text{a}}$, III, VI, VIII). It is worth noting that the studied
molecules display large variability in molecular size and, more importantly,
amplitude of the solvatochromic shifts, which range from 0.01 eV (VII) to
-1.56 eV (VIII).

We now focus on compounds IX, X, XI (see Tabs. S19-S24 in the SM for the
corresponding raw data). For IX, all embedding approaches and the reference
QMw method report a slight bathochromic shift. This is not in agreement with
the small hypsochromic shift of -0.032 eV expected from our extrapolated
linear fit. However, it is worth noting that, as previously reported, Momeni
and Brown (2015) TD-DFT fails to reproduce excitation energies of such a
system. Therefore, the systematic failure of all embedding models to reproduce
even the qualitative nature of the solvatochromic shift for IX may be probably
due to such an incorrect description of the electronic transition provided by
TD-DFT. Moving to X, the experimental hypsochromic shift is correctly
described by FQ2, FQ3 and FQF$\mu$, but not by the other approaches (see Sec.
S5.11 on the SM). Remarkably, FQ2, FQ3 and FQF$\mu$ are able to grasp the sign
of the shift, although the absolute value is underestimated, of about 1.23 eV
(FQ2), 0.91 (FQ3) and 0.760 eV (FQF$\mu$). TIP3P, FQ1 and all DRF
parametrizations incorrectly predict the solvatochromic shift, both in
absolute value and sign. A similar trend is also reported for XI: FQ3 and FQFμ
perform much better than the other embedding approaches, by underestimating
the fitted shift of 0.047 eV and 0.036 eV, respectively (see Sec. S5.12 on the
SM). In this case, only FQ1 incorrectly describes the absolute value and sign
of the experimental shift.

###  IV.3 Discussion

The quality of the results reported in the previous section, and therefore the
reliability and accuracy of the different embedding approaches, can be further
quantified by means of the following statistical estimators (SE):

Mean Absolute Error (MAE)
$\displaystyle=\frac{\sum_{i}|x_{i}-\bar{x}_{i}|}{N}$ (10) Mean Error (ME)
$\displaystyle=\frac{\sum_{i}x_{i}-\bar{x}_{i}}{N}$ (11) Mean Error unsigned
(MEu) $\displaystyle=\frac{\sum_{i}|x_{i}|-|\bar{x}_{i}|}{N}$ (12) Max MEu
(MMEu) $\displaystyle=\text{max}\left(|x_{i}|-|\bar{x}_{i}|\right)$ (13)

where $x_{i}$ and $\bar{x}_{i}$ are the computed and the reference
solvatochromic shifts (either the experimental values or the computed QMw
ones) of the $i$-th transition (see Fig. 4), respectively, whereas $N$ is the
total number of transitions considered.

In the following, we report computed SE values for the whole set of
investigated systems by taking as reference experimental shifts (see Tab. 3)
and calculated QMw data (see Tab. 4).

| QMw | TIP3P | FQ1 | FQ2 | FQ3 | DRF1 | DRF2 | DRF3 | FQF$\mu$
---|---|---|---|---|---|---|---|---|---
MAE | 0.29 | 0.33 | 0.37 | 0.32 | 0.32 | 0.30 | 0.33 | 0.33 | 0.36
ME | 0.25 | 0.25 | 0.27 | 0.24 | 0.13 | 0.27 | 0.28 | 0.28 | 0.08
MEu | 0.07 | -0.11 | -0.17 | -0.11 | 0.13 | -0.02 | -0.06 | -0.06 | 0.22
MMEu | -1.16 | -1.24 | -1.18 | -1.23 | -0.91 | -1.17 | -1.14 | -1.14 | 0.78
succ % | 92% | 83% | 75% | 92% | 92% | 83% | 83% | 83% | 92%
Table 3: Calculated MAE, ME, MEu, MMEu (eV) for all systems, as described by the various embedding approaches. Reported SE are computed with respect to reference experimental shifts. The success rate in percentage (succ %) of each model in predicting the correct experimental sign of the shift is also reported. | TIP3P | FQ1 | FQ2 | FQ3 | DRF1 | DRF2 | DRF3 | FQF$\mu$
---|---|---|---|---|---|---|---|---
MAEQ | 0.18 | 0.25 | 0.18 | 0.13 | 0.10 | 0.14 | 0.14 | 0.20
MEQ | 0.00 | 0.02 | -0.01 | -0.12 | 0.02 | 0.03 | 0.03 | -0.17
ME${}_{u}^{\mathrm{Q}}$ | -0.18 | -0.24 | -0.18 | 0.06 | -0.08 | -0.12 | -0.12 | 0.15
MME${}_{u}^{\mathrm{Q}}$ | -0.57 | -0.75 | -0.53 | 0.31 | -0.31 | -0.45 | -0.45 | 0.71
succ % | 92 % | 83% | 100 % | 100 % | 92 % | 92 % | 92 % | 100 %
Table 4: Calculated MAE, ME, MEu, MMEu (eV) for all systems, as described by
the various embedding approaches. Reported SE are computed with respect to
calculated QM/QMw/FQF$\mu$ shifts. The success rate in percentage (succ %) of
each model in reproducing the sign of the computed QM/QMw/FQF$\mu$ shift is
also reported.

By first focusing on Tab. 3, we first note that MAE values range from 0.30 eV
(DRF1) to 0.37 eV (FQ1) with 0.29 eV for QMw, indicating that all models
similarly behave as far as absolute errors are considered. ME values confirm
the best performance of FQ3 and FQF$\mu$. However, low ME values may also
reveal error cancellations, i.e. compensation of over- and underestimated
data. In line with QMw results, FQF$\mu$ and FQ3 are the only that predict
positive MEu values: in particular the largest MEu value is reported for
FQF$\mu$ and both models overestimate experimental solvatochromic shifts. The
opposite is instead valid for all other embedding approaches, which tend to
underestimate reference values. Among all methods, the lowest MEu values are
associated to the DRF family: in particular for DRF1 ME${}_{u}=0.02$ eV, i.e.
lower than the chemical accuracy (1 kcal/mol). Low MEu values are certainly
related to a good average reproduction of the experimental solvatochromic
shifts. However, the MMEu estimator reveals that large maximum errors can be
reported by all methods, including DRF. In this case, FQF$\mu$ (followed by
FQ3) gives the lowest MMEu. Moreover, it is worth noting that only FQF$\mu$,
FQ2, and FQ3 can correctly predict the sign of the experimental solvatochromic
shifts for all molecules (except for IX), reporting the same success rate of
QMw calculations (see succ % in Tab. 3). To further rationalize our findings,
we now discuss the SE computed by taking as a reference the QMw shifts
(indicated with the superscript “Q”, see Tab. 4). As stated above, QMw
accounts for effects which are not considered in our QM/MM calculations
(mainly charge transfer between solute and solvent molecules, and Pauli
repulsion). We first note that MAEQ is bound to lower values, ranging from
0.10 eV (DRF1) to 0.25 eV (FQ1), as compared to MAE (see Tab. 3). The lowest
ME${}_{u}^{\mathrm{Q}}$ are reported for FQ3 with 0.06 eV and DRF1 with -0.08
eV, thus highlighting that, on average, FQ3 overestimates the reference QMw
shifts, whereas DRF1 underestimates them. Notably, as shown by some of the
present authors in Ref. 14, the inclusion of Pauli repulsion in polarizable
QM/MM calculations, might reduce, in absolute value, the computed
solvatochromic shift and thus may compensate possible overestimation of the
shifts with respect to QMw calculations. For this reason, we can argue that
the inclusion of Pauli repulsion in QM/MM calculations might improve FQ3 and
FQF$\mu$ performance and worsen those of all other models. Indeed, this can be
justified by the fact that FQ3 and FQF$\mu$ parametrizations aim to reproduce
electrostatic and polarization solute-solvent interactions only.Giovannini
_et al._ (2019a, b) Therefore, overestimation, in absolute value, of
solvatochromic shifts is expected. Finally, an overall reduction of the
MME${}^{Q}_{u}$ indicator compared to the values in Tab. 3 is obtained, being
the lowest values reported once again for FQ3 (0.31 eV) and DRF1 (-0.31 eV).

The results discussed in terms of the statistical estimators can be better
rationalized by focusing on two specific systems, i.e. quinolinium (III) and
Brooker’s merocyanine (X), characterized by an experimental hypsochromic
shift. Computed absorption spectra, together with the main MOs involved in the
transition (HOMO and LUMO), are reported in Fig. 5 left panel for molecule III
and right panel for molecule X. The reported hypsochromic shift in polar
solvents Reichardt (2004) of III is related to larger stabilization of the
dipolar charge-separated ground state with respect to the non-zwitterionic
excited state Murugan (2011); Kumoi  _et al._ (1970). A similar explanation
has also been reported to explain the large hypsochromic shift of X when
dissolved in polar solventsWada, Nakano, and Sato (2014); Morley  _et al._
(1997); Wada, Nakano, and Sato (2014); Abdel-Halim and Awad (2005); Tanaka
_et al._ (2017).

Figure 5: Calculated absorption spectra of III (left) and X (right) by
exploiting the different QM/MM approaches. The calculated spectrum in gas-
phase (black) and the corresponding HOMO and LUMO orbitals are also
graphically depicted.

The left panel of Fig. 5 shows that all embedding approaches correctly predict
the hypsochromic shift, with an associated enhancement of the absorption
signal. In particular, the larger the computed solvatochromic shift, the
higher enhancement is reported. Oppositely, only FQ2, FQ3 and FQF$\mu$
correctly reproduce the correct sign of the experimental shift of X, also
reporting a significant hypochromic effect (see Fig. 5 right panel). Such
findings can be directly related to the physico-chemical differences of the
embedding approaches. In fact, by first focusing on the FQ family, we recall
that the FQ force field is defined in terms of $\chi$ and $\eta$, the former
($\chi$) being the source (together with the QM density) of charge
redistribution, the latter ($\eta$) defining charge self-interaction. Thus,
lower $\eta$ values and larger $\chi$ differences between different MM atoms
imply that, at equilibrium, larger charges arise on the FQ atoms. As a
consequence, the QM/MM electrostatic interaction is larger. Thus, larger
distortion of the QM MOs, with an associated increased dipole moment, is
expected. The reported FQ trends can then be rationalized in terms of computed
ground state dipole moments (see Tab. S25 in the SM). In fact, FQ3 enhances
QM/MM site-specific interactions, possibly deforming GS orbitals more than
FQ2, and in turn FQ1.

By moving to the DRF family, we note that DRF2 and DRF3 give almost the same
results, which are in agreement with previous studies Jensen, Van Duijnen, and
Snijders (2003a). DRF1 reports slightly larger solvatochromic shift, because
it exploits lower oxygen polarizability ($\alpha$) and higher hydrogen
$\alpha$. Then, the electrostatic QM/MM interaction is overall increased (see
Tab. S26 in the SM), thus resulting in higher shifts with respect to other DRF
parametrizations. Nevertheless, estimated vacuo-to-water shifts show only
minor differences among different DRF models. This is in line with Ref. 16,
and confirms that the dissection of water polarizability into atomic
contributions does not affect much computed transition energies. Thus, the
values of the fixed charges assigned to DRF atoms crucially determine vertical
transition energies.

FQF$\mu$ parametrization is characterized by lower oxygen and hydrogen
polarizabilities as compared to DRF1. However, FQF$\mu$ also considers
polarizable charges. As a consequence, it provides the most intense QM/MM
interaction, as confirmed by computed dipole moments (see Tab. S25 and Tab.
S26 in the SM). It is also worth noting that stronger QM/MM interaction is
also reflected in a larger computed full-width-at-half-maximum (FWHM), which
is directly related to a broader distribution of excitation energies (see also
Sections S5.2 – S5.12 in the SM). The discussed trends can also be
rationalized in terms of the average dipole moment of water in the liquid
phase as modeled by exploiting the various embedding approaches (see Sec.
S5.14 and Tab. S26 in the SM).

This analysis of molecules III and X, allows us to conclude that the
contribution of charges substantially affects the electronic response. This is
not surprising and has been previously reported in other contexts.Giovannini
_et al._ (2019b); Mao  _et al._ (2017) This rationalizes the rather good
performance of the non-polarizable TIP3P in terms of statistical estimators
(see Tab. 3 and Tab. 4). However, TIP3P yields the largest MMEu and
MME${}_{u}^{Q}$, thus highlighting the importance of polarization effects.

In order to obtain a physico-chemical rationalization of the discrepancies in
the estimated blue-shifts for X and XI, we can exploit a simple model to
theoretically explain the measured hypsochromic effects. Let us assume that
both GS and ES are a linear combination of two states $|A\rangle$ and
$|B\rangle$, i.e. $|\psi\rangle=c^{\psi}_{A}|A\rangle+c^{\psi}_{B}|B\rangle$,
where $|\psi\rangle$ is either GS or ES. Let us now assume that $|A\rangle$ is
a charge separated state, whereas $|B\rangle$ is a state diffused over the
molecule. We consider two different cases:

  1. 1.

In vacuo, $c^{GS}_{A}\gg c^{GS}_{B}$ and $c^{ES}_{A}\ll c^{ES}_{B}$. Thus, GS
is charge separated, whereas ES is mostly diffused on the whole molecule. We
then expect large GS dipole moment ($\mu_{GS}$) and transition charge transfer
(CT) indices. For example, in vacuo $\mu_{GS}=6.42$ Debye, $\Delta r=1.96$ Å
for III, and $\mu_{GS}=14.47$ Debye, $\Delta r=4.48$ Å for VIII. When the
molecule is dissolved in water, GS charge separation, together with
$\mu_{GS}$, increases. Overall, electronic transitions of this type display a
substantial blue shift (I${}_{n\to\pi^{*}}$, III, VI and VIII). For these
molecules, FQ3 and FQF$\mu$ give the largest errors and consistently
overestimate the shift. On the contrary, DRF1 gives the best values and the
other DRF models yield errors of $\sim$ 10 % higher than DRF1 for these
transitions. Therefore, embedding models that stabilize the GS are needed to
grasp the blue-shifting nature of these transitions. However, models which are
parameterized to provide strong QM/MM electrostatic interactions, such as FQ3
and FQFμ, might lead to unphysical GS overstabilization.

  2. 2.

In vacuo, $c^{GS}_{A}\approx c^{GS}_{B}$ and $c^{ES}_{A}\approx c^{ES}_{B}$.
Thus, GS and ES are a superposition of both $|A\rangle$ and $|B\rangle$
states. Due to the similarities between GS and ES, we expect a small CT index
in vacuo. For example, in vacuo $\mu_{GS}=15.60$ Debye, $\Delta r=0.70$ Å for
X, and $\mu_{GS}=5.61$ Debye, $\Delta r=0.64$ Å for XI. In aqueous solution,
charge separation increases, and also $\mu_{GS}$ and $c^{GS}_{A}$. If the
solvent removes the degeneracy between $|A\rangle$ and $|B\rangle$, GS and ES
are expected to have substantially different charge distributions. In this
case, the CT index is large, and the transition is blue-shifted due to a
strong GS stabilization. This is the case for X (see Tab. S25 in the SM). On
the contrary, if both $c^{GS}_{A}$ and $c^{GS}_{B}$ increase in a similar
manner, blue-shift is expected if $c^{GS}_{A}>c^{ES}_{A}$, i.e. if GS is more
stabilized than ES. In this case, the solvatochromic shift is small and the CT
index remains similar to its value in vacuo. A transition of this type is
reported for XI (see Tab. S25 in the SM). For both transition types, FQF$\mu$
gives the lowest errors, followed by FQ3. This suggests that, for this class
of transitions, models that provide the strongest QM/MM interaction are needed
to grasp the correct electronic reorganization of the states involved.

Although qualitative, this analysis allows us to rationalize the behaviour of
blue-shifting compounds in terms of the $\Delta r$-index calculated in vacuo.
Indeed, large $\Delta r$-index calculated in vacuo (i.e. $\Delta r>1.5$ Å)
tracks a transition with a strong CT character Guido  _et al._ (2013), thus
indicating that GS and ES have substantial different charge distributions. In
Fig. 6 the absolute deviation of the computed shift from the estimated shift
($|\Delta E^{calc}-\Delta E^{exp}|$) is plotted with respect to calculated
$\Delta r$-index in vacuo for each embedding model. Absolute errors less than
zero indicate that the shift is underestimated with respect to the experiment.
In Fig. 6 a) the results obtained for the blue shifting compounds are
depicted. If gas-phase $\Delta r<1.5$ Å, all models underestimate the shifts
but FQF$\mu$ and FQ3 report the lowest absolute errors. On the other hand, if
$\Delta r>1.5$ Å, FQF$\mu$ and FQ3 consistently overestimate the shifts,
whilst FQ1 reports the largest underestimation for each transition. In this
case, DRF1 reports the lowest absolute errors. In Fig. 6 b), the results
obtained for red shifting compounds are instead considered. Interestingly, no
particular correlation between the absolute estimation error and the $\Delta
r$-index can be found. We notice however that for all molecules, DRF models
compute similar shifts, FQF$\mu$ and FQ3 yield the largest shifts, whereas FQ1
gives the smallest shifts.

Therefore, our results show that none of the tested embedding model
outperforms the others, because the agreement with experiments is strongly
dependent on the nature of the transition under consideration. This is clearly
depicted in Fig. 6 a) and finds a theoretical explanation in the two-states
model described above.

Figure 6: Correlation between gas-phase $\Delta r$ indices and absolute errors
predicted by all embedding methods for blue-shifting (panel a) and red-
shifting (panel b) transitions. Note that negative errors indicate
underestimation of the shift with respect to reference experimental values.

##  V Conclusions

In this work we have applied diverse QM/MM approaches and related
parametrizations to the calculation of vacuo-to-water solvatochromic shifts.
The selected QM/MM methods are based on the non-polarizable TIP3P, and
polarizable FQ, DRF, and FQF$\mu$ force fields. For each model we have
considered various parameterizations and computed shifts have been compared to
experimentally available data and to computed QM/QMw/FQF$\mu$ values.

For I–VIII, all models are able to grasp the correct nature of the shift,
however reporting a very large MMEu. FQ3 and FQF$\mu$ consistently
overestimate the shift (except for I${}_{\pi\to\pi^{*}}$ and II), which is
consistent with the fact that such models have been parameterized by
accounting for solute-solvent electrostatic and polarization interactions
only. For IX, X and XI, FQ3 and FQF$\mu$ instead yield the best agreement with
experimental values, and remarkably, for X, only FQ2, FQ3 and FQF$\mu$ can
describe the experimentally measured hypsochromism. Moreover, when compared
with QMw computed shifts, FQ3 and DRF1, report the best performance in terms
of statistical estimators. However, we argue that the inclusion of Pauli
repulsion in QM/MM calculations would reduce in absolute value the computed
shifts, thus improving FQ3 and FQF$\mu$ results and worsening those of all
other models.Giovannini, Lafiosca, and Cappelli (2017)

The different trends have been analyzed in terms of the physico-chemical
description provided by the employed QM/MM approaches, in terms of the
predicted strength of QM/MM electrostatic interactions and QM dipole moments,
and by means of a simplified two-state model for the blue-shifting
transitions.

The results of this study reveal that the correct description of
solvatochromic shifts is a delicate task. Polarizable QM/MM approaches are
nowadays becoming a golden-standard for condensed phase simulations, however,
their use as a black-box may potentially yield completely wrong predictions as
for molecule X. Indeed, computed values not only depend on the QM/MM approach
which is exploited, but also on the specific parameterization and on the
nature of the solute’s electronic transition to be described.

## Supplementary Material

Theoretical aspects related to the QM/MM interaction integral; details on FQ
model; details on DRF model; details on FQF$\mu$ model; convergence of
QM/QMw/FQF$\mu$ excitation energy as a function of the number of frames;
computed results for I–XI; computed properties for the blue-shifting
compounds; average dipole moment for liquid water.

###### Acknowledgements.

We gratefully acknowledge the Center for High Performance Computing (CHPC) at
SNS for providing the computational infrastructure.

## Conflict of Interest

The authors have no conflicts to disclose.

## Data Availability

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

## References

## References

  * Warshel and Levitt (1976) A. Warshel and M. Levitt, “Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme,” J. Mol. Biol. 103, 227–249 (1976).
  * Senn and Thiel (2009) H. M. Senn and W. Thiel, “Qm/mm methods for biomolecular systems,” Angew. Chem. Int. Ed. 48, 1198–1229 (2009).
  * Giovannini, Egidi, and Cappelli (2020a) T. Giovannini, F. Egidi, and C. Cappelli, “Molecular spectroscopy of aqueous solutions: a theoretical perspective,” Chem. Soc. Rev. 49, 5664–5677 (2020a).
  * Cappelli (2016) C. Cappelli, “Integrated qm/polarizable mm/continuum approaches to model chiroptical properties of strongly interacting solute–solvent systems,” Int. J. Quantum Chem. 116, 1532–1542 (2016).
  * Gao (1997) J. Gao, “Energy components of aqueous solution: Insight from hybrid qm/mm simulations using a polarizable solvent model,” J. Comput. Chem. 18, 1061–1071 (1997).
  * Lin and Truhlar (2007) H. Lin and D. G. Truhlar, “Qm/mm: what have we learned, where are we, and where do we go from here?” Theor. Chem. Acc. 117, 185–199 (2007).
  * Mennucci and Corni (2019) B. Mennucci and S. Corni, “Multiscale modelling of photoinduced processes in composite systems,” Nat. Rev. Chem. 3, 315–330 (2019).
  * Bondanza  _et al._ (2020) M. Bondanza, M. Nottoli, L. Cupellini, F. Lipparini, and B. Mennucci, “Polarizable embedding qm/mm: the future gold standard for complex (bio) systems?” Phys. Chem. Chem. Phys. 22, 14433–14448 (2020).
  * Curutchet and Mennucci (2017) C. Curutchet and B. Mennucci, “Quantum chemical studies of light harvesting,” Chem. Rev. 117, 294–343 (2017).
  * Cupellini  _et al._ (2019) L. Cupellini, M. Corbella, B. Mennucci, and C. Curutchet, “Electronic energy transfer in biomacromolecules,” WIREs Comput. Mol. Sci. 9, e1392 (2019).
  * Boulanger and Harvey (2018) E. Boulanger and J. N. Harvey, “Qm/mm methods for free energies and photochemistry,” Curr. Opin. Struct. Biol. 49, 72–76 (2018).
  * Giovannini, Lafiosca, and Cappelli (2017) T. Giovannini, P. Lafiosca, and C. Cappelli, “A general route to include pauli repulsion and quantum dispersion effects in qm/mm approaches,” J. Chem. Theory Comput. 13, 4854–4870 (2017).
  * Giovannini  _et al._ (2019a) T. Giovannini, P. Lafiosca, B. Chandramouli, V. Barone, and C. Cappelli, “Effective yet reliable computation of hyperfine coupling constants in solution by a qm/mm approach: Interplay between electrostatics and non-electrostatic effects,” J. Chem. Phys. 150, 124102 (2019a).
  * Giovannini, Ambrosetti, and Cappelli (2019) T. Giovannini, M. Ambrosetti, and C. Cappelli, “Quantum confinement effects on solvatochromic shifts of molecular solutes,” J. Phys. Chem. Lett. 10, 5823–5829 (2019).
  * Olsen  _et al._ (2015) J. M. H. Olsen, C. Steinmann, K. Ruud, and J. Kongsted, “Polarizable density embedding: A new qm/qm/mm-based computational strategy,” J. Phys. Chem. A 119, 5344–5355 (2015).
  * Jensen, Van Duijnen, and Snijders (2003a) L. Jensen, P. T. Van Duijnen, and J. G. Snijders, “A discrete solvent reaction field model within density functional theory,” J. Chem. Phys. 118, 514–521 (2003a).
  * Boulanger and Thiel (2012) E. Boulanger and W. Thiel, “Solvent boundary potentials for hybrid qm/mm computations using classical drude oscillators: a fully polarizable model,” J. Chem. Theory Comput. 8, 4527–4538 (2012).
  * Boulanger and Thiel (2014) E. Boulanger and W. Thiel, “Toward qm/mm simulation of enzymatic reactions with the drude oscillator polarizable force field,” J. Chem. Theory Comput. 10, 1795–1809 (2014).
  * Loco  _et al._ (2016) D. Loco, É. Polack, S. Caprasecca, L. Lagardère, F. Lipparini, J.-P. Piquemal, and B. Mennucci, “A qm/mm approach using the amoeba polarizable embedding: from ground state energies to electronic excitations,” J. Chem. Theory Comput. 12, 3654–3661 (2016).
  * Olsen and Kongsted (2011) J. M. H. Olsen and J. Kongsted, “Molecular properties through polarizable embedding,” in  _Advances in quantum chemistry_ , Vol. 61 (Elsevier, 2011) pp. 107–143.
  * Giovannini  _et al._ (2019b) T. Giovannini, A. Puglisi, M. Ambrosetti, and C. Cappelli, “Polarizable qm/mm approach with fluctuating charges and fluctuating dipoles: the qm/fqf$\mu$ model,” J. Chem. Theory Comput. 15, 2233–2245 (2019b).
  * Wu  _et al._ (2017) X. Wu, J.-M. Teuler, F. Cailliez, C. Clavaguera, D. R. Salahub, and A. de la Lande, “Simulating electron dynamics in polarizable environments,” J. Chem. Theory Comput. 13, 3985–4002 (2017).
  * Loco  _et al._ (2021) D. Loco, L. Lagardère, O. Adjoua, and J.-P. Piquemal, “Atomistic polarizable embeddings: energy, dynamics, spectroscopy, and reactivity,” Acc. Chem. Res. 54, 2812–2822 (2021).
  * Jensen and Van Duijnen (2005) L. Jensen and P. T. Van Duijnen, “The first hyperpolarizability of p-nitroaniline in 1, 4-dioxane: A quantum mechanical/molecular mechanics study,” J. Chem. Phys. 123, 074307 (2005).
  * Caricato  _et al._ (2013) M. Caricato, F. Lipparini, G. Scalmani, C. Cappelli, and V. Barone, “Vertical electronic excitations in solution with the eom-ccsd method combined with a polarizable explicit/implicit solvent model,” J. Chem. Theory Comput. 9, 3035–3042 (2013).
  * Giovannini  _et al._ (2019c) T. Giovannini, M. Macchiagodena, M. Ambrosetti, A. Puglisi, P. Lafiosca, G. Lo Gerfo, F. Egidi, and C. Cappelli, “Simulating vertical excitation energies of solvated dyes: From continuum to polarizable discrete modeling,” Int. J. Quantum Chem. 119, e25684 (2019c).
  * Segatta  _et al._ (2019) F. Segatta, L. Cupellini, M. Garavelli, and B. Mennucci, “Quantum chemical modeling of the photoinduced activity of multichromophoric biosystems: Focus review,” Chem. Rev. 119, 9361–9380 (2019).
  * Nottoli  _et al._ (2021) M. Nottoli, L. Cupellini, F. Lipparini, G. Granucci, and B. Mennucci, “Multiscale models for light-driven processes,” Ann. Rev. Phys. Chem. 72, 489–513 (2021).
  * Reinholdt, Nørby, and Kongsted (2018) P. Reinholdt, M. S. Nørby, and J. Kongsted, “Modeling of magnetic circular dichroism and uv/vis absorption spectra using fluctuating charges or polarizable embedding within a resonant-convergent response theory formalism,” J. Chem. Theory Comput. 14, 6391–6404 (2018).
  * Prioli and Kongsted (2021) S. Prioli and J. Kongsted, “Modeling environmental effects in two-photon circular dichroism calculations,” Theor. Chem. Acc. 140, 1–12 (2021).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, “Comparison of simple potential functions for simulating liquid water,” J. Chem. Phys. 79, 926–935 (1983).
  * Rick, Stuart, and Berne (1994) S. W. Rick, S. J. Stuart, and B. J. Berne, “Dynamical fluctuating charge force fields: Application to liquid water,” J. Chem. Phys. 101, 6141–6156 (1994).
  * Chen and Martínez (2009) J. Chen and T. J. Martínez, “Charge conservation in electronegativity equalization and its implications for the electrostatic properties of fluctuating-charge models,” J. Chem. Phys. 131, 044114 (2009).
  * Chen, Hundertmark, and Martínez (2008) J. Chen, D. Hundertmark, and T. J. Martínez, “A unified theoretical framework for fluctuating-charge models in atom-space and in bond-space,” J. Chem. Phys. 129, 214113 (2008).
  * Jensen, Van Duijnen, and Snijders (2003b) L. Jensen, P. T. Van Duijnen, and J. G. Snijders, “A discrete solvent reaction field model for calculating molecular linear response properties in solution,” J. Chem. Phys. 119, 3800–3809 (2003b).
  * Curutchet  _et al._ (2009) C. Curutchet, A. Muñoz-Losa, S. Monti, J. Kongsted, G. D. Scholes, and B. Mennucci, “Electronic energy transfer in condensed phase studied by a polarizable qm/mm model,” J. Chem. Theory Comput. 5, 1838–1848 (2009).
  * List, Olsen, and Kongsted (2016) N. H. List, J. M. H. Olsen, and J. Kongsted, “Excited states in large molecular systems through polarizable embedding,” Phys. Chem. Chem. Phys. 18, 20234–20250 (2016).
  * Ren and Ponder (2002) P. Ren and J. W. Ponder, “Consistent treatment of inter-and intramolecular polarization in molecular mechanics calculations,” J. Comput. Chem. 23, 1497–1506 (2002).
  * Ren and Ponder (2003) P. Ren and J. W. Ponder, “Polarizable atomic multipole water model for molecular mechanics simulation,” J. Phys. Chem. B 107, 5933–5947 (2003).
  * Stone (2013) A. Stone,  _The theory of intermolecular forces_ (oUP oxford, 2013).
  * Geerlings, De Proft, and Langenaeker (2003) P. Geerlings, F. De Proft, and W. Langenaeker, “Conceptual density functional theory,” Chem. Rev. 103, 1793–1874 (2003).
  * Ambrosetti  _et al._ (2021) M. Ambrosetti, S. Skoko, T. Giovannini, and C. Cappelli, “Quantum Mechanics/Fluctuating Charge Protocol to Compute Solvatochromic Shifts,” J. Chem. Theory Comput. 17, 7146–7156 (2021).
  * Carnimeo, Cappelli, and Barone (2015) I. Carnimeo, C. Cappelli, and V. Barone, “Analytical gradients for mp 2, double hybrid functionals, and td-dft with polarizable embedding described by fluctuating charges,” J. Comput. Chem. 36, 2271–2290 (2015).
  * Giovannini  _et al._ (2019d) T. Giovannini, R. R. Riso, M. Ambrosetti, A. Puglisi, and C. Cappelli, “Electronic transitions for a fully polarizable qm/mm approach based on fluctuating charges and fluctuating dipoles: linear and corrected linear response regimes,” J. Chem. Phys. 151, 174104 (2019d).
  * Casida  _et al._ (1998) M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, “Molecular excitation energies to high-lying bound states from time-dependent density-functional response theory: Characterization and correction of the time-dependent local density approximation ionization threshold,” J. Chem. Phys. 108, 4439–4449 (1998).
  * Caricato  _et al._ (2006) M. Caricato, B. Mennucci, J. Tomasi, F. Ingrosso, R. Cammi, S. Corni, and G. Scalmani, “Formation and relaxation of excited states in solution: A new time dependent polarizable continuum model based on time dependent density functional theory,” J. Chem. Phys. 124, 124520 (2006).
  * List, Jensen, and Kongsted (2016) N. H. List, H. J. A. Jensen, and J. Kongsted, “Local electric fields and molecular properties in heterogeneous environments through polarizable embedding,” Phys. Chem. Chem. Phys. 18, 10070–10080 (2016).
  * Rick  _et al._ (1995) S. W. Rick, S. J. Stuart, J. S. Bader, and B. Berne, “Fluctuating charge force fields for aqueous solutions,” J. Mol. Liq. 65, 31–40 (1995).
  * Giovannini, Egidi, and Cappelli (2020b) T. Giovannini, F. Egidi, and C. Cappelli, “Theory and algorithms for chiroptical properties and spectroscopies of aqueous systems,” Phys. Chem. Chem. Phys. 22, 22864–22879 (2020b).
  * Mortier, Van Genechten, and Gasteiger (1985) W. J. Mortier, K. Van Genechten, and J. Gasteiger, “Electronegativity equalization: application and parametrization,” J. Am. Chem. Soc. 107, 829–835 (1985).
  * Sanderson (1951) R. Sanderson, “An interpretation of bond lengths and a classification of bonds,” Science 114, 670–672 (1951).
  * Giovannini  _et al._ (2019e) T. Giovannini, L. Grazioli, M. Ambrosetti, and C. Cappelli, “Calculation of ir spectra with a fully polarizable qm/mm approach based on fluctuating charges and fluctuating dipoles,” J. Chem. Theory Comput. 15, 5495–5507 (2019e).
  * Marrazzini  _et al._ (2020) G. Marrazzini, T. Giovannini, F. Egidi, and C. Cappelli, “Calculation of linear and non-linear electric response properties of systems in aqueous solution: A polarizable quantum/classical approach with quantum repulsion effects,” J. Chem. Theory Comput. 16, 6993–7004 (2020).
  * Aidas  _et al._ (2008) K. Aidas, A. Møgelhøj, E. J. Nilsson, M. S. Johnson, K. V. Mikkelsen, O. Christiansen, P. Söderhjelm, and J. Kongsted, “On the performance of quantum chemical methods to predict solvatochromic effects: The case of acrolein in aqueous solution,” J. Chem. Phys. 128, 194503 (2008).
  * Millefiori  _et al._ (1977) S. Millefiori, G. Favini, A. Millefiori, and D. Grasso, “Electronic spectra and structure of nitroanilines,” Spectrochim. Acta A 33, 21 – 27 (1977).
  * Kovalenko  _et al._ (2000) S. Kovalenko, R. Schanz, V. Farztdinov, H. Hennig, and N. Ernsting, “Femtosecond relaxation of photoexcited para-nitroaniline: solvation, charge transfer, internal conversion and cooling,” Chem. Phys. Lett. 323, 312 – 322 (2000).
  * Novaki and El Seoud (1996) L. P. Novaki and O. A. El Seoud, “Solvatochromism in pure solvents: Effects of the molecular structure of the probe,” Ber. Bunsenges. Phys. Chem. 100, 648–655 (1996).
  * Soujanya, Krishna, and Samanta (1992) T. Soujanya, T. Krishna, and A. Samanta, “The nature of 4-aminophthalimide-cyclodextrin inclusion complexes,” J. Phys. Chem 96, 8544–8548 (1992).
  * Politzer  _et al._ (1989) I. Politzer, K. Crago, D. Kiel, and T. Hampton, “The effects of $\beta$-cyclodextrin on the fluorescence, uv absorption and solubility of selected bimanes in aqueous solutions,” Anal. Lett. 22, 1567–1580 (1989).
  * Abdel-Halim and Awad (2009) S. T. Abdel-Halim and M. K. Awad, “Solvatochromism, molecular and electronic structures of trans and cis isomers of a typical styryl pyridinium cyanine dye,” J. Mol. Struct. 920, 332–341 (2009).
  * Zakerhamidi  _et al._ (2014) M. Zakerhamidi, M. Johari-Ahar, S. S. Ahmadian, and R. Kian, “Photo-physical behavior of some antitumor anthracycline in solvent media with different polarity,” Spectrochim. Acta, Part A 130, 257–262 (2014).
  * Reichardt and Welton (2011) C. Reichardt and T. Welton,  _Solvents and solvent effects in organic chemistry_ (John Wiley & Sons, 2011).
  * Arroyo  _et al._ (2009) I. J. Arroyo, R. Hu, G. Merino, B. Z. Tang, and E. Pena-Cabrera, “The smallest and one of the brightest. efficient preparation and optical description of the parent borondipyrromethene system,” J. Org. Chem. 74, 5719–5722 (2009).
  * Abdel-Halim and Awad (2005) S. T. Abdel-Halim and M. K. Awad, “Absorption, fluorescence, and semiempirical ased-mo studies on a typical brooker’s merocyanine dye,” J. Mol. Struct. 754, 16–24 (2005).
  * Jacques (1986) P. Jacques, “On the relative contributions of nonspecific and specific interactions to the unusual solvtochromism of a typical merocyanine dye,” J. Phys. Chem 90, 5535–5539 (1986), https://doi.org/10.1021/j100280a012 .
  * Martínez-Fernández  _et al._ (2016) L. Martínez-Fernández, A. Pepino, J. Segarra-Martí, A. Banyasz, M. Garavelli, and R. Improta, “Computing the absorption and emission spectra of 5-methylcytidine in different solvents: A test-case for different solvation models,” J. Chem. Theory Comput. 12, 4430–4439 (2016).
  * Tomasi, Mennucci, and Cammi (2005) J. Tomasi, B. Mennucci, and R. Cammi, “Quantum mechanical continuum solvation models,” Chem. Rev. 105, 2999–3094 (2005).
  * Guido  _et al._ (2013) C. A. Guido, P. Cortona, B. Mennucci, and C. Adamo, “On the metric of charge transfer molecular excitations: a simple chemical descriptor,” J. Chem. Theory Comput. 9, 3118–3126 (2013).
  * Guido, Cortona, and Adamo (2014) C. A. Guido, P. Cortona, and C. Adamo, “Effective electron displacements: A tool for time-dependent density functional theory computational spectroscopy,” J. Chem. Phys. 140, 104101 (2014).
  * (70) “Adf 2021.1, scm, theoretical chemistry, vrije universiteit, amsterdam, the netherlands,” .
  * Te Velde  _et al._ (2001) G. t. Te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. van Gisbergen, J. G. Snijders, and T. Ziegler, “Chemistry with adf,” J. Comput. Chem. 22, 931–967 (2001).
  * (72) “Ams 2021.1, scm, theoretical chemistry, vrije universiteit, amsterdam, the netherlands,” .
  * Akinaga and Ten-no (2008) Y. Akinaga and S. Ten-no, “Range-separation by the yukawa potential in long-range corrected density functional theory with gaussian-type basis functions,” Chem. Phys. Lett. 462, 348–351 (2008).
  * Seth and Ziegler (2012) M. Seth and T. Ziegler, “Range-Separated Exchange Functionals with Slater-Type Functions,” J. Chem. Theory Comput. 8, 901–907 (2012).
  * Van Lenthe and Baerends (2003) E. Van Lenthe and E. J. Baerends, “Optimized Slater-type basis sets for the elements 1–118,” J. Comput. Chem. 24, 1142–1156 (2003).
  * Mark and Nilsson (2001) P. Mark and L. Nilsson, “Structure and dynamics of the tip3p, spc, and spc/e water models at 298 k,” J. Phys. Chem. A 105, 9954–9960 (2001).
  * Reichardt (1994) C. Reichardt, “Solvatochromic dyes as solvent polarity indicators,” Chem. Rev. 94, 2319–2358 (1994).
  * Sok  _et al._ (2011) S. Sok, S. Y. Willow, F. Zahariev, and M. S. Gordon, “Solvent-induced shift of the lowest singlet $\pi$→ $\pi$* charge-transfer excited state of p-nitroaniline in water: An application of the tddft/efp1 method,” J. Phys. Chem. A 115, 9801–9809 (2011).
  * Momeni and Brown (2015) M. R. Momeni and A. Brown, “Why Do TD-DFT Excitation Energies of BODIPY/Aza-BODIPY Families Largely Deviate from Experiment? Answers from Electron Correlated and Multireference Methods,” J. Chem. Theory Comput. 11, 2619–2632 (2015).
  * Reichardt (2004) C. Reichardt, “Pyridinium n-phenolate betaine dyes as empirical indicators of solvent polarity: Some new findings,” Pure Appl. Chem. 76, 1903–1919 (2004).
  * Murugan (2011) N. A. Murugan, “Modeling solvatochromism of a quinolinium betaine dye in water solvent using sequential hybrid qm/mm and semicontinuum approach,” J. Phys. Chem. B 115, 1056–1061 (2011).
  * Kumoi  _et al._ (1970) S. Kumoi, K. Oyama, T. Yano, H. Kobayashi, and K. Ueno, “Spectrophotometric determination of water in organic solvents with solvatochromic dyes,” Talanta 17, 319–327 (1970).
  * Wada, Nakano, and Sato (2014) T. Wada, H. Nakano, and H. Sato, “Solvatochromic shift of brooker’s merocyanine: Hartree–fock exchange in time dependent density functional calculation and hydrogen bonding effect,” J. Chem. Theory Comput. 10, 4535–4547 (2014).
  * Morley  _et al._ (1997) J. O. Morley, R. M. Morley, R. Docherty, and M. H. Charlton, “Fundamental studies on brooker’s merocyanine,” J. Am. Chem. Soc. 119, 10192–10202 (1997).
  * Tanaka  _et al._ (2017) Y. Tanaka, Y. Kawashima, N. Yoshida, and H. Nakano, “Solvatochromism and preferential solvation of brooker’s merocyanine in water–methanol mixtures,” J. Comput. Chem. 38, 2411–2419 (2017).
  * Mao  _et al._ (2017) Y. Mao, Y. Shao, J. Dziedzic, C.-K. Skylaris, T. Head-Gordon, and M. Head-Gordon, “Performance of the amoeba water model in the vicinity of qm solutes: a diagnosis using energy decomposition analysis,” J. Chem. Theory Comput. 13, 1963–1979 (2017).