Electronic Transitions for a Fully Polarizable QM/MM Approach Based on Fluctuating Charges and Fluctuating Dipoles: Linear and Corrected Linear Response Regimes

TL;DR

This paper extends a fully polarizable QM/MM method to calculate excitation energies in solvated systems, effectively capturing solvation effects and charge redistribution with good experimental agreement.

Contribution

The work introduces a new extension of the QM/FQFμ approach for excited states, incorporating corrected linear response to better model charge equilibration in solvation.

Findings

Accurately predicts solvatochromic shifts for various molecules.

Shows good agreement with experimental excitation energies.

Demonstrates the effectiveness of the cLR regime in excited-state calculations.

Abstract

Fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles, named QM/FQF{\mu} (J. Chem. Theory Comput. 2019, 15 2233-2245), is extended to the calculation of vertical excitation energies of solvated molecular systems. Excitation energies are defined within two different solvation regimes, i.e. linear response (LR), where the response of the MM portion is adjusted to the QM transition density, and corrected-Linear Response (cLR) in which the MM response is adjusted to the relaxed QM density, thus being able to account for charge equilibration in the excited state. The model, which is specified in terms of three physical parameters (electronegativity, chemical hardness, and polarizability) is applied to vacuo-to-water solvatochromic shifts of aqueous solutions of para-nitroaniline, pyridine and pyrimidine. The results show a good agreement with their experimental counterparts, thus highlighting the potentialities of this approach.

Tables

Table 1. Table 1: pNA (CAM-B3LYP/aug-cc-pVDZ), Pyridine and Pyrimidine (M06/6-311+G(2df,2p) computed vertical excitation energies in gas phase (E vert ) given in eV. The character of the transition is reported, together with ground ( μ g subscript 𝜇 𝑔 \mu_{g} ) and excited ( μ e ​ x subscript 𝜇 𝑒 𝑥 \mu_{ex} ) molecular dipole moments (Debye). Experimental excitation energies (eV) are also shown.

Molecule	Excitation	${\mu }_g$	${\mu }_{ex}$	$E_{vert}$	Exp.
pNA	$\pi \to {\pi }^{\ast }$	7.4	12.9	4.33	4.25${}^{\text{a}}$
pyridine	$n\to {\pi }^{\ast }$	2.3	0.4	4.67	4.59,${}^{\text{b}}$ 4.63,${}^{\text{c}}$ 4.74${}^{\text{d}}$
pyrimidine	$n\to {\pi }^{\ast }$	2.4	0.6	4.17	4.18,${}^{\text{e}}$ 4.19${}^{\text{f}}$

Table 2. Table 2: pNA, Pyridine and Pyrimidine QM/PCM, QM/FQ and QM/FQF μ 𝜇 \mu ω 0 subscript 𝜔 0 \omega_{0} , LR and cLR vertical excitation energies (E vert ) and vacuo-to-water solvatochromic shifts ( Δ Δ \Delta E). "Best results" (see text) are in bold. Calculated QM/FQ and QM/FQF μ 𝜇 \mu standard errors are also reported. They are calculated as σ / N 𝜎 𝑁 \sigma/\sqrt{N} , where σ 𝜎 \sigma is the standard deviation and N 𝑁 N is the number of the snapshots used to obtain the average property. Experimental data taken from the literature are reported for comparison’s sake. All data are given in eV.

Molecule	Method	QM/PCM		QM/FQ		QM/FQF$\mu$		Exp.
		Evert	$\Delta$E	Evert	$\Delta$E	Evert	$\Delta$E	Evert	$\Delta$E
pNA	${\omega }_0$	3.92	0.41	3.74$\pm$0.03	0.59$\pm$0.03	3.47$\pm$0.03	0.86$\pm$0.03	3.25${}^{\text{a}}$	0.99
	LR	3.81	0.52	3.70$\pm$0.03	0.63$\pm$0.03	3.31$\pm$0.03	1.02$\pm$0.03
	cLR	3.84	0.49	3.71$\pm$0.03	0.62$\pm$0.03	3.37$\pm$0.03	0.96$\pm$0.03
pyridine	${\omega }_0$	4.84	-0.17	5.00$\pm$0.01	-0.33$\pm$0.01	5.21$\pm$0.01	-0.54$\pm$0.01	4.94${}^{\text{b}}$	-0.31
	LR	4.84	-0.17	5.00$\pm$0.01	-0.33$\pm$0.01	5.20$\pm$0.01	-0.53$\pm$0.01
	cLR	4.79	-0.12	4.97$\pm$0.01	-0.30$\pm$0.01	5.12$\pm$0.01	-0.45$\pm$0.01
pyrimidine	${\omega }_0$	4.31	-0.14	4.52$\pm$0.01	-0.34$\pm$0.01	4.67$\pm$0.01	-0.49$\pm$0.01	4.57${}^{\text{c}}$	-0.35 $\pm$ 0.06${}^{\text{d}}$
	LR	4.30	-0.13	4.51$\pm$0.01	-0.34$\pm$0.01	4.66$\pm$0.01	-0.48$\pm$0.01
	cLR	4.28	-0.11	4.49$\pm$0.01	-0.31$\pm$0.01	4.63$\pm$0.01	-0.45$\pm$0.01

Table 3. Table 3: pNA, Pyridine, Pyrimidine vacuo, QM/PCM, QM/FQ and QM/FQF μ 𝜇 \mu D CT values (Å) for the studied transitions. The values are calculated by using both the unrelaxed density 𝐏 K T subscript superscript 𝐏 𝑇 𝐾 \mathbf{P}^{T}_{K} and relaxed 𝐏 K Δ subscript superscript 𝐏 Δ 𝐾 \mathbf{P}^{\Delta}_{K} densities .

Density	Molecule	Transition	Vacuo	QM/PCM	QM/FQ	QM/FQF$\mu$
${𝐏}_K^T$	pNA	$\pi \to {\pi }^{\ast }$	3.14	3.11	2.60	3.35
	pyridine	$n\to {\pi }^{\ast }$	1.60	1.51	1.48	1.40
	pyrimidine	$n\to {\pi }^{\ast }$	0.99	0.99	1.00	1.00
${𝐏}_K^{\Delta }$	pNA	$\pi \to {\pi }^{\ast }$	2.05	2.26	1.66	2.22
	pyridine	$n\to {\pi }^{\ast }$	0.80	1.00	0.82	0.85
	pyrimidine	$n\to {\pi }^{\ast }$	0.54	0.73	0.58	0.64

Table 4. Table 4: pNA, Pyridine and Pyrimidine QM/PCM, QM/FQ and QM/FQF μ 𝜇 \mu GS( μ g subscript 𝜇 𝑔 \mu_{g} ) and ES ( μ e ​ x subscript 𝜇 𝑒 𝑥 \mu_{ex} ) electric dipole moments and their squared difference ( Δ ​ μ 2 = | μ e ​ x − μ g | 2 Δ superscript 𝜇 2 superscript subscript 𝜇 𝑒 𝑥 subscript 𝜇 𝑔 2 \Delta\mu^{2}=\left|{\mu_{ex}-\mu_{g}}\right|^{2} ). ( 2 ​ μ g , e ​ x 2 2 subscript superscript 𝜇 2 𝑔 𝑒 𝑥 2\mu^{2}_{g,ex} ), with μ g , e ​ x subscript 𝜇 𝑔 𝑒 𝑥 \mu_{g,ex} being the transition dipole moment, is also reported. All dipole data are given in Debye.

Molecule	QM/PCM				QM/FQ				QM/FQF$\mu$
	${\mu }_g$	${\mu }_{ex}$	$\Delta {\mu }^2$	$2{\mu }_{g,ex}^2$	${\mu }_g$	${\mu }_{ex}$	$\Delta {\mu }^2$	$2{\mu }_{g,ex}^2$	${\mu }_g$	${\mu }_{ex}$	$\Delta {\mu }^2$	$2{\mu }_{g,ex}^2$
pNA	10.5	17.6	50.4	72.5	10.3	15.0	22.1	38.1	12.9	19.0	37.2	68.0
pyridine	3.1	0.2	8.4	0.5	3.6	0.9	7.3	0.4	4.4	1.6	7.8	0.5
pyrimidine	3.3	0.8	6.2	0.7	3.6	1.9	2.9	0.6	4.4	2.7	2.9	0.7

Equations

$$\mathcal{E}(\mathbf{q},\mathbf{\mu })$$

$$+\sum _i\sum _{j\ne i}{q}_i{\text{T}}_{ij}^{q\mu }{\mathbf{\mu }}_j+\frac{1}{2}\sum _i\sum _{j\ne i}{\mathbf{\mu }}_i^{†}{\mathbf{T}}_{ij}^{\mu \mu }{\mathbf{\mu }}_j-\frac{1}{2}\sum _i{\mathbf{\mu }}_i^{†}{\alpha }_i^{-1}{\mathbf{\mu }}_i$$

$$\mathcal{E}(\mathbf{q},\mathbf{\mu })$$

$$+\sum _i\sum _j{q}_i{\text{T}}_{ij}^{q\mu }{\mathbf{\mu }}_j^{†}+\sum _i{q}_i{\chi }_i=$$

$$=\frac{1}{2}{\mathbf{q}}^{†}{\mathbf{T}}^{qq}\mathbf{q}+\frac{1}{2}{\mathbf{\mu }}^{†}{\mathbf{T}}^{\mu \mu }\mathbf{\mu }+{\mathbf{q}}^{†}{\mathbf{T}}^{q\mu }\mathbf{\mu }+{\mathbf{\chi }}^{†}\mathbf{q}$$

$$\mathcal{E}\left(\mathbf{q},\mathbf{\mu },\mathbf{\lambda }\right)$$

$$=\frac{1}{2}\sum _{i\alpha }\sum _{j\beta }{q}_{i\alpha }{\text{T}}_{i\alpha ,j\beta }^{qq}{q}_{j\beta }+\frac{1}{2}\sum _i\sum _j{\mathbf{\mu }}_{i\alpha }^{†}{\mathbf{T}}_{i\alpha ,j\beta }^{\mu \mu }{\mathbf{\mu }}_{j\beta }+$$

$$+\sum _i\sum _j{q}_{i\alpha }{\text{T}}_{i\alpha ,j\beta }^{q\mu }{\mathbf{\mu }}_{j\beta }^{†}+\sum _{i\alpha }{q}_{i\alpha }{\chi }_{i\alpha }+$$

$$+\sum _{\alpha }{\lambda }_{\alpha }\left[\sum _i{q}_{\alpha i}-Q_{\alpha }\right]=$$

$$=\frac{1}{2}{\mathbf{q}}^{†}{\mathbf{T}}^{qq}\mathbf{q}+\frac{1}{2}{\mathbf{\mu }}^{†}{\mathbf{T}}^{\mu \mu }\mathbf{\mu }+{\mathbf{q}}^{†}{\mathbf{T}}^{q\mu }\mathbf{\mu }+{\mathbf{\chi }}^{†}\mathbf{q}+{\mathbf{\lambda }}^{†}\text{q}$$

$$\left(\begin{array}[]{ccc}\mathbf{T}^{qq}&\mathbf{1}_{\bm{\lambda}}&\mathbf{T}^{q\mu}\\ \mathbf{1}^{\dagger}_{\bm{\lambda}}&\mathbf{0}&\mathbf{0}\\ -\mathbf{T}^{q\mu^{\dagger}}&\mathbf{0}&\mathbf{T}^{\mu\mu}\end{array}\right)\left({\begin{array}[]{c}\mathbf{q}\\ \bm{\lambda}\\ \bm{\mu}\end{array}}\right)=\left(\begin{array}[]{c}-\bm{\chi}\\ \mathbf{Q}\\ \mathbf{0}\end{array}\right)\Rightarrow\mathbf{D}\mathbf{L}_{\lambda}=-\mathbf{C}_{Q}$$

$$\mathcal{E}(\mathbf{P},\mathbf{q},\mathbf{\mu },\mathbf{\lambda })$$

$$+{\mathbf{q}}^{†}{\mathbf{T}}^{q\mu }\mathbf{\mu }+{\mathbf{\chi }}^{†}\mathbf{q}+{\mathbf{\lambda }}^{†}\text{q}+{\mathbf{q}}^{†}\mathbf{V}(\mathbf{P})-{\mathbf{\mu }}^{†}\mathbf{E}(\mathbf{P})$$

$${\tilde{F}}_{\mu \nu }=\frac{\partial \mathcal{E}}{\partial P_{\mu \nu }}=h_{\mu \nu }+G_{\mu \nu }(\mathbf{P})+{\mathbf{V}}_{\mu \nu }^{†}\mathbf{q}-{\mathbf{E}}_{\mu \nu }^{†}\mathbf{\mu }$$

$$\displaystyle\left(\begin{array}[]{ccc}\mathbf{T}^{qq}&\mathbf{1}_{\bm{\lambda}}&\mathbf{T}^{q\mu}\\ \mathbf{1}^{\dagger}_{\bm{\lambda}}&\mathbf{0}&\mathbf{0}\\ -\mathbf{T}^{q\mu^{\dagger}}&\mathbf{0}&\mathbf{T}^{\mu\mu}\end{array}\right)\left({\begin{array}[]{c}\mathbf{q}\\ \bm{\lambda}\\ \bm{\mu}\end{array}}\right)$$

$$\mathbf{D}{\mathbf{L}}_{\lambda }$$

$$\left(\begin{array}[]{cc}\tilde{\mathbf{A}}&\tilde{\mathbf{B}}\\ \tilde{\mathbf{B}}^{*}&\tilde{\mathbf{A}}^{*}\\ \end{array}\right)\left(\begin{array}[]{c}\mathbf{X}\\ \mathbf{Y}\\ \end{array}\right)=\omega\left(\begin{array}[]{cc}\bm{1}&\bm{0}\\ \bm{0}&-\bm{1}\\ \end{array}\right)\left(\begin{array}[]{c}\mathbf{X}\\ \mathbf{Y}\\ \end{array}\right)$$

$${\tilde{A}}_{ai,bj}=({ϵ}_a-{ϵ}_i){\delta }_{ab}{\delta }_{ij}+(ai|bj)-c_x(ab|ij)+c_l{f}_{ai,bj}^{xc}+C_{ai,bj}^{pol}$$

$${\tilde{B}}_{ai,bj}=(ai|bj)-c_x(aj|ib)+C_{ai,bj}^{pol}$$

$$C_{ai,bj}^{FQ}=\sum _P^{N_q}\left({\int }_{{\mathbb{R}}^3}{\varphi }_a(\mathbf{r})\frac{1}{\left|\mathbf{r}-{\mathbf{r}}_p\right|}{\varphi }_i(\mathbf{r})\mathrm{d}\mathbf{r}\right)\cdot q_P^T({\varphi }_b,{\varphi }_i)$$

$$C_{ai,bj}^{FQF\mu }$$

$$-\sum _P^{N_{\mu }}\left({\int }_{{\mathbb{R}}^3}{\varphi }_a(\mathbf{r})\frac{(\mathbf{r}-{\mathbf{r}}_p)}{{\left|\mathbf{r}-{\mathbf{r}}_p\right|}^3}{\varphi }_i(\mathbf{r})\mathrm{d}\mathbf{r}\right)\cdot {\mathbf{\mu }}_P^T({\varphi }_b,{\varphi }_i)$$

$$\displaystyle\left(\begin{array}[]{ccc}\mathbf{T}^{qq}&\mathbf{1}_{\bm{\lambda}}&\mathbf{T}^{q\mu}\\ \mathbf{1}^{\dagger}_{\bm{\lambda}}&\mathbf{0}&\mathbf{0}\\ -\mathbf{T}^{q\mu^{\dagger}}&\mathbf{0}&\mathbf{T}^{\mu\mu}\end{array}\right)\left({\begin{array}[]{c}\mathbf{q}\\ \bm{\lambda}\\ \bm{\mu}\end{array}}\right)$$

$$\mathbf{D}{\mathbf{L}}_{\lambda }^T$$

$$V({\mathbf{P}}_K^T)$$

$$\text{E}({\mathbf{P}}_K^T)$$

$${\omega }_{K,FQ}^{cLR}={\omega }_K^0+\frac{1}{2}\sum _P^{Nq}{q}_p({\mathbf{r}}_p;{\mathbf{P}}_K^{\Delta })V({\mathbf{r}}_p,{\mathbf{P}}_K^{\Delta })$$

$${\omega }_{K,FQF\mu }^{cLR}={\omega }_K^0$$

$$-\frac{1}{2}\sum _P^{N_{\mu }}{\mathbf{\mu }}_p({\mathbf{r}}_p;{\mathbf{P}}_K^{\Delta })\mathbf{E}({\mathbf{r}}_p,{\mathbf{P}}_K^{\Delta })$$

$${\mathbf{P}}_K^{\Delta }={\mathbf{P}}_K^T+{\mathbf{Z}}_K$$

$${\epsilon }_i(\tilde{\nu })=\frac{\sqrt{\pi }{e}^2{N}_A}{1000\ln (10)c^2{m}_e}\frac{f_i}{\sigma }\exp \left[-{\left(\frac{\tilde{\nu }-{\tilde{\nu }}_i}{\sigma }\right)}^2\right]$$

Full text

††preprint: AIP/123-QED

I Introduction

Excited-state phenomena play a crucial role in many application fields, as for instance photocatalysis, optical information storage and solar cells. In the past decades, theoretical modeling of excited-state properties of molecules in the gas-phase has become a widespread strategy of investigation,Helgaker et al. (2012) giving precious information on, for instance, the nature of the electronic excitation,Le Bahers, Adamo, and Ciofini (2011); Guido et al. (2013); Savarese et al. (2017) nuclei-electron coupling effectsBaiardi, Bloino, and Barone (2016); Bloino, Baiardi, and Biczysko (2016); Improta et al. (2009) and excited state electron dynamics.Barbatti (2011); Barbatti et al. (2010); González, Escudero, and Serrano-Andrés (2012)

However, for electronic phenomena taking place in the condensed phase,Cannelli et al. (2017); Carlotti et al. (2018); Barone et al. (2012); Reichardt (1994); Mennucci, Cammi, and Tomasi (1998); Jacquemin, Mennucci, and Adamo (2011); Corni et al. (2005); Cupellini et al. (2019); Mennucci and Corni (2019) the interplay between the molecule and its environment can substantially alter the electronic response to external electromagnetic fields. Therefore, any accurate modeling of excited states of solvated systems asks for reliable theoretical approaches to include the effects of the environment at all levels of the excitation phenomenon.

Most of the currently available approaches to describe the effects of the external environment on molecular properties belong to the class of the so-called focused models;Warshel and Levitt (1976); Warshel and Karplus (1972); Miertuš, Scrocco, and Tomasi (1981); Tomasi and Persico (1994); Orozco and Luque (2000); Tomasi, Mennucci, and Cammi (2005) the attention is focused on the molecule and the environment is treated a lower level of sophistication as it modifies, but not determines, the molecular response to the external radiation. In order to keep the atomistic description of the environment, thus substantially overcoming well-known and amply used continuum solvent descriptions,Tomasi, Mennucci, and Cammi (2005); Tomasi and Persico (1994); Mennucci and Cammi (2007); Mennucci (2010, 2013); Lipparini and Mennucci (2016) multiscale QM/Molecular Mechanics (MM) approaches have been developed.Senn and Thiel (2009); Lin and Truhlar (2007) In such models, the molecule (solute) is treated at the QM level, whereas the environment (solvent) is modelled by means of classical MM force fields (FF). The interaction between the QM and MM portions is usually described in terms of electrostatic forces, although approaches to include non-electrostatic QM/MM interactions have been proposed.Jensen and Gordon (1996, 1998); Ben-Nun and Martínez (1998); Giovannini, Lafiosca, and Cappelli (2017); Curutchet et al. (2018); Giovannini et al. (2019a); Nåbo et al. (2016); Reinholdt, Kongsted, and Olsen (2017) In order to fully capure the physics of solute-solvent interactions, mutual polarization effects between the QM and MM moieties need to be accounted for. Therefore, several polarizable QM/MM approaches have been proposed, based on distributed multipoles,Day et al. (1996); Kairys and Jensen (2000); Mao et al. (2016); Loco et al. (2016); Loco and Cupellini (2018) induced dipoles,Thole (1981); Curutchet et al. (2009); Olsen and Kongsted (2011); Steindal et al. (2011); Jurinovich, Curutchet, and Mennucci (2014) Drude oscillatorsBoulanger and Thiel (2012), Fluctuating Charges (FQ)Rick, Stuart, and Berne (1994); Rick and Berne (1996); Cappelli (2016) and the recently developed approach based on both Fluctuating Charges and Fluctuating Dipoles (FQF$\mu$).Giovannini et al. (2019b, 2019)

The latter can be seen as a refinement of QM/FQ, previously developed by some of us, Giovannini, Ambrosetti, and Cappelli (2018); Cappelli (2016); Lipparini, Cappelli, and Barone (2012); Lipparini et al. (2012); Lipparini, Cappelli, and Barone (2013); Giovannini et al. (2017); Giovannini, Olszowka, and Cappelli (2016); Giovannini et al. (2018); Egidi et al. (2019); Puglisi et al. (2019); Di Remigio et al. (2019) where the MM portion is described by means of electric charges, which vary as a function of differences in MM atomic electronegativities and as a response to the electric potential generated by the QM density. Therefore, in QM/FQ only monopoles, i.e. zeroth order of the electrostatic Taylor expansion, are taken into consideration. This means that the intrinsic anisotropy and out-of-plane contributions of molecule-environment interactions are not explicitly taken into account. To overcome this problem, the electrostatic description given by FQ has been refined by including fluctuating dipoles (F$\mu$), thus resulting in the QM/FQF$\mu$ approach. Notice that similar approaches have been developed in other contexts,Stern et al. (1999); Naserifar et al. (2017); Oppenheim, Naserifar, and Goddard III (2018); Mayer (2007); Jensen and Jensen (2009); Rinkevicius et al. (2014) however they are not specifically intended to model molecular systems/properties in solution.

The potentialities of QM/FQF$\mu$ at predicting electrostatic interaction energies of organic molecules in solution have been previously highlighted.Giovannini et al. (2019b) In this work, this model is further extended to the calculation of vertical excitation energies, by resorting to linear response (LR) theory for Self Consistent Field (SCF) methods. Excitation energies are determined by solving the so-called Casida equations,Casida (1995) for the QM portion, appropriately modified so to include extra terms due to QM/MM mutual polarization effects.Lipparini, Cappelli, and Barone (2012); Loco et al. (2016); Rinkevicius et al. (2014); Steindal et al. (2011) The LR approach may reliably describe electronic transitions involving excited states associated with large transition dipole moments. However, it lacks any accounting of the relaxation of the classical portion as a response to charge equilibration following the excitation of the QM density. For this reason, several state specific (SS) approaches have been proposed in the literature both in case of continuum modelsCaricato et al. (2006); Marenich et al. (2011); Marenich, Cramer, and Truhlar (2015); Budzak et al. (2014); Improta et al. (2006); Corni et al. (2005); Caricato (2014); Guido et al. (2018); Duchemin et al. (2018); Schröder and Schwabe (2018); Guido and Caprasecca (2019); Guido et al. (2017, 2015) and polarizable QM/MM approaches.Zeng and Liang (2015); Loco et al. (2016) In this work, the perturbative SS corrected linear response (cLR) schemeCaricato et al. (2006) originally developed for the Polarizable Continuum Model (PCM) is exploited and extended to both QM/FQ and QM/FQF$\mu$ transition energies for the first time. In this way, the classical MM portion is adjusted to the relaxed density matrix of a given QM excited state. Notice that a similar approach has been applied to polarizable QM/AMOEBA calculationsLoco et al. (2016).

To show both the potentialities of QM/FQF$\mu$ at describing excitation energies and the differences arising from resorting to LR and cLR regimes, the method is applied to the calculation of vacuo-to-water solvatochromic shifts of para-nitroaniline (pNA), pyridine and pyrimidine. From such an analysis, it will be possible to (i)highlight differences in the description of the excitation phenomenon moving from continuum PCM to QM/FQ and QM/FQF$\mu$ fully atomistic approaches; (ii) directly quantify the separate role fluctuating dipoles and fluctuating charges, and their coupling; (iii) evaluate effects arising from the specification of LR or cLR regimes.

The manuscript is organized as follows. In the next section, the FQF$\mu$ force field is presented and its coupling with a QM description at the SCF level (QM/FQF$\mu$) is detailed. LR and cLR working equations are then presented and discussed. After a brief section discussing on the computational protocol which is adopted, numerical results are presented. Some drawn conclusions and a discussion on the future perspectives of the approach end the manuscript.

II Theory

II.1 QM/FQF$\mu$ Approach

In the FQF$\mu$ force field each MM atom is endowed with both a charge ($q$) and an atomic dipole ($\bm{\mu}$), that are not fixed but can vary according to the external electric potential/ electric field.

The total energy $\mathcal{E}$ associated with a gaussian distribution of charges and dipoles reads: Mayer (2007)

[TABLE]

where $\chi$ is the atomic electronegativity, $\eta$ the chemical hardness and $\alpha$ the atomic polarizability. $\text{T}_{ij}^{qq}$, $\mathbf{T}_{ij}^{q\mu}$ and $\mathbf{T}_{ij}^{\mu\mu}$ are the charge-charge, charge-dipole and dipole-dipole interaction kernels, respectively. Their functional form can be found in Refs.69; 54.

In order to collect all the charge quadratic terms, the diagonal elements of $\textbf{T}^{qq}$ and $\textbf{T}^{\mu\mu}$ can be imposed to be related to atomic chemical hardness $\eta$ and atomic polarizability $\alpha$, respectively.Giovannini et al. (2019b) Therefore, the parameters entering the definition of FQF$\mu$ are limited to electronegativity, chemical hardness and polarizability of each atom type. As a result, Eq. 1 can be formally rewritten as:

[TABLE]

where a vector notation has been adopted. Similarly to QM/FQ, the equilibrium condition of Eq. 2 is reached when the Electronegativity Equalization Principle (EEP) is satisfied, i.e. when each atom has the same electronegativity.Sanderson (1951) If each MM moiety is constrained to assume a fixed, total charge Qα, Eq. 2 can be written by exploiting a set of Lagrangian multipliers ($\lambda_{\alpha}$), whose number is equal to the total number of moieties in the MM portion:

[TABLE]

where $\alpha$ and $\beta$ run over MM moieties. By minimizing $\mathcal{E}$ with respect all variables, the following linear system is defined:Giovannini et al. (2019b)

[TABLE]

where $\textbf{1}_{\bm{\lambda}}$ is a rectangular matrix, which accounts for the Lagrangians. $\mathbf{C}_{Q}$ is a vector containing atomic electronegativities and total charge constraints, whereas $\mathbf{L}_{\lambda}$ is a vector containing charges, dipoles and Lagrange multipliers.

FQF$\mu$ FF can be effectively coupled to a QM SCF description, within a QM/MM framework:

[TABLE]

where $\mathbf{h}$ and $\mathbf{G}$ are the usual one- and two-electron matrices, and P is the density matrix. $\mathbf{q}^{\dagger}\mathbf{V}(\mathbf{P})$ and $\bm{\mu}^{\dagger}\mathbf{E}(\mathbf{P})$ are the electrostatic couplings between the QM density and the FQs and F$\mu$s, respectively. The effective Fock matrix, i.e. the derivative of the energy with respect to the density matrix, reads:

[TABLE]

Charges and dipoles are obtained by imposing the global functional to be stationary with respect to charges, dipoles and Lagrangian multipliers:

[TABLE]

Notice that, with respect to Eq. 4, a new term, $\mathbf{R}(\mathbf{P}$), appears, which collects QM polarization sources. $\mathbf{L}_{\lambda}$ is the vector containing charges, dipoles and Lagrangian multipliers.

II.2 QM/FQF$\mu$ Electronic Transition Energies

II.2.1 Linear Response Regime

LR for polarizable QM/MM approaches has been already discussed in the previous literature.Lipparini, Cappelli, and Barone (2012); Rinkevicius et al. (2014); Steindal et al. (2011) Electronic excitation energies $\omega$ for SCF Hamiltonians can be obtained by solving Casida’s equations:Casida (1995)

[TABLE]

$\tilde{A}$ and $\tilde{B}$ matrices are defined as:

[TABLE]

[TABLE]

where $(pq|rs)$ are two electron integrals, and $\varepsilon$ are molecular orbital (MO) energies. $c_{x}$ and $c_{l}$ are coefficients of which the definition depends on the SCF level adopted ($c_{x}$ =1, $c_{l}$ = 0 for HF wavefunctions, $c_{x}$ =0, $c_{l}$ = 1 for DFT).

In case of polarizable QM/MM approaches, $\tilde{A}$ and $\tilde{B}$ in Eqs. 22 and 23 account for an extra term, $C^{pol}$, which is specified according to a particular method. In particular, for QM/FQ it reads:

[TABLE]

where, $q^{T}$ are the perturbed fluctuating charges adjusted to the transition density $\mathbf{P}^{T}_{K}=\mathbf{X}_{K}+\mathbf{Y}_{K}$.Lipparini, Cappelli, and Barone (2012) For QM/FQF$\mu$, an extra term appears n the definition of $C^{pol}$, accounting for the presence of fluctuating dipoles:

[TABLE]

Perturbed charges ($q^{T}$) and perturbed dipoles ($\bm{\mu}^{T}$) are calculated by solving the following system of equations:

[TABLE]

where:

[TABLE]

The right hand side of Eq. 35 contains both the electric potential and field due to the perturbed density matrix $\mathbf{P}^{T}_{K}$. Notably, the constant term in Eq. 19 due to atomic electronegativities vanishes. Therefore, polarization arises from the QM density, and not from differences in electronegativity, that similarly to the QM/FQ approach.Lipparini, Cappelli, and Barone (2012); Giovannini et al. (2019)

II.2.2 Corrected Linear Response Regime

LR and SS approaches provide the same result in case of gas-phase systems.Casida (1995) However, when polarizable QM/classical methods are exploited, the two approaches differ.Corni et al. (2005) This is essentially due to the non-linear term introduced in the Hamiltonian (see e.g. Eq. 5). Different approximations to solve this problem have been proposed.Caricato et al. (2006); Improta et al. (2006); Marenich et al. (2011)In this work, the SS cLR approach, firstly proposed in for QM/PCM Hamiltoniands and successfully extended to the polarizable QM/AMOEBA methodLoco et al. (2016), is exploited. To the best of our knowledge, this the first time that the extension of cLR to QM/FQ and QM/FQF$\mu$-like methods is proposed.

The common way to discuss the physical differences between LR and cLR regimes in condensed phase is to resort to a two-step protocol, in which the first step is the same in both approaches and consists in the electronic excitation of the solute to the $K$-th excited state by keeping frozen the MM response to the solute ground state. From a computational point of view, such a picture can be achieved by imposing the extra $C^{pol}$ term in Eqs. 22 and 23 to be equal to zero. This means that MM contributions are only included in the specification of MOs through Eq. 6. The excitation energy obtained in this approximation is refereed to as $\omega_{0}^{K}$. In the second step of the process, MM electronic degrees of freedom adjust to the $K$-th excited state density. Due to the fact that in the LR approach the MM variables are determined as a response of the whole transition densities, LR is unable to catch the energy differences associated to the relaxation of the QM electron density, whereas it only accounts for a correction due to dynamic solute-environment interactions. QM/FQ and QM/FQF$\mu$ cLR excitation energies from ground to $K$-th excited state can be written as:

[TABLE]

[TABLE]

$\mathbf{P}^{\Delta}_{K}$ is the so-called relaxed-density matrix, computed through the so-called Z-vector approach, as:

[TABLE]

where $\mathbf{P}_{K}^{T}$ is the unrelaxed density matrix adopted in the LR approach, whereas $\mathbf{Z}_{K}$ is the Z-vector contribution which accounts for orbital relaxations. Fluctuating Charges and Fluctuating Dipole in Eq. II.2.2 are calculated as in Eq. 35, by substituting $\mathbf{P}^{T}_{K}\rightarrow\mathbf{P}^{\Delta}_{K}$.

III Computational Details

The following computational protocol was exploited:Giovannini et al. (2018, 2019); Giovannini, Ambrosetti, and Cappelli (2018); Giovannini et al. (2019a)

Definition of the system: The solute (pNA, pyridine and pyrimidine) was put at the center of a cubic box containing a number of water molecules large enough to account for solute-solvent interactions (vide infra). pNA, pyridine and pyrimidine geometries were optimized at the B3LYP/aug-cc-pVDZ level of theory, and RESP charges were calculated by including aqueous solvent effects by means of the PCM approach.Tomasi, Mennucci, and Cammi (2005). 2. 2.

Classical Molecular Dynamics (MD) runs: For each system, MD simulations were performed by imposing periodic boundary conditions (PBC) on the cubic box. Each MD run was preceded by an equilibration step. From each MD run, a set of 100 uncorrelated snapshots was extracted. 3. 3.

Definition of the the two-layer scheme: For each snapshot, a solute-centered sphere was cut. The radius of the sphere was chosen so to describe all specific water-solute interactions (vide infra). 4. 4.

QM/MM calculations: QM/FQ and QM/FQF$\mu$ vertical excitation energies were calculated on the spherical frames obtained at the previous step. The results obtained for each frame were convoluted with a gaussian band-shape and averaged to produce the final spectrum:

[TABLE]

where, $\varepsilon$ is the molar absorptivity in L mol*-1cm-1*, $\tilde{\nu}_{i}$ is the excitation energy of the $i-th$ state, $f_{i}$ is its oscillator strength, whereas $\sigma$ is the standard deviation of the gaussian convolution. $e$, $N_{A}$, $c$ and $m_{e}$ are the electron charge, the Avogadro number, the speed of light, and the electron mass, respectively.

MD simulations were performed using GROMACSAbrahama et al. (2015), by exploiting GAFF Oostenbrink et al. (2004) to describe both intra- and inter-molecular interactions. RESP charges were used to account for electrostatic interactions, and TIP3P FF was used to describe water molecules Mark and Nilsson (2001). A single molecule of different solutes was dissolved in a cubic box containing 5781, 5342 and 5324 water molecules in case of pNA, pyridine and pyrimidine, respectively. Pyridine and pyrimidine structures were kept fixed during all the steps of the MD run because they exhibit a single conformer in aqueous solution. PNA geometry was kept fixed during the equilibration step only, whereas it was free to move in the production run. The different solutes were brought to 0 K with the steepest descent minimization procedure and then heated to 298.15 K in an NVT ensemble using the velocity-rescalingBussi, Donadio, and Parrinello (2007) method with an integration time step of 1.0 fs and a coupling constant of 0.1 ps for 100 ps. NPT simulations (using the Parrinello-Rahman barostat and a coupling constant of 1.0 ps) for 1 ns were performed to obtain a uniform distribution of molecules in the box. 5 ns production runs in the NVT ensemble were finally carried out, fixing the fastest internal degrees of freedom, i.e. hydrogen atoms, by means of the LINCS algorithm ($\delta$t=2.0 fs)Hess et al. (1997). Electrostatic interactions were treated by using particle-mesh Ewald (PME) Darden, York, and Pedersen (1993) method with a grid spacing of 1.2 Å and a spline interpolation of order 4. Intramolecular interactions between atom pairs separated up to three bonds were excluded. A snapshot every 50 ps was extracted in order to obtain a total of 100 uncorrelated snapshots for each system.

For each snapshot a solute-centered sphere with a radius of 15 Å was cut (containing approximately 400 water molecules). On each droplet QM/FQ and QM/FQF$\mu$ vertical excitation energies were calculated by exploiting both LR and cLR regimes. In case of QM/FQ, water molecules were modeled by exploiting polarizable FQ SPC parametrization proposed by Carnimeo, Cappelli, and Barone (2015), whereas in case of QM/FQF$\mu$ calculations the parametrization proposed by Giovannini et al. (2019b) was adopted. For the sake of comparison, additional QM calculations on the chromophores treated at the PCM level were performed. PNA was described at the CAM-B3LYP/aug-cc-pVDZ level of theory in agreement with previous studies.Egidi et al. (2014) Pyridine and Pyrimidine were instead treated at the M06/6-311+G(2df,2p) level, according to Ref. 75. All computed absorption spectra were convoluted with a gaussian band-shape with a Full Width at Half Maximum (FWHM) of 0.3 eV.

All QM/FQ and QM/FQF$\mu$ LR and cLR calculations were performed by using a locally modified version of the Gaussian 16 packageFrisch et al. (2016).

IV Numerical Results

In this section, computed UV-Vis spectra for aqueous solutions of pNA, pyridine and pyrimidineSok et al. (2011); Kosenkov and Slipchenko (2010); Eriksen et al. (2013); Sneskov et al. (2011); Olsen, Aidas, and Kongsted (2010); Eriksen et al. (2012); Frutos-Puerto, Aguilar, and Fdez. Galv’an (2013); DeFusco et al. (2011); Mennucci (2002); Pagliai et al. (2017); Biczysko et al. (2012); Cossi and Barone (2001); Marenich et al. (2011); Mason (1959); Millefiori et al. (1977); Cai and Reimers (2000); Schreiber et al. (2008); da Silva et al. (2010); Prabhumirashi and Kunte (1986); Kovalenko et al. (2000); de Almeida et al. (2001); Silva-Junior et al. (2010) (see Fig. 1 for their structures and atom labelling) are discussed. First, excitation energies of the three isolated molecules are discussed and compared to experimental values. Then, QM/FQF$\mu$, QM/FQ and QM/PCM LR and cLR solvatochromic shifts are presented and compared to experimental findings.

IV.1 Gas-Phase Vertical Excitation Energies

Computed gas phase vertical excitation energies for the three systems are reported in Table 1, together with experimental values, taken from the literature. Calculated ground and excited dipole moments are also reported.

The molecular orbitals (MOs) involved in the transitions are graphically depicted in Figure 2. The investigated excitation is a pure $\pi\rightarrow\pi^{*}$ transition in case of pNA, whereas in case of both pyridine and pyrimidine we focus our attention to the dark $n\rightarrow\pi^{*}$ transition. pNA $\pi\rightarrow\pi^{*}$ transition is often considered to have a Charge Transfer (CT) character, due to the fact that it involves a huge variation of the dipole moment between the ground and excited state; this is confirmed by the computed values, see Table 1. The same also occurs of pyridine and pyrimidine, which differently from pNA show a decrease of the electric dipole moment in their excited states with respect to their ground state.

The comparison between computed and experimental vertical excitation energies in Table 1, confirms that our choice of DFT functionals/basis sets give computed data that are nicely in agreement with experiments for all systems.

IV.2 MD Analysis

In this section, MD runs for pNA, pyridine and pyrimidine in aqueous solution are analyzed by using the TRAVIS package, Brehm and Kirchner (2011) so to extract hydration patterns and analyze the role of hydrogen bonding (HB) interactions, which are crucial to rationalize the trends of vertical excitation energies (vide infra). In particular, the radial distribution function $g(r)$, together with the running coordination number (RCN) between the solute atoms highlighted in Fig. 1 and solvent atoms (HW and OW) are studied. the RCN for a specific site is defined as the integral of the $g(r)$ as a function of the distance. The coordination number of a specific hydration shell is the RCN at the corresponding $g(r)$ minima. In addition, in case of pNA, the dihedral distribution function (ddf) is also discussed.

$g(r)$ and RCN for the three studied molecules are reported in Figs. 3 and 4. In case of both pyridine and pyrimidine, only the nitrogen atom(s) can be involved in HB interactions with water molecules. For this reason, in Fig. 3, panels b and c, the $g(r)$ for N1$\cdots$OW (and N2$\cdots$OW in case of pyrimidine) are reported. As it can be evinced, the $g(r)$ of both pyridine and pyrimidine present a sharp peak related to the first solvation shell at about 2 Å. The coordination numbers associated with the first solvation shell are approximately 2 and 1, for pyridine and pyrimidine, respectively. This indicates that two water molecules are involved in HB interactions with the nitrogen atom(s) of this two solutes.

In case of pNA (see Fig. 3, panel a), H1$\cdots$OW and O1$\cdots$HW $g(r)$ are both characterized by a sharp peak at about 2 Å, whereas N1 and N2 are not involved in HB interactions with the solvent. The coordination numbers of the first shells are around 1, showing that both the push and pull groups are involved in HB interactions with one water molecule, on average. Notice that pNA is symmetric, thus two water molecules undergo HB interactions with both push and pull groups of the solute.

To end the discussion, the ddf of $\gamma$, $\delta$ and $\tau$ pNA dihedral angles are reported in Fig. 5 (see Fig. 1 for atoms involved in the definition of the dihedral angles). As it can be seen, all the considered dihedral angles fluctuate around 0 degrees (or 180 degrees in case of $\tau$), thus indicating that pNA oscillates around its planar configuration during the MD simulation. Specifically, out-of-plane motions of the benzene ring were monitored by analyzing $\delta$ dihedral angle, which only marginally varies during all MD simulations, thus indicating an almost rigid structure of this ring. In addition, the analysis of $\gamma$ and $\theta$ dihedral angles shows that the oxygen atoms of the nitro group can pass the rotational barrier (e.g. they can exchange their positions), whereas the same does not apply to the hydrogen atoms of the amine group. Indeed, they can oscillate in the molecule’s plane, but they do not exchange their positions during the MD runs.

IV.3 Vacuo-to-Water Solvatochromism

QM/FQ and QM/FQF$\mu$ vertical excitation energies of pNA (CAM-B3LYP/aug-cc-pVDZ), pyridine and pyrimidine (M06/6-311+G(2df,2p) in aqueous solution were computed on 100 uncorrelated snapshots extracted from MD runs.

In Figure 6, pNA, pyridine and pyrimidine cLR absorption intensities as a function of the excitation energy, obtained by exploiting QM/FQ (top) and QM/FQF$\mu$ (bottom), for all 100 snapshots extracted from the MD simulations are plotted. Notice that stick spectra are obtained by plotting the raw data extracted from QM/MM calculations performed on each snapshot. $\omega_{0}$ and LR stick spectra are reported in Figs. S1 and S2 in the Supplementary Material (SM).

In all cases, the inspection of Fig. 6 shows that QM/FQF$\mu$ data (sticks) present a larger variability with respect to QM/FQ, which results in a greater broadening of the transition band, and in a different description of vertical excitation energies. To quantify such a variability, the standard deviation of both energies and intensities was computed (see Table S1 given as SM). As expected, QM/FQF$\mu$ energy standard deviations are always larger than QM/FQ ones. Such findings are due to the inclusion of fluctuating dipoles in the MM portion, which increase the molecular dipole moment of MM water molecules. As a consequence, QM/FQF$\mu$ convoluted spectra in Fig. 6 show a larger inhomogeneous broadening with respect to QM/FQ. To further investigate the origin of QM/FQF$\mu$ larger inhomogeneous broadening, in Fig. 7 pNA QM/FQ (top) and QM/FQF$\mu$ (bottom) excitation energies as a function of $\gamma$ and $\tau$ dihedral angles (see Fig. 1) are reported. Such a plot permits to understand whether the observed trend in broadening is related to the solute’s geometry fluctuations or to solvent fluctuations, i.e. to the different description of the electrostatic interactions given by QM/FQ and QM/FQF$\mu$. As it is clearly shown by Fig. 7, the two approaches give almost identical excitation energies as a function of the dihedral angles. Therefore, it is confirmed that the different broadening is only due to a different description of the QM/MM electrostatic coupling.

PNA, pyridine and pyrimidine $\omega_{0}$, LR and cLR vertical excitation energies in aqueous solution, as calculated by exploiting QM/PCM, QM/FQ and QM/FQF$\mu$, and the corresponding vacuo-to-water solvatochromic shifts are reported in Tab. 2. The associated spectra are graphically depicted in Fig. 8, where gas-phase excitation energies are also given as stick bars. By moving from a continuum description (QM/PCM) to an atomistic approach (QM/FQ and QM/FQF$\mu$) a gradual increase of the solvatochromic shift is reported, independently from which vertical excitation energy ($\omega_{0}$, LR or cLR) is taken in consideration (see Tab. 2 and Fig. 8). In particular, in case of pNA, QM/FQ redshifts with respect to QM/PCM of 0.11 eV, whereas QM/FQF$\mu$ of almost 0.50 eV, thus doubling the computed solvatochromism with respect to continuum QM/PCM values. In case of pyridine and pyrimidine, a shift in the opposite direction with respect to pNA is reported by all the solvation models. In particular, QM/FQ and QM/FQF$\mu$ predict an increase of vacuo-to-water solvatochromism of 0.18/0.21 and 0.33/0.35 eV with respect to QM/PCM, respectively. Notice that such an increase corresponds to a relative difference with respect to QM/PCM of almost 300 and 500 %, in case of QM/FQ and QM/FQF$\mu$, respectively. This huge blueshift can be related to the dynamical atomistic description of the solvent water molecules given by QM/MM approaches coupled to MD simulations. The difference between QM/FQ and QM/FQF$\mu$ results reported by all the selected molecular systems, shows instead the relevance of the inclusion of fluctuating dipoles in addition to fluctuating charges.

By deepening the analysis of Tab. 2, the different pictures given by LR and cLR approaches emerge. First, let’s focus on pNA $\omega_{0}$, LR and cLR results. The continuum QM/PCM and atomistic QM/FQ and QM/FQF$\mu$ LR and cLR excitation energies are very similar, and shifted with respect to the frozen density approximation ($\omega_{0}$) of 0.11, 0.04 and 0.16 eV, respectively. The similarity between LR and cLR energies indicates that in this particular case, the relaxation of the solvent (Which is considered in the cLR approach), is as relevant as the dynamical solute-solvent interactions taken into account by the LR regime. This is not surprising, by considering that the studied excitation is a pure $\pi\rightarrow\pi^{*}$ transition, for which the LR regime is generally successfull.Guido et al. (2015) Furthermore, the major differences with respect to $\omega_{0}$ are reported for LR energies, which are therefore selected as “best” results in case of pNA (and reported in bold in Tab. 2).

Differently from pNA, the $n\rightarrow\pi^{*}$ transition is studied in case of pyridine and pyrimidine. For such a transition, which is associated to a large variation of the dipole moment together with a transition dipole moment close to zero (i.e. the transition is dark), cLR regime has been previously reported to be essential in order to correctly account for the relaxation of the solvent both in QM/PCMCaricato et al. (2006) and polarizable QM/MM calculations.Loco et al. (2016) This is confirmed by the data reported in Tab. 2, where $\omega_{0}$ and LR excitation energies are predicted to be almost the same by all investigated solvation approaches, whereas a shift in case of cLR is reported. In particular, QM/PCM cLR corrections shift $\omega_{0}$ (and LR) energies of about 0.05 and 0.03 eV, for pyridine and pyrimidine, respectively; this corresponds to 33% deviation in both cases. The magnitude of the QM/FQ shift reported for both molecules, and of QM/FQF$\mu$ shift for pyrimidine are very similar to QM/PCM values. This is not surprising, because it has been recently reported to also apply to QM/AMOEBA cLR absorption energies Loco et al. (2016). QM/FQF$\mu$ cLR correction is particularly relevant for pyridine in aqueous solution, for which a shift of 0.09 (20%) with respect to the frozen density approximation $\omega_{0}$ (and LR) is obtained.

Best estimates of QM/PCM, QM/FQ and QM/FQF$\mu$ vacuo-to-water solvatochromic shifts for the three systems are highlighted in Tab. 2. They are obtained by exploiting the LR regime for pNA and cLR for both pyridine and pyrimidine. In case of QM/FQ and QM/FQF$\mu$, standard errors are also given, which are calculated as $\sigma/\sqrt{N}$, where $\sigma$ is the standard deviation and $N$ is the number of the snapshots used to obtain the average values of the shifts. We note that the continuum QM/PCM is generally unable to recover the experimental vacuo-to-water solvatochromism, that reasonably due to the imporper account of specific solute-solvent interactions and to the "static" description of the solvation phenomenon.Giovannini et al. (2017); Loco et al. (2018, 2016) As a consequence, the error on pNA solvatochromism is of 50%, whereas in case of pyridine and pyrimidine QM/PCM only recovers 35% of the experimental shift. A much better agreement with the experimental values is achieved by exploiting the polarizable QM/FQ coupled to MD simulations; in fact, this approach is able to recover the experimental vacuo-to-water solvatochromism in case of pyridine and pyrimidine, whereas it fails at reproducing the value for pNA.pNA is correctly described by QM/FQF$\mu$, which, however, overestimates of about 0.10 eV the values for pyridine and pyrimidine. Remarkably, in case of QM/FQF$\mu$, cLR shifts $\omega_{0}$ and LR vertical excitation energies towards their experimental counterparts (see Tab. 2).

To investigate the reasons of the overestimation of vacuo-to-water solvatochromic shifts reported for QM/FQF$\mu$ (which, as stated above, is particularly evident for pyridine and pyrimidine) the charge transfer (CT) nature of the considered transitions was characterized by a quantitative index, denoted as $D_{\mathrm{CT}}$.Le Bahers, Adamo, and Ciofini (2011); Egidi et al. (2018) The barycenters of the positive and negative density distributions are calculated as the difference of the ground state (GS) and excited state (ES) densities. The CT length index ($D_{\mathrm{CT}}$ ) is defined as the distance between the two charge barycenters. $D_{\mathrm{CT}}$ indexes, as calculated by considering both the unrelaxed $\mathbf{P}^{T}_{K}$ and relaxed densities $\mathbf{P}^{\Delta}_{K}$, are reported in Tab. 3.

The data in Tab. 3, confirm that all the considered transitions are associated with an intramolecular CT. In particular, the largest $D_{\mathrm{CT}}$ values are reported in case of pNA (both in gas phase and in solution), independently from the specific solvation model. Notice that the CT nature of each transition is different (lower) if the calculation is performed by exploiting the relaxed density ($\mathbf{P}^{\Delta}_{K}$) instead of the unrelaxed one ($\mathbf{P}^{T}_{K}$), in line with what reported in previous studies.Maschietto et al. (2018) Such an observation supports the necessity of taking into account the relaxation of the solvent polarization, i.e. to adopt a cLR regime rather than the common LR approach for the solvated systems. However, in order to correctly decide the best regime for our systems, molecular dipole moments need to be analyzed. In particular, differences in cLR and LR vertical excitation energies can be related to a different description of the change in the dipole moment moving from GS to ES ($\Delta\mu=\left|{\mu_{ex}-\mu_{g}}\right|$) and the transition dipole moment ($\mu_{g,ex}$). In fact, it has been shown that if $\Delta\mu^{2}$ is smaller than $2\mu^{2}_{g,ex}$, the LR excitation energy is smaller than cLR one, and vice versa.Caricato et al. (2006) Therefore, the relevant values are reported in Table 4 for the three investigated solvation approaches.

We notice that moving from GS to ES the dipole moment decreases for both pyridine and pyrimidine, whereas it increases for pNA. Remarkably, this is related to the negative solvatochromism reported in case of pyridine and pyrimidine, and positive in case of pNA.Marenich, Cramer, and Truhlar (2015) Also, the results obtained in Tab. 2 are coherent with the fact that for pNA $2\mu^{2}_{g,ex}$ is larger than $\Delta\mu^{2}$, whereas the opposite applies to both pyridine and pyrimidine. As a consequence, LR excitation energies are smaller than cLR for pNA, whereas cLR ones are smaller than LR for the other two systems (see Table 2). Furthermore, the analysis based on dipole moments (Tab. 4) suggest that in case of pNA, QM/FQF$\mu$ outperforms QM/FQ because it predicts a larger difference (of almost 25%) between ES and GS dipole moments, thus resulting in a larger vacuo-to-solvent solvatochromism. The same does not apply to pyridine and pyrimidine, because QM/FQ and QM/FQF$\mu$ predict almost the same difference in dipole moments. The computed larger solvatochromism for QM/FQF$\mu$ is thus probably due to the fact that both GS ans ES dipole moments are increased of almost 50% with respect to QM/FQ. Howevef, it is worth noticing that molecular dipole moments are highly affected by non-electrostatic interactions, in particular repulsion energy terms, as recently reported by some of the present authors.Giovannini, Lafiosca, and Cappelli (2017) The inclusion of such contributions can in principle reduce the molecular dipole moments of both GS and ES. As a consequence, non-electrostatic interactions should reasonably shift QM/FQF$\mu$ values towards the experimental data. On the other hand, in light of such considerations, the good performance of QM/FQ for pyridine and pyrimidine might be ascribed to error cancellations between the description of electrostatic effects and the lack of explicit consideration of exchange-repulsion effects.Giovannini, Lafiosca, and Cappelli (2017); Curutchet et al. (2018); Giovannini et al. (2019a)

V Summary and Conclusions

In this work, we have extended the fully polarizable QM/FQF$\mu$ approach to the evaluation of excited state energies. The peculiarity of QM/FQF$\mu$ stands in the fact that each MM atom is endowed with both a charge and a dipole, which can vary as a response of the external electric potential and field. QM/FQF$\mu$ is therefore a refinement of QM/FQ, in which only a fluctuating charge is placed on each MM atom.

The extension of QM/FQF$\mu$ to vertical excitation energies, is achieved by first resorting to linear response (LR) theory for SCF methods. Second, a correction, namely corrected Linear Response (cLR), is introduced for the first time for both QM/FQ and QM/FQF$\mu$. The necessity of such a correction is due to the intrinsic limitations of the LR regime when applied to polarizable embedding approaches. In fact, LR adequately describes environmental effects in case of transitions involving bright states associated with large transition dipole moments, but it is not able to catch the relaxation of the environment as a response to charge equilibration of the QM density on the specific excited state. cLR can instead describe such a relaxation, that is particularly relevant in case of transitions with low (or null) transition dipole moments, and with a large difference between the ground and excited dipole moments.

To test our approach, we selected three molecules, namely para-nitroaniline, pyridine and pyrimidine, in aqueous solution, of which we studied $\pi\rightarrow\pi^{*}$ transition for pNA, and dark $n\rightarrow\pi^{*}$ transitions in case of both pyridine and pyrimidine. Such transitions were selected because vacuo-to-water solvatochromic experimental data are available in the literature.Marenich, Cramer, and Truhlar (2015) For all the studied systems, both QM/FQ and QM/FQF$\mu$ over perform the implicit QM/PCM approach, due to the fact that QM/PCM cannot properly describe specific solute-solvent interactions. QM/FQF$\mu$ correctly reproduces the experimental solvachromic shift in case of pNA, whereas it overestimates the experimental value for pyridine and pyrimidine. The analysis of ground and excited state molecular dipole moments shows that such a discrepancy is related to the overestimation of molecular dipole moments, probably related to the electrostatic-only description of the solvation phenomenon which is here adopted.

It is worth noticing that the comparison between QM/MM results and experiments can be biased by experimental conditions and error bar. As a matter of fact, QM/classical results might be compared with full-QM calculation on clusters (i.e. structures in which few solvent molecules are treated at the QM level). However, full-QM super-molecule approaches would include non-electrostatic interactions, such as Pauli repulsion, solute-solvent and/or solvent-solvent Charge Transfer, and dispersion (in case correlation is taken into account). Therefore, a good agreement between QM/classical approaches and such reference values could be affected by error cancellation. The inclusion of non-electrostatic interactions, and also vibronic effects,Santoro and Jacquemin (2016) which can in principle increase the agreement between computed and experimental data, will be the topic of future communications.

VI Supplementary Material

QM/FQ and QM/FQF$\mu$ standard deviations on vertical excitations and oscillator strengths for pNA, pyridine and pyrimidine in aqueous solution. QM/FQ and QM/FQF$\mu$ $\omega_{0}$ and cLR stick and convoluted spectra of pNA, pyridine and pyrimidine in aqueous solution.

VII Acknowledgment

We are thankful for the computer resources provided by the high performance computer facilities of the SMART Laboratory (http://smart.sns.it/). CC gratefully acknowledges the support of H2020-MSCA-ITN-2017 European Training Network “Computational Spectroscopy In Natural sciences and Engineering” (COSINE), grant number 765739. TG acknowledges funding from the Research Council of Norway through its grant TheoLight (grant no. 275506)

VIII References