On the Calculation of IR Spectra with a Fully Polarizable QM/MM Approach Based on Fluctuating Charges and Fluctuating Dipoles

TL;DR

This paper extends a fully polarizable QM/MM approach to compute nuclear gradients and IR spectra, demonstrating its effectiveness on various molecules in aqueous solution and comparing it with other models.

Contribution

It introduces analytical equations for energy derivatives in the QM/FQFμ approach and applies them to calculate IR spectra in condensed phases.

Findings

Accurate IR spectra of molecules in solution were obtained.

The approach outperforms or matches existing models like QM/PCM and QM/FQ.

Analytical derivatives enable efficient spectrum calculations.

Abstract

The 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 evaluation of nuclear gradients and the calculation of IR spectra of molecular systems in condensed phase. To this end, analytical equations defining first and second energy derivatives with respect to nuclear coordinates are derived and discussed. The potentialities of the approach are shown by applying the model to the calculation of IR spectra of Methlyoxirane, Glycidol and Gallic Acid in aqueous solution. The results are compared with the continuum QM/PCM and the polarizable QM/FQ, which is based on Fluctuating Charges only.

Equations

$${\rho }_{q_i}(\text{r})$$

$${\rho }_{{\mathbf{\mu }}_i}(\text{r})$$

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

$$+\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$$

$${\text{T}}_{ij}^{qq}$$

$${\text{T}}_{ij}^{q\mu }$$

$${\mathbf{T}}_{ij}^{\mu \mu }$$

$$+\frac{4}{\sqrt{\pi }{R}_{{\mu }_i-{\mu }_j}^3}\frac{{\text{r}}_{ij}\otimes {\text{r}}_{ij}}{|{\text{r}}_{ij}{|}^2}\exp {\left(-\frac{|{\text{r}}_{ij}|}{R_{{\mu }_i-{\mu }_j}}\right)}^2$$

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

$$=\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 },{\lambda }_{\alpha }\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)\qquad\Rightarrow\qquad\mathbf{D}\mathbf{L}_{\lambda}=-\mathbf{C}_{Q}$$

$$E_{QM/MM}=\sum _{i}V[{\rho }_{QM}]({\text{r}}_i)q_i-{\mathbf{\mu }}_i^{†}\mathbf{E}[{\rho }_{QM}]({\text{r}}_i)$$

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

$$+{\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 }$$

$$\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}\\ \textbf{Q}_{\text{tot}}\\ \mathbf{0}\end{array}\right)+\left(\begin{array}[]{c}-\textbf{V}(\textbf{P})\\ \textbf{0}\\ \textbf{E}(\textbf{P})\end{array}\right)\qquad\Rightarrow\qquad\mathbf{D}\mathbf{L}_{\lambda}=-\mathbf{C}_{Q}-\mathbf{R}(\mathbf{P})$$

$${\mathcal{E}}^x(\mathbf{P},\mathbf{q},\mathbf{\mu },\mathbf{\lambda })=\frac{\partial \mathcal{E}}{\partial x}+\frac{\partial \mathcal{E}}{\partial \mathbf{P}}\frac{\partial \mathbf{P}}{\partial x}+\frac{\partial \mathcal{E}}{\partial \mathbf{q}}\frac{\partial \mathbf{q}}{\partial x}+\frac{\partial \mathcal{E}}{\partial \mathbf{\mu }}\frac{\partial \mathbf{\mu }}{\partial x}+\frac{\partial \mathcal{E}}{\partial \mathbf{\lambda }}\frac{\partial \mathbf{\lambda }}{\partial x}$$

$$\frac{\partial \mathcal{E}}{\partial x}=\text{mbox}tr{\mathbf{h}}^x\mathbf{P}+\frac{1}{2}\text{mbox}tr{\mathbf{G}}^{(x)}(\mathbf{P})\mathbf{P}+{\mathbf{q}}^{†}{\mathbf{V}}^{(x)}(\mathbf{P})-{\mathbf{\mu }}^{†}{\mathbf{E}}^{(x)}(\mathbf{P})$$

$$V_i^{(x)}(\mathbf{P})$$

$$=$$

$$E_i^{(x)}(\mathbf{P})$$

$$-\mathbf{P}\tilde{\mathbf{F}}\mathbf{P}{\mathbf{S}}_{oo}^x=-\tilde{\mathbf{W}}{\mathbf{S}}_{oo}^x$$

$${\mathcal{E}}^x(\mathbf{P},\mathbf{q},\mathbf{\lambda })=\text{mbox}tr{\mathbf{h}}^x\mathbf{P}+\frac{1}{2}\text{mbox}tr{\mathbf{G}}^{(x)}(\mathbf{P})\mathbf{P}+{\mathbf{q}}^{†}{\mathbf{V}}^{(x)}(\mathbf{P})-{\mathbf{\mu }}^{†}{\mathbf{E}}^{(x)}(\mathbf{P})-\text{mbox}tr\mathbf{W}{\mathbf{S}}_{oo}^x$$

$${\mathcal{E}}^{xy}$$

$$\mathbf{D}{\mathbf{L}}^x=-{\mathbf{R}}^{(x)}(\mathbf{P})-\mathbf{R}({\mathbf{P}}^x)$$

$${\tilde{F}}_{\mu \nu }^x={\tilde{F}}_{\mu \nu }^{(x)}+G_{\mu \nu }({\mathbf{P}}^x)-{\mathbf{R}}_{\mu \nu }^{†}{\mathbf{D}}^{-1}\mathbf{R}({\mathbf{P}}^x)$$

$${\tilde{F}}_{\mu \nu }^{(x)}=h_{\mu \nu }^x+G_{\mu \nu }^{(x)}(\mathbf{P})+{\mathbf{L}}^{†}{\mathbf{R}}_{\mu \nu }^x+{\mathbf{R}}_{\mu \nu }^{†}{\mathbf{L}}^{(x)}$$

$${ϵ}_i{P}_{ia}^x-{ϵ}_a{P}_{ia}^x=-{\tilde{Q}}_{ia}$$

$${\tilde{A}}_{ia,jb}=({ϵ}_a-{ϵ}_i){\delta }_{ij}{\delta }_{ab}+⟨aj||ib⟩-{\mathbf{R}}_{ia}^{†}{\mathbf{D}}^{-1}{\mathbf{R}}_{jb}$$

$${\tilde{B}}_{ia,jb}=⟨ab||ij⟩-{\mathbf{R}}_{ia}^{†}{\mathbf{D}}^{-1}{\mathbf{R}}_{bj}$$

$$\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)=\left(\begin{array}[]{c}\mathbf{Q}\\ \mathbf{Q}^{*}\\ \end{array}\right)$$

$${\tilde{A}}_{ia,jb}=({ϵ}_a-{ϵ}_i){\delta }_{ij}{\delta }_{ab}+⟨aj||ib⟩-{\mathbf{R}}_{ia}^{†}{\mathbf{D}}^{-1}{\mathbf{R}}_{jb}$$

Full text

1 Introduction

Vibrational spectroscopy, in particular infrared spectroscopy, is one of the most common techniques to study structural and dynamical features of molecular systems. Experimental spectra can be affected by a combination of effects, ranging from anharmonicity to solvent effects, the latter playing a relevant role because most experiments are conducted in the condensed phase.1, 2, 3, 4, 5, 6 In this work, we especially focus on the development of a method to account for the mutual interaction between a molecular systems and its environment and its effect on the prediction of IR spectra. In fact, the presence of the environment can alter the electronic response of the target molecule to the external electric field and the vibrational frequecies associated with the normal modes. Therefore, approaches able to accurately describe environmental effects are required to obtain computed spectra directly comparable with experiments.

In the computational practice, the effects of the environment on a given spectral property are usually included by resorting to focused models,7, 8, 9, 10, 11, 12 which are based on the assumption that the spectral signal is essentially due to the target molecule (i.e. the solute in case of solutions) and the environment (i.e. the solvent) only modifies but not determines it.

Besides the commonly used continuum solvation approaches 13, 14 the family of QM/MM methods,15, 16, 7, 8 may be exploited. Their quality is connected to the specific force field which is exploited to treat the MM portion, and on the approach which is used to define the QM/MM interaction. The latter can be modeled by means of the basic mechanical embedding approach or by resorting to the so-called electrostatic embedding, where the coupling between QM and MM portions is described in terms of the Coulomb law, i.e. the electrostatic interaction between the potential of the QM density and the fixed charges which are placed at MM atoms. This picture is refined in the so-called polarizable QM/MM approaches, in which the mutual QM/MM polarization is taken into account; it can be modeled in different ways, e.g. by resorting to distributed multipoles,17, 18, 19, 20, 21 induced dipoles,22, 23, 24, 25 Drude oscillators26, Fluctuating Charges (FQ)27, 28, 29 and the recently developed approach based on both Fluctuating Charges and Fluctuating Dipoles (FQF$\mu$).30

In QM/FQF$\mu$ the MM portion is described in terms of both fluctuating charges and fluctuating dipoles, which are place at MM atoms positions and can vary as a response to the QM electric potential (the FQs) and MM atomic electronegativities and QM electric field (the F$\mu$s), respectively. Such an approach is a pragmatical extension of the QM/FQ approach, previously developed by some of the present authors, 31, 29, 32, 31, 33, 34, 35, 36 in which the MM portion is described by means of electric charges which can be polarized by the QM density and viceversa. As a consequence, QM/FQF$\mu$ takes into account both the out-of-plane and anisotropy contributions to polarization thanks to the inclusion of the electric dipoles in the MM portion. It is worth remarking that similar approaches have been proposed in other contexts,37, 38, 39, 40, 41, 42 however they are not based on a variational formalism and therefore are not specifically intended to model molecular properties/spectra. Also, to the best of our knowledge, they have never been extended to the calculation of molecular properties/spectra determined by the nuclear response to external fields.

The quality of QM/FQF$\mu$ at predicting electrostatic interaction energies has been recently discussed 30 and some of the present authors have recently extended this approach to the calculation of electronic vertical transition energies of organic molecules in solution at the TD-DFT level.43 In this work, QM/FQF$\mu$ is further extended to the calculation of IR spectra, through its extension to energy nuclear derivatives. Remarkably, other QM/MM approaches have been extended to the calculation of energy gradients,44, 45, 46, 47, 48, 49, 50, 51 however the only previous polarizable QM/MM approach extended to vibrational spectroscopy is the QM/FQ method developed by some of us.36, 35, 34, 52

The manuscript is organized as follows. In the next section, the FQF$\mu$ force field is presented and its coupling with a QM wavefunction defined at the SCF level (QM/FQF$\mu$) is detailed. Equations for analytical first and second energy derivatives are then presented and discussed. After a brief section discussing on the computational protocol which is adopted, numerical results are presented. In particular, QM/FQF$\mu$ is challanged against the description of IR spectra of three organic molecules in aqueous solution, namely methyloxirane, glycidol and gallic acid. IR spectra are computed by exploiting a hyerarchy of polarizable embedding approaches, namely QM/PCM, QM/FQ and QM/FQF$\mu$. Computed spectra are compared to their experimental counterparts, which are taken from the literature.53, 54, 55 From the comparison between the different polarizable approaches, differences arising from the moving from continuum PCM to QM/FQ and QM/FQF$\mu$ fully atomistic approaches are highlighted, and the role of fluctuating charges, and their coupling with fluctuating dipoles is also appreciated. Some drawn conclusions and the discussion on future perspectives of this approach end the manuscript.

2 Theoretical Model

2.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 can vary according to the external electric potential and electric field. Both charges and dipoles are described as s-type gaussian distribution functions:

[TABLE]

where $R_{q_{i}}$ and $R_{\bm{\mu}_{i}}$ are the width of the Gaussian distributions $\rho_{q_{i}}$ and $\rho_{\bm{\mu}_{i}}$, respectively. $\hat{\textbf{n}}_{i}$ is a unit vector pointing to the dipole direction $\bm{\mu}_{i}$.

The total energy $E$ associated with a distribution of charges and dipoles can be written as: 40, 30

[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. If the gaussian distributions in Eq. 1 are adopted, the functional form of the interaction kernels provided by Mayer 40 can be exploited (see also Ref. 30):

[TABLE]

In order to collect all quadratic terms in the charges, the diagonal elements of $\textbf{T}^{qq}$ can be imposed to be equal to the atomic chemical hardness $\eta$, so that the width of the charge distribution $R_{q}$ is defined without the need of any parametrization.30 The same holds for the diagonal elements of $\textbf{T}^{\mu\mu}$ and the dipole distribution $R_{\mu}$, which can be defined in terms of the atomic polarizabilities ($\alpha$).30 The definition of the gaussian width $R_{q_{i}}$ and $R_{\mu_{i}}$ in terms of $\eta_{i}$ and $\alpha_{i}$ limits the number of parameters which enter the definition of FQF$\mu$ to electronegativity, chemical hardness and polarizability for each atom type. Therefore, Eq. 2 can be rewritten as:

[TABLE]

where a matrix notation has been adopted. In Eq. 6, charges are not forced to any value by external constraints and the equilibrium condition is reached when the Electronegativity Equalization Principle (EEP) is satisfied. If each MM molecule is constrained to assume a fixed, total charge equal to Qα, Eq. 6 can be written by exploiting a set of Lagrangian multipliers ($\lambda_{\alpha}$), whose number is equal to the total number of molecules in the MM portion.

[TABLE]

where $\alpha$ and $\beta$ run over molecules and the constraints $\lambda_{\alpha}$ are meant to preserve the total charge $Q_{\alpha}$ of each molecule. Therefore, the constrained minimum is found by imposing all derivatives of $\mathcal{F}$ with respect all variables to be equal to zero, thus resulting in the following linear system30:

[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$ can be effectively coupled to a QM SCF wavefunction in a QM/MM framework. The QM density interacts as a classical density of charge with both charges and dipoles:

[TABLE]

where $V[\rho_{QM}](\textbf{r}_{i})$ and $\mathbf{E}[\rho_{QM}](\textbf{r}_{i})$ are the electric potential and electric field, respectively, calculated at the $i$-th charge and $i$-th dipoles placed at $\textbf{r}_{i}$. Finally, the global QM/MM energy functional for a SCF-like description of the QM portion finally reads:

[TABLE]

where $\mathbf{h}$ and $\mathbf{G}$ are the usual one- and two-electron matrices. The effective Fock matrix is defined as the derivative of the energy with respect to the density matrix:

[TABLE]

where the interaction of the electron density with both charges and dipoles are included through the coupling electrostatic terms. 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. 8, a new source term $\mathbf{R}(\mathbf{P}$), due to the coupling of both charges and dipole with the SCF density, arises. Again, $\mathbf{L}_{\lambda}$ is a vector containing charges, dipoles and Lagrangian multipliers.

2.2 Analytical Energy Derivatives

In this section QM/FQF$\mu$ analytical first and second energy derivatives are presented and discussed. The following equations are defined in the so-called Partial Hessian Vibrational Approach (PHVA), 56, 57, 58 which has been amply exploited to treat vibrational phenomena of complex systems.52, 35, 34, 36 Within such a framework, it is assumed that the geometrical perturbation only acts on the QM portion of the system, whereas MM atoms are unaffected. Remarkably, the PHVA is fully consistent with a focused approach. For the sake of completeness, however, equations for first and second derivatives of FQF$\mu$ MM atoms are given in the Appendix section. The following derivation directly follows what already reported by some of the present authors in case of QM/FQ.52 This allows to directly identify the additional terms which depend on the presence of fluctuating dipoles in the MM portion.

2.2.1 Energy first derivatives

The energy first derivative of Eq. 10 with respect to the $x$ nuclear displacement can be expressed by means of the chain rule:59, 52

[TABLE]

The last three terms vanish because of the stationarity conditions. The first term, which is the partial derivative of the energy with respect to the position of a QM nucleus, reads:

[TABLE]

where

[TABLE]

The term involving the derivatives of the density matrix can be computed starting from the idempotency condition:

[TABLE]

where the subscript $oo$ denotes the occupied–occupied block of the matrix in the MO basis, and $\mathbf{W}$ is the energy-weighted density matrix contribution. By collecting all the terms:

[TABLE]

Notice that the term ($\mathbf{q}^{\dagger}\mathbf{V}^{(x)}(\mathbf{P})$) is the same as computed in the QM/FQ approach.52 Therefore, the inclusion of fluctuating dipoles gives rise to the additional term $\bm{\mu}^{\dagger}\mathbf{E}^{(x)}(\mathbf{P})$.

2.2.2 Energy second derivatives

The energy second derivative with respect to nuclear displacements $x$ and $y$ is obtained by differentiating eq. 17 and by exploiting once again the chain rule:

[TABLE]

Thus, the derivatives of the off-diagonal blocks of the density matrix and charges/dipoles need to be calculated. Density matrix derivatives can be obtained through a Coupled Perturbed Hartree–Fock (CPHF) or Kohn-Sham (CPKS) procedure.

FQF$\mu$ charge and dipole derivatives can be calculated by differentiating eq. 12:

[TABLE]

The Fock matrix derivative is defined as:

[TABLE]

where $\widetilde{F}_{\mu\nu}^{(x)}$, which collects all explicit derivatives of the Fock matrix, reads:

[TABLE]

By rearranging the terms, the CPHF/CPKS equations are obtained (MO basis):

[TABLE]

By taking the adjunct equation and introducing the following matrices (we assume orbitals to be real):

[TABLE]

[TABLE]

the following equation is obtained:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Therefore, the derivatives of the density matrix and FQF$\mu$ charge/dipole derivatives with respect to QM nuclear positions are obtained by solving Eqs. 24 and 19, respectively. Notice that such a derivation is coherent with what has been reported for linear response in the zero-frequency limit.43

To summarize, FQF$\mu$ contributions to analytical second derivatives can be grouped into three categories:

explicit terms:

[TABLE] 2. 2.

contributions to Fock matrix derivatives:

[TABLE] 3. 3.

additional terms to the CPHF/CPKS matrix:

[TABLE]

Notice that only the last term is required in case of electric perturbations.Also, similarly to what has already been discussed for energy first derivatives, QM/FQF$\mu$ second derivatives differ from QM/FQ ones because additional terms depending on fluctuating dipoles need to be included.52

3 Computational Details

R-methyloxirane (MOXY), (S)-Glycidol (GLY) and Gallic Acid (GA) geometries were optimized at the B3LYP/ aug-cc-pVDZ. Solvent effects on solutes’ geometries were included through the Polarizable Continuum Model (PCM). MOXY and GLY Molecular Dynamics (MD) simulations were carried out as detailed in previous works by some of the present authors60, 36 by using GROMACS.61

A 25 ns MD simulation of GA in aqueous solution was performed by using a similar computational protocol. A GA molecule was placed at the center of a cubic box and solvated with 6025 TIP3P water molecules.62 The parameters used to describe the GA inter-/intra-molecular interactions were taken from the GAFF force field63 by using the ANTECHAMBER package.64 GA molecule was keeped fixed during all the steps of the simulation. To model intermolecular solute-solvent interactions, ga partial charges were computed by relying on the Hirshfeld population analysis,65 as previously done by some of the present authors.66, 31 Partial charges were computed at the B3LYP/6-311++G** level of theory by including solvent effects by means of PCM. Two subsequent 100 ps runs were performed for equilibration purposes by using NVT and NPT ensembles, respectively. A 25 ns NVT MD simulation was then performed to sample the configurartional space at time steps of 1 fs and by saving coordinates every 10 ps. The system was simulated by using three-dimensional periodic boundary conditions; non-bonded interactions cutoff was set to 10 Å. A particle mesh Ewald (PME) correction for the long-range electrostatics was applied and the temperature was maintained at 300 K by using the velocity rescale algorithm.

A total of 200 uncorrelated snapshots were extracted from the last 20 ns, 50 ns and 25 ns of MD simulations in case of MOXY, GLY and GA, respectively. For each snapshot, a sphere centered at the solute’s geometric center was cut.A cutting radius of 12 Å was used for MOXY and GLY, whereas a cutting radius of 15 Åwas used for GA.

QM/FQ and QM/FQF$\mu$ partial geometry optimization for the solute on each snapsho was performed according to the default settings of Gaussian16,67 by keeping all water molecules fixed. Finally, infrared (IR) spectra were calculated on each partially optimized snapshot with the QM/FQ and QM/FQF$\mu$ models at the B3LYP/aug-cc-pVDZ level (MOXY and GLY) and B3LYP/6-311++G** level (GA). SPC parameters for water 27 were used for FQ calculations. The set of parameters recently developed by som of us for QM/FQF$\mu$ calculations in aqueous solutions were exploited.30. IR data were averaged to obtain final spectra for the three solutes; peaks were convoluted with a Lorentzian lineshape, with Full Width at Half Maximum (FWHM) of 4 $\mathrm{cm}^{-1}$. For the sake of comparison, QM/PCM IR spectra were also computed.

All QM/FQ and QM/FQF$\mu$ calculations were performed by using a locally modified version of the Gaussian 16 package,67 where QM/FQF$\mu$ analytical energy first derivatives were implemented.

4 Numerical Results

In this section, the results obtained by applying QM/FQF$\mu$ to the calculaiton of IR spectra of MOXY, GLY and GA in aqueous solutions are reported (see Fig. 1, panels a-c for molecular structures). MOXY is a widely exploited test system for computational models68, 69, 70, 71, 72, 73, 74. GLY, which bears an additional hydroxil group, has been previously studied both theoretically and experimentally.36, 54 In particular, it has been shown that eight different GLY conformers exist in aqueous solution, thus its theoretical modeling is challenging.75, 54, 36 GA is an organic acid characterized by the presence of three hydroxil groups linked to the aromatic ring. The modeling of IR spectra of GA in aqueous solution is also challenging, because it has been previously reported by one of the present authors that the implicit PCM fails at reproducing the experimental IR spectrum, and that experimental spectral features can be recovered by adopting a supramolecular approach which includes eight explicit water molecules in the QM portion.53 Therefore, it appears to be an ideal test-bed for QM/FQF$\mu$.

QM/FQF$\mu$ spectra will be compared to QM/FQ, which only considers fluctuating charges on MM atoms. The rerason for such a comparison is twofold: (i) QM/FQF$\mu$ is formally an extension of QM/FQ, and (ii) QM/FQ has been successfully applied to vibrational spectra of molecular systems in aqueous solution.35, 34, 36 Thus, the comparison between the two approaches can directly quantify the relevance of fluctuating dipoles in the description of vibrational spectra and report on the performance of the novel QM/FQF$\mu$ method.

4.1 Methyloxirane in Aqueous Solution

Figure 2 compares QM/FQ (top) and QM/FQF$\mu$ (bottom) stick and convoluted spectra for 200 snapshots extracted from the MD simulation. Such a number of snapshots is enough to reach convergence.34, 36 Stick spectra are obtained by plotting raw data extracted from each frequency calculation. Figure 2 clearly shows that both QM/FQ and QM/FQF$\mu$ IR wavenumbers and dipole strengths depend on the snapshot, i.e. on the arrangement of water molecules around the solute. As compared to QM/FQ, QM/FQF$\mu$ exhibits a larger spread, in particular in the vibrational wavenumbers. The largest variability of QM/FQF$\mu$ sticks occurs in the region between 1450 and 1500 cm*-1*, which is associated to methyl and CH bending modes (see Figures S1 in the Supporting Information (SI) for a pictorial view of the normal modes for a randomly chosen snapshot).

Clearly, band broadening is automatically considered in both QM/FQ and QM/FQF$\mu$ approaches coupled with the dynamical description given by the MD simulation, which samples over the solute-solvent phase-space. Therefore, solvent inhomogeneous broadening (due to the fluctuations of the solvent molecules) does not need to be artificially considered by imposing a pre-defined (and arbitrary) band-width, which is instead necessary when other static approaches, such as PCM, are used. The lorentzian convolution obtained by using a FWHM of 4 cm*-1* is plotted for both QM/FQ and QM/FQF$\mu$ spectra. Notice that QM/FQ IR spectrum is identical to what we reported in case of the three layer QM/FQ/PCM approach, where the PCM is used as third layer and coupled to both QM and FQ portions.Such a similarity means that water molecules explicitly included in snapshots are able to also account for bulk water effects.

We now move to the comparison between computed and experimental spectra (see Fig.3). Notice that the experimental spectrum was measured for neat liquid MOXY rather than aqueous solution;55 therefore, some discrepancies with our computed results need to be expected.

The computed and experimental IR spectra are dominated by an intense band at about 850 cm*-1*. This signal is given by the symmetric stretching of the C-O bond of the epoxyl group. QM/FQF$\mu$ band is blueshifted with respect to QM/FQ, thus the inclusion of fluctuating dipoles increases solute-solvent interactions. In fact, assuming that the normal modes are equally predicted by both QM/FQ and QM/FQF$\mu$ approaches, an increase in the solute-solvent interactions is reflected in a decrease in the solute intramolecular bond strengths, therefore, resulting in a blueshifted. Similar considerations apply also to the regions between 900-1100 cm*-1*, where the normal modes involve vibrations of the MOXY oxygen atom. In the other regions of the spectra, such a blueshift is instead not recorded. This can be explained by the fact that the normal modes do not involve vibrations of the oxygen atom, which is the only MOXY atom potentially exhibiting an Hydrogen Bond with the solvent molecules.

Further differences between QM/FQ and QM/FQF$\mu$ are computed for the inhomogeneous band broadening, which is almost absent in case of QM/FQ, whereas it affects almost all QM/FQF$\mu$ bands. This is not unexpected if the raw data depicted in Fig. 2 are considered. In fact, as already pointed out, QM/FQF$\mu$ generally spreads a larger energy region for each vibrational normal mode.

The major discrepancies with the experiment are reported in case of QM/FQ and QM/FQF$\mu$ in case of the inhomogenous band broadening, which is almost absent in the experiment. Thus, the main feature added by aqueous solution seem to be a larger broadening of vibrational bands.

4.2 Glycidol in Aqueous Solution

Similarly to MOXY, QM/FQ and QM/FQF$\mu$ IR spectra of GLY in aqueous solution were calculated by averaging over 200 snapshots extracted from the MD simulation.36 QM/FQ and QM/FQF$\mu$ raw data are graphically plotted in Fig. 4, together with their lorentzian convolution. As stated before, the case of GLY in aqueous solution is far more complicated than MOXY, because GLY is a flexible molecule, i.e. it exists in different conformations. This is reflected in the stick spectra depicted in Fig. 4, which show a larger variability both in intensities and wavenumbers as compared to MOXY (see Fig. 2); this applies to both QM/FQ and QM/FQF$\mu$ calculations with the exception of the region around 1110 cm*-1*. There, QM/FQ shows a substantial variation in intensity, whereas QM/FQF$\mu$ spreads a larger wavenumber range. The spreading of the intensities is instead larger for QM/FQF$\mu$ in the region between 1400-1600 cm*-1*. Despite such discrepancies between QM/FQ and QM/FQF$\mu$, inhomogeneous broadening is described by both approaches.

QM/PCM, QM/FQ and QM/FQF$\mu$ convoluted IR spectra are compared to experiments in Fig. 5. The experimental IR spectrum i reproduced from Ref. 54. Normal modes for a randomly chosen snapshot extracted from the MD are depicted in Figure S2 in the S for the region 700-1800 cm*-1* .

Computed and experimental IR spectra are dominated by an intense band at about 1050 cm*-1*, assigned to a diffuse stretching/bending normal mode, involving the hydroxyl group. QM/FQF$\mu$ spectrum is generally blueshifted with respect to both QM/FQ and QM/PCM, probably due to the fact it predicts larger solute-solvent interactions.

Moving to the comparison with experimental data, both the experimental and computed spectra are characterized by a two-peak-shaped band between 1200 and 1300 cm*-1*, which is assigned to the C-OH and C-CH bending modes (at about 1230 and 1270 cm*-1*, respectively). Furthermore, above 1400 cm*-1* a three-peak-shaped band can be identified, due to the C-OH bending (1395 cm*-1*), a diffuse C-CH bending (1440 cm*-1*) and a CH2 bending (1465 cm*-1*).

QM/FQ, QM/FQF$\mu$ and the experimental IR spectru are nicely in agreement. In fact, most of relative intensities and the band broadening are correctly reproduced. In particular, QM/FQ accurately predicts the two-peak band between 1200 and 1300 cm*-1*, whereas QM/FQF$\mu$ is able to catch the inhomogeneity of the three-peak-shaped band between 1400 and 1500 cm*-1*. Some discrepancies are reported in case of the peak at about 1100 cm*-1* (due to the CH scissoring), which is predicted to have a very low in intensity by both atomistic QM/MM approaches, whereas it is the second most intense peak in the experimental spectrum. Notice that it has been reported that the relative intensity of this peak can be correctly reproduced if a supermolecule approach is adopted,, i.e. if water molecules are included in the QM region.54 These findings, together with the results obtained by adopting our QM/classical modeling, suggest that the inclusion of non-electrostatic interactions, which are considered in the full QM supermolecule approach, can play a relevant role to improve the quality of the computed spectrum in this region. Overall, it is worth noticing that the continuum QM/PCM approach cannot reproduce the experimental spectrum (see top of Fig. 5), thus remarking once again the huge potentialities of our atomistic QM/FQ and QM/FQF$\mu$ approaches to model vibrational spectra of solutes strongly interacting with the aqueous environment.

To end the discussion on the IR spectrum of GLY in aqueous solution, we point out that the broad band measured between 1600-1700 cm*-1* in the experiment is not reproduced by any of the selected QM/classical approaches. As already reported by some of the present authors 35, 36 and in Refs. 76, 77, 54, such band is due to the OH bending mode of water molecules linked to GLY; therefore, it cannot be modeled by our approaches, in which the normal modes/frequencies of the environment are not computed.

4.3 Gallic Acid in Aqueous Solution

4.3.1 MD Analysis

before discussing GA IR spectra (see Fig. 6 for atom labelling), in this section the MD trajectory is examined in terms of both radial distribution functions (rdfs) and running coordination numbers (RCNs). In particular, in order to obtain a description of the solvent local structure and to analyze hydrogen bonding patterns between GA and water molecules, intermolecular H(GA)$\cdots$OW and O(GA)$\cdots$HW rdfs and the corresponding running coordination numbers were calculated and are reported in Figs. 8 and 7 respectively. The coordination number of a specific site is defined by combining the distance at which the first/second minimum of rdf occurs with the corresponding distance in the running coordination number.

The oxygen atoms of the three hydroxyl groups, i.e. O1, O3 and O4 (see Fig. 6), present a radial distribution function with a similar shape, consisting of a broad peak at about 3.4 Å, which is associated to the second solvation shell. The coordination number of these sites is equal to 16.7, 12.8 and 12.2, respectively. Remarkably, O3 rdf shows a peak at about 2.0 Å, due to the fact that this is the only one among the hydroxyl groups that can form intermolecular HB with water molecules. O1 and O4 are involved in an intramolecular HB. Moving to the acid group, O5 rdf presents a sharp peak related to the first solvation shell, with a maximum at 2.0 Å  and a coordination number corresponding to 1.6; the plot for O2 is very similar to what has already been discussed in case of O1 and O4.

The specific arrangement of the hydroxyl groups in the meta and para positions (with respect to the acid) of GA induces a weak but detectable anisotropy on the rdfs of both hydrogen and oxygen atoms. In fact, both H4 and H5 rdfs present a similar shape, with a first peak at about 2.1 and 2.2 Å, respectively. They correspond to coordination numbers of 0.8 and 0.7, respectively. A second peak related to the second solvation shell is placed at 3.68 Å , which corresponds to coordination numbers of 7.5 and 6.3, respectively. The H3 and H6, instead, show a pronounced first narrow peak at about the same distance (1.73 Å), with a coordination number of 1.0 and a second peak around 3.97 Å  with coordination numbers of 15.2 and 13.8, respectively. Again, the differences between H3 and H4/H5 are due to the fact that H4 and H5 are involved in the intramolecular HB, whereas H6 is free to form intermolecular HB with water molecules. Remarkably, the results here discussed are similar to the findings previously reported byone of the present authors.53

4.3.2 IR spectrum of GA in Aqueous Solution

QM/FQ and QM/FQF$\mu$ IR spectra of GA were obtained by averaging over 200 snapshots extracted from the MD run. Similarly to MOXY and GLY, we checked that such a number of snapshots produce converged spectra. The raw data extracted from the calculations, i.e. stick spectra are reported in Fig. 9 for the region 1000-1800 cm*-1* (i.e. the region of interest for the experimental investigation, vide infra); lorentzian convolution (FWHM = 4 cm*-1*) is also depicted. Both QM/FQ and QM/FQF$\mu$ stick spectra show a large spreading in intensities and frequencies. Such a feature is reported for most computed bands, in particular in the region between 1100-1400 cm*-1*, in which single bands are not easily detectable (see for comparison MOXY and GLY raw spectra in Figs. 2 and 4).

The comparison between QM/PCM, QM/FQ, QM/FQF$\mu$ and the experimental IR spectrum 78 shows that the latter is dominated by a three-band broad structure between 1200 and 1500 cm*-1*, which probably involve more than one normal mode. A well-separated peak is present at 1000 cm*-1* and it is associated to the C-OH bending (see Fig. S3 in the SI for a graphical depiction of the normal modes). Moreover, three small bands of the same intensity are reported in the region 1500-1800 cm*-1*, which are mainly due to composite C-OH bending modes of the hydroxyl groups and the acidic C=O stretching.

The QM/PCM spectrum is dominated by two peaks, placed at about 1170 cm*-1* and 1750 cm*-1*, respectively. Such peaks are related to a composite C-OH bending and to the C=O stretching. A similar spectrum is predicted by adopting the atomistic QM/FQ approach, in which the most intense peak is predicted in case of the C=O stretching at about 1750 cm*-1*, whereas the peaks at about 1200 and 1400 cm*-1* have almost the same intensity. It is worth noticing that in QM/FQ spectra most bands present an inhomogenous broadening that is related once again to the dynamical picture given by the sampling of the phase-space through MD. In addition, similarly to MOXY and GLY, most of the computed QM/FQ bands are blueshifted with respect to their QM/PCM counterparts, thus reflecting the stronger solute-solvent interaction.

Most QM/FQF$\mu$ bands are blueshifted with respect to QM/PCM, whereas they are redshifted with respect to QM/FQ values, thus highlighting the different electrostatic description given by the two explicit approaches. Remarkably, similarly to GLY in aqueous solution, vibrational frequencies are not completely in agreement with the experimental ones, that probably due to the lack of anharmonicity effects and the use of DFT. Moreover, inhomogenous band broadening is correclty repoduced by both atomistic approaches, thus resulting in a very good agreement with the experiments, in particular for the experimentally most intense band at about 1350 cm*-1*.

By further deepening the analysis of computed spectra, we see that QM/FQF$\mu$ IR spectrum is dominated by three bands at about 1200, 1350 and 1700 cm*-1*, which are predicted almost with the same intensity. This is a big improvement with respect to both QM/PCM and QM/FQ approaches, because the most intense band in the experimental spectrum is correctly predicted only by QM/FQF$\mu$. It is worth noticing that a correct reproduction of the intensity of this peak was achieved by some of the present authors by resorting to a supermolecule approach (here called QM/QM${}_{\text{w}}$/PCM), i.e. by including 8 QM water molecules in the definition of the QM solute in QM/PCM calculations (see Fig. 11). Such an observation indeed indicates that an explicit solvation approach is needed to recover the experimental features in this region. In addition, due to the fact that QM/FQF$\mu$ appropriately reproduces the most intense band of the experimental spectrum, in a similar way as the supermolecule approach, we can conclude that the electrostatic description of the HB interaction is the most relevant solute-solvent contribution, and that the electrostatic description given by QM/FQF$\mu$ in this case overcomes that modeled by the QM/FQ. The largest discrepancies between QM/FQF$\mu$ and QM/QM${}_{\text{w}}$/PCM are reported for normal mode frequencies, which are better reproduced by the supramolecule approach, in particular in the region between 1200 and 1400 cm*-1*. Such an improvement can be related to the fact that the supermolecule accounts for non-electrostatic interactions (in particular Pauli repulsion) which can therefore play a relevant role in the determination of vibrational frequencies.

To conclude the discussion on GA, it is worth noticing that all the considered computational approaches predict very large intensities for the three modes in the 1500-1800 cm*-1* region, even the supramolecule QM/QM${}_{\text{w}}$/PCM. This is probably due to both the lack of vibrational anharmonicity, which has been reported to affect not only frequencies but also intensities.1, 2, 3

5 Summary and Conclusions

In this work, the fully polarizable QM/FQF$\mu$ approach, recently developed by some of the present authors, has been extended to the evaluation of IR spectra, though the development of analytical energy first and second derivatives. In QM/FQF$\mu$ both a charge and a dipole, which can vary as a response to the external electric field and potential, are placed on each atom of the MM portion. Such a model can be viewed as a refinement of the QM/FQ approach, in which only fluctuating charge are considered.

Our approach has been tested against the reproduction of IR spectra of three systems in aqueous solution, namely methyloxirane, glycidol and gallic acid. The selected molecules can interact with water through strong solute-solvent interactions; this is reflected in the computed IR spectra by the fact that the atomistic QM/FQ and QM/FQF$\mu$ generally over-perform the implicit QM/PCM. In particular, and as expected, the bands which are mostly affected by the atomistic description of the environment, are those involving the polar moieties of the investigated molecular systems.

In case of both methyloxirane and glycidol in aqueous solution, QM/FQF$\mu$ predicts similar spectra with respect to QM/FQ, whereas for gallic Acid the inclusion of anisotropic terms in the MM modeling, i.e. the inclusion of fluctuating dipoles, permits to obtain a better agreement with the experimental data. In the latter case, we also noticed that some bands can probably be affected by anharmonicity79 and non-electrostatic solute-solvent interactions. The extension of QM/FQF$\mu$ so to include anharmonicity and non-electrostatic interactions, for instance by extending the method already developed by some of the present authors,80, 81, 82 might be beneficial and will be the subject of future communications. Moreover, it is worth pointing out that the development and implementation of analytical first energy derivatives, i.e. energy gradients, is not only the basic ingredient for computing vibrational spectra, but it also allows for a further extension of the model to QM/MM MD, as already reported for other kinds of polarizable QM/MM approaches.48, 49, 50, 51

6 Appendix

6.1 FQF$\mu$ Energy First Derivatives with respect to MM coordinates

The derivative of the energy with respect to the position of an MM atom, which we will denote with the superscript $\xi$, can be obtained by using the chain rule. Only explicit contributions arise, as the overlap matrix does not depend on MM atom positions. In fact,

[TABLE]

where the derivative of the QM/MM interaction potential is equal to the electric field produced by the QM density acting on the charges, whereas the derivative of the QM/MM interaction field is the electric field gradient acting on the dipoles.

The derivatives of the interactions kernels $\mathbf{T}^{\xi}_{qq}$, $\mathbf{T}^{\xi}_{q\mu}$ and $\mathbf{T}^{\xi}_{\mu\mu}$ can be obtained by differentiating Eqs. 3, 4 and 5:

[TABLE]

where $\mathbf{K}=\textbf{r}_{i,j}\otimes\textbf{r}_{i,j}$ and $\mathbf{K}^{(x_{i})}$ is the derivative of $\mathbf{K}$ with respect to the $x$ component of the $i$-th element. For the sake of clarity, in Eq. 31 $R_{\mu_{i}-\mu_{j}}$ was substituted by $R$.

6.2 FQF$\mu$ Energy Second Derivatives with respect to MM coordinates

For the sake of completeness, in this appendix the formulas for the full Hessian matrix, i.e. including also the QM/MM and MM/MM blocks, are given. Derivatives with respect to MM atom coordinates will be denoted by the superscripts $\xi,\eta$. The QM-MM block of the Hessian can be obtained by differentiating once the forces on the MM portion with respect to the position of a QM nucleus:

[TABLE]

where the last term vanishes. Substituting $\mathcal{E}^{\xi}$ with Eq 28:

[TABLE]

The derivatives of the density matrix and of the FQs can be obtained by solving the CPHF equations described in Section 2.2.2: therefore, there is no need to enlarge the CPHF system of equations to calculate the derivatives of the density matrix with respect to the positions of the MM atoms. This is, however, unavoidable when the MM-MM block of the Hessian is to be calculated. By differentiating Eq. 28 with respect to the position of MM atoms:

[TABLE]

The last term vanishes once again; to calculate the derivatives of the charges and the density, a new set of CPHF equations needs to be solved. Differentiation of the Liouville equation (as the overlap does not depend on the positions of the MM atoms, we can work in the MO basis) and projection onto the o-v block gives:

[TABLE]

The Fock matrix derivatives have no contributions arising from the one- and two-electron matrices, but only from the density and FQF$\mu$ derivatives:

[TABLE]

By differentiating the FQF$\mu$ equations we obtain:

[TABLE]

By putting everything together, a Casida-like system of equations is obtained, where the matrices are defined in Eq. 22 and the right-hand side reads:

[TABLE]

7 Supporting Information

Graphical depiction of normal modes of MOXY, GLY and GA in aqueous solution.

8 Acknowledgments

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.