Development, Validation, and Pilot Application of a Generalized Fluctuating Charge Model for Computational Spectroscopy in Solution.

A general approach enforcing nonperiodic boundary conditions for the computation of spectroscopic properties in solution has been improved including an effective description of charge-transfer contributions and coordination number adjustment for explicit solvent molecules. Both contributions are obtained from a continuous description of intermolecular hydrogen bonds, which has been employed also for an effective clustering of molecular dynamics trajectories. Fine tuning of the model has been performed for several water clusters, and then its efficiency and reliability have been demonstrated by computing the absorption spectra of different creatinine tautomers in aqueous solution.

Introduction

The computation of spectroscopic parameters in solution can be formally separated in two steps. The first step is the generation of a sufficient number of representative structures by means of stochastic simulations (e.g., Monte Carlo or molecular dynamics), which are then used in the second step to compute averaged values of the sought property. While in principle both steps can be performed at the same computational level,[1] the accuracy required in the second step is usually much higher than that needed in the generation of representative snapshots. As a consequence, quite simple force fields are now available for describing solvent dynamics and also solute–solvent interactions can be modeled quite effectively by nonpolarizable force fields or, at most, by quite simple QM-MM approaches. The situation is different for the computation of spectroscopic parameters where solvent polarization can have a non-negligible effect. Polarizable MM and QM/MM approaches can be based on distributed multipoles,[2] induced dipoles,[3,4] Drude oscillators,[5] or fluctuating charges (FQ).[6] We have followed the FQ route in view of its remarkable efficiency, close connection with the polarizable continuum description of bulk solvent effects, and straightforward improvement by means of virtual sites. In this connection, we have developed in recent years an overall computational protocol enforcing nonperiodic boundary conditions, which usually provides remarkably accurate results.[7−9] Some not fully satisfactory aspects are, however, still present, which concern a more effective treatment of reaction field effects by bulk solvent together with the lack of charge transfer and proper coordination number for solvent molecules near the outer boundary. In the present paper, we will be concerned with the inclusion of both contributions in the stage of property evaluation, whereas a companion paper is devoted to the stochastic simulation step. In particular, reaction field effects are treated by the very powerful ddCOSMO model,[10] whereas charge transfer and coordination number saturation are included by a new recipe described in detail in the next section. Then, after validation of the whole procedure for water clusters, the absorption spectra of amine and imine tautomers of creatinine in aqueous solution are used to show the potentialities of the new procedure. In this connection, effective strategies for clustering the snapshots generated by the stochastic simulation are discussed, which allow to strongly reduce the number of expensive QM/MM computations of spectroscopic parameters.

Theory and Computational Details

Theory

In the following, a generalized FQ model will be presented in order to deal with charge transfer and outer boundary effects in polarizable QM/MM calculations. The basic idea behind such developments is that the geometry of the molecular systems under study allows the estimate of the intermolecular bonding by means of a suitable function (in the present case, the hydrogen bond function proposed by Pagliai et al.[11]), from which the coordination number of the fragments can be estimated and subsequently used in order to balance (i) the truncation effects occurring at the outer boundary of the cluster and (ii) the charge transfer effects occurring between the fragments.

The first problem (i) has been solved by endowing each atom with a reference charge (q0), whereas the second issue (ii) has been faced by imposing that the total charge on each fragment (Q) depends on the geometry of the system. Noted is that the charge transfer effects have been included only in the subsystem treated with classical approaches, while no charge transfer has been taken into account across the QM and MM boundaries.

Energy and Charges

The variation of the energy of an isolated atom with respect to its value for a reference charge q0 is approximated by a Taylor series truncated at the second order:[12]The total electrostatic energy () of a system of N charges can be partitioned into self-energy () and interaction () contributions,[12,13] with the former term being just the sum of atomic energies defined according to eq and the latter term being an effective charge–charge interaction:The interaction potential J(R) is a function of the module of the interatomic distance , depending on some atomic parameters of the ith and jth atoms.

By substituting eq in the last equation, we obtain

Next, we impose that the functions J(R) satisfy the condition:employing the Ohno functional form.[14] Then expansion of (q – q0)2 leads to the following expression for the energy: is just the sum of the energies e0 related to the q charges and can be set to zero by choosing the proper reference for the energy. The functional in eq can be viewed as a generalization of the previously implemented FQ method,[13] which can be obtained in the limit of for all atoms.

The set of FQ charges q can be obtained by a constrained minimization of eq with respect to the charges, imposing the constraint of the charge conservation on each fragment (molecule). In order to do that, it is convenient to split the i index running over the atoms (i = 1,..., N), into two indices, one (I) running over the fragments (I = 1,..., Nfrag) and the other (ι) running over the atoms within the Ith fragment (ι = 1,..., N),so that a generic summation over i becomesThe charge conservation constraint on each fragment can be easily written with this formalism aswhere Q is the total charge on the fragment I. One possibility is that the Q’s are assigned a priori at the beginning of the calculation and then kept fixed (e.g., by imposing that all the molecules stay neutral or have a fixed predefined charge), but, more generally, the fragment charges can be a function of the geometry of the system (Q(R)) as we will show in the next sections.

The Lagrangian associated with is then defined asThe derivatives of with respect to the variational parameters (charges and Lagrange multipliers) are then:where the index K runs over the fragments and κ over the atoms belonging to each fragment. Imposing the minimum condition for , the resulting set of equations isThis problem can be recast by using the formalism of ref (13), so that the charges can be obtained by solvingwhere in the right-hand side the atomic electronegativities (χ), the total charge on each fragment (Q), and the vector obtained by the Hadamard product of the arrays containing the atomic hardness and reference charge (η ◦ 0) are included.

Whenever the FQ system interacts with a QM charge density, the electrostatic interaction between the charges and the QM charge density must be added to eq . In terms of D, the one-particle QM density matrix in the atomic orbital (AO) representation, the electrostatic potential V at the coordinates of the kth atom iswhere Z is the atomic number of the atom A of the QM subsystem, the μ,ν indices run over the atomic orbital basis functions , the σ index runs over the spin states (α and β), and the integrals are defined asThe FQ charges are obtained by solving

A Strength Function for Hydrogen Bonds

Now, restricting the treatment to the MM subsystem only, we want to define a function of the geometry of the system which is able to count the intermolecular interactions occurring among the FQ fragments, and estimate the ”strength” of such interactions.

Hence, indicating with {R⃗, k = 1, ..., N} the position vectors of the FQ atoms, for each pair of k and l atoms, the function f(R, ϑ) can be defined, where R = |R⃗ – R⃗| and is the valence angle involving a third atom. As detailed in the original work,[11] such a function of two variables is defined as a product of two functions of one variable, and gives values in the interval where the empirical parameters are R0 = 1.85 Å, σ = 0.20, ϑ = 10°, and σ = 10°.

The function provides a continuous estimate of the relative strength of hydrogen bonds, which vanishes, for a given pair of atoms, in the absence of any interaction and goes to the asymptotic value of 1 when the hydrogen bond reaches its maximum strength.

The dependence on ϑ can be better understood by looking at Figure , where a trimer of water molecules is shown.

Figure 1

Three hydrogen bonded water molecules.

Oxygen 1 is involved in two hydrogen bonds with atoms 5 and 9. For example, the strength of the first one depends on the distances R15 through A(R15) and on the position of atom 4—the first bonded neighbor of atom 5—through the function , so that the strength function is .

The functions are symmetric with respect to the permutation of k and l indices, and once they are evaluated for all the intermolecular atomic pairs, the result is a very sparse lower-triangular matrix (F) containing all the H-bond strengths. For the system in Figure , assuming that the two hydrogen bonds have strengths of, say, 0.8 and 0.7, the F matrix has the following nonvanishing elements: F(1,5) = 0.8 and F(1,9) = 0.7.

A Model for Outer Boundary Effects from the Strength Function Analysis

As shown in a previous work,[15] the FQ atoms located at the edge of the system suffer from truncated connectivity, resulting from nonphysical coordination numbers. Such a problem was solved by endowing the more external water molecules with fixed charges (FX) obtained from the TIP3P-FB model.[16]

While this approach was very effective to correctly reproduce the charge distribution of the FQ subsystem, it can become unpractical when it has to be applied to many snapshots, since it requires the manual identification of the atoms close to the outer boundary and the reassignment to them of the sought fixed charges. A more viable solution can be based on the observation that for systems whose primary intermolecular interactions are hydrogen bonds (e.g., aqueous solutions), the number of hydrogen bonds formed by each molecule can be used to define an effective and continuous coordination number and the strength function can be employed to automatically detect the atoms located at the edge of the molecular system which are affected by boundary effects. In particular, we defined a threshold value of 0.01 (1%) for the hydrogen bond function f(R,ϑ) above which a (possibly weak) hydrogen bond does occur. Then, for an oxygen atom belonging to a water molecule in the bulk phase, the ideal coordination number is 2, since two lone pairs are available for two intermolecular bonds with the hydrogen atoms of other two water molecules. Analogously, for a hydrogen atom belonging to a bulk water molecule, the ideal coordination number is 1, since it can in principle bind one oxygen from another water molecule. Then, if the coordination number for a given water molecule is smaller than the above values, we assume that such a molecule is located near the outer boundary, so that it can interact with an artificially reduced number of other solvent molecules and the electrostatic potential experienced by such a molecule is wrong. Hence, a correction must be included in order to restore the correct electrostatics.

We decided to face such an issue by increasing the absolute value of the q0 charge associated with the selected atoms as follows:where NH is a signed adimensional integer array containing the number of hydrogen bonds formed by each solvent atom and α is an empirical parameter that provides the correct physical dimensions. Noted is that the sign of the correction depends on whether the atom is the donor or the acceptor in the hydrogen bond under examination. For example, in the molecular system shown in Figure , the integer array NH has dimension 9 and nonvanishing elements NH(1) = 2, NH(5) = NH(9) = −1, indicating that atom 1 is involved in 2 hydrogen bonds as a donor, atoms 5 and 9 are involved in 1 hydrogen bond each as acceptors, while all the other atoms are not involved in any hydrogen bonding. Thus, all the three water molecules suffer from nonphysical coordination due to boundary effects, and the q0 charges correcting for such an effect are q0(2) = q0(3) = q0(6) = q0(8) = -α and q0(4) = q0(7) = 2α, whereas q0(1) = q0(5) = q0(9) = 0.0.

Thus, with such a model, the correction for the boundary effects can be automatically included in the calculations without any user modification of the input file, and this approach is particularly suitable when a large number of snapshots from a MD simulation need be analyzed.

A Charge Transfer Model from the Strength Function Analysis

We assume that, in a system ruled by intermolecular hydrogen bonds, the intermolecular charge transfer on one atom is proportional to the strength of all the hydrogen bonds occurring between that atom and the others. Assigning opposite signs when the atom behaves as donor or acceptor respectively, the total charge (Q) on the fragment J is proportional to the sum of all the atomic charge transfers coming from the atoms in that fragment (ι ∈ I) and all the others.

The hydrogen bond function can be used to define effective charges over the fragments (Q), to be included into the charge constraints of eq , thus providing an estimate of the charge transfer between the fragments. Since the original function takes values between 0 and 1, a simple approximation is to employ a single scaling factor to obtain fragment charges Q:where s is a sign function with values +1 and −1 for donor and acceptor atoms, respectively, whereas a is the scaling factor, which allows us to go from the adimensional values of f to the charges, providing also the correct physical dimensions. For the molecular system of Figure , the Q array issuggesting that an escaped charge proportional to +1.5 is found on atom 1, which is compensated by charges of −0.8 on atoms 5 and −0.7 on atom 9. Of course, the convention adopted for the signs enforces the conservation of the total charge of the system.

It is worth noting that such a scheme is able to take into account the CT among the FQ fragments, but it cannot account for any correction at the QM-FQ boundary, due to the intrinsic limitations of the QM/MM scheme.

QM/MM Energy

When the classical system is the low level part of a QM/MM calculation, the total energy is[4,13,15,17−20]where the functional form of EQM depends on the specific choice of the QM method employed for the treatment of the high level region, while EMM contains all the intra- and intermolecular terms of the MM force field of the low-level region. When the QM and MM subsystems are not covalently bonded and charge transfer effects between the two regions are neglected, the interaction energy term Eint can be partitioned into electrostatic, polarization, repulsion and dispersion contributions.[15] The first two contributions are treated by the approach sketched in the previous sections, whereas the latter two terms can be approximated by simple Lennard–Jones (LJ) functions, so that the QM/MM interaction energy can be written aswhere ELJ includes the dispersion-repulsion energy contributions between the QM and MM subsystems, while the second term is the electrostatic interaction between the QM and MM subsystems. In order to focus on the new developments of the FQ model, it is convenient to drop out from EMM the electrostatic interaction between the charges of the MM subsystem, so that Etot readswhere EMM* includes all the bonded (stretching, bending, torsion...) and nonbonded interactions between the MM atoms, with the exception of the electrostatic interaction.

Since the q charges do not enter the QM/MM electrostatic interaction energy, the Fock and coupling matrix elements related to eq are analogous to the ones of the FQ model previously implemented,[15,18] with the obvious modification that the FQ charges are obtained from eq .

Although not used in the present paper, we implemented also first derivatives with respect to nuclear coordinates. Since the q charges are onsite corrections, they do not introduce any change in the gradient formulas, but this is not the case for the Q vector. In fact, eq is variational with respect to q and λ but not with respect to the Q. As long as the Q are just fixed numbers, we can assume thatis identically zero, sinceHowever, this is no longer true when Q = Q(R) and a new term appears in the total gradient, which for the charge transfer contribution to the x-component of the force on a generic atom l, reads:The derivatives ∂/∂x are analytical and can be easily obtained from eqs and 19 (see the Appendix).

Molecular Dynamics Simulations and Spectra Computation

The TIP3P-FB[16] rigid water molecule was used in all the simulations, whereas the topology and initial GAFF2[21] parameters of creatinine have been obtained with the PrimaDORAC[22] web interface. The geometries of amine and imine tautomers were optimized at the B3LYP/6-31+G(d) level taking into account bulk solvent effects by means of the polarizable continuum model (PCM),[23] and point atomic charges were assigned performing a CM5 population analysis[24] of the corresponding Kohn–Sham orbitals, following the protocol adopted with success for the study of several other systems.[8,25−28] Molecular structures, atomic labels, and charges for the two tautomers are shown Figure .

Figure 2

Structure of the imine and amine tautomers of creatinine with labels and atomic charges of atoms involved in hydrogen bonds. A TIP3P-FB water molecule with its corresponding atomic charges is also shown for the sake of completeness.

All the simulations were performed employing the Proxima library[29] to build spherical boxes filled with water molecules at the experimental density and, in the case of creatinine aqueous solution, placing a creatinine molecule at the center of the box and deleting all molecules within the sum of covalent radii from any creatinine atom. In order to avoid any convergence issue, MD simulations were performed using large cavities with radii of 20 and 21 Å (containing 1116 and 1170 water molecules) for pure water and creatinine solution, respectively. Then, QM/MM computations of charges and, possibly, spectra were performed for the four polarizable models taking the water molecules occupying the innermost 12 Å of the spherical cavities employed in the MD simulations.

Following an established protocol,[30,31] we carried out MD simulations under spherical nonperiodic boundary conditions (NPBC) for both tautomers employing a locally modified version of the Gaussian[32] suite of programs. The RVV1 integrator[26] was used in all simulations, with a convergence criterion ϵ = 10–9 for the calculation of quaternion derivatives. The “rough walls”[26] boundary condition and the Bussi–Donadio–Parrinello[33] thermostat were used in all the NVT/NPBC simulations. The equilibration of the system involved an initial minimization with the conjugate gradient method and a subsequent simulation for 1000 ps with a small integration step of 0.5 fs and temperature of 298.15 K. The production run was then initiated at 298.15 K and continued for 25 ns with an integration time step of 4.0 fs. During the simulations the creatinine geometry was kept frozen. Snapshots were saved every 2 ps. Discarding an initial part of both simulations of 1 ns a total sampling of 9 ns was used for all the subsequent analyses.

The calculation of structural properties was carried out with standard procedures, taking into account proper normalization for NPBCs.

The selection of snapshots chosen in order to provide a well converged description of the trajectories was carried out with an in house implementation of the density peaks method.[34] To estimate the local density we used a Gaussian kernel with σ = 0.02N, where N is the number of points in the data set. As a feature space for the snapshot selection of the creatinine trajectories we used the FHB hydrogen bond definition as a function of time for the atoms Hr, Oc, Hd, Nd (imine form) and Hd1, Hd2, Oc, Nr (amine form), i.e., a four dimensional feature space.

For the selected snapshots, TD-DFT computations were performed to simulate the UV–vis spectra employing the B3LYP density functional in conjunction with the 6-31+G(d) basis set. Environmental effects were taken into account by QM/MM computations employing four different polarizable methods, which starting from the standard fluctuating charge approach (FQ), include charge-transfer contributions (FQ-CT) or outer boundary corrections (FQ-BC) or both terms (FQ-All).

Results and Discussion

The χ and η parameters entering eq have been optimized in ref (18) in order to reproduce the dipole moment of an isolated water molecule and the TIP3P-FB charges[16] for bulk water. The same results are obtained in all the FQ formulations for water molecules forming either 0 or 4 hydrogen bonds. The FQ-CT variant approaches the parent FQ model as β tends to zero restoring the original Q = 0 charge constraints (see eq ), whereas it leads to different results for intermediate numbers of hydrogen bonds when β ≠ 0, with β = 0.1 leading to charge transfers close to the reference CM5 values for the six clusters shown in Figure . For instance, the charge transfer estimated from CM5 charges on top of B3LYP/6-31+G(d) computations for the isolated water dimer is 0.09, which becomes 0.11 when including bulk solvent effects by means of PCM, with both values comparing well with the FQ-CT estimate of 0.1. Coming to the FQ-BC and FQ-All variants, α = 0.2 (see eq ) leads to charges close to the limiting TIP3P-FB values for solvent molecules near the outer boundary.

Figure 3

Cluster models with atom labeling.

For purposes of illustration, the oxygen atomic charges of a representative configuration issued from a full MD simulation of pure water (see Figure ) obtained employing the four different polarizable models (FQ, FQ-CT, FQ-BC, and FQ-All) and treating the central water molecule either at MM or QM levels are shown in Figure .

Figure 4

Snapshot of a water MD simulation (left) and atom labeling of the five central water molecules (right).

Figure 5

Difference (Δq°) between MM and QM/MM charges on water oxygen atoms.

Noted is that MM charges are intrinsically larger than the QM counterparts due to the need of reproducing dipole moments by a three-site point charge model. The results displayed in Figure show that, when present, the central QM molecule mainly perturbs the hydrogen bonded water molecules, while increasing the distance from the central oxygen atom (O1) the QM/MM and MM systems have almost the same atomic charge values.

From a quantitative point of view, the FQ-CT model provides charges slightly smaller, in absolute value, than the original FQ model without changing the spread between the smallest and largest value, whereas the FQ-BC model increases both the absolute values and the spread. Finally, the FQ-All model fulfills the sought requirement of decreasing the spread without reducing the smallest value of the FQ charges.

Analogous trends are apparent in Figure for the whole MM droplet.

Figure 6

Difference (Δq°) between polarizable and TIP3P-FB charges on water oxygen atoms.

Creatinine in Aqueous Solution

Hydrogen Bond Structure and Dynamics

The trajectories obtained from molecular dynamics simulations have been analyzed to characterize the interactions of imine and amine tautomers of creatinine with the aqueous environment. The solute–solvent interactions are essentially ruled by the formation of intermolecular hydrogen bonds involving the creatinine atoms indicated in Figure .

The first structural information related to the hydrogen bond interaction is retrieved by computing the pair radial distribution functions between creatinine atoms involved in hydrogen bonds with water molecules for both tautomers (see Figure ). Although the average interaction between the creatinine tautomers with solvent is similar to that reported in the study by Léon at al.,[35] the analysis of molecular dynamics trajectories provides not only structural, but also dynamical information.

Figure 7

Radial distribution functions related to the creatinine/water interactions. Full and dashed lines refer to radial distribution function, g(r), and coordination number, n(r), respectively. In red and blue are the results for imine and amine tautomers, respectively.

Comparison of the g(r) related to the interaction between the carbonyl oxygen of both tautomers with the hydrogen atoms of water (Oc···H in Figure ) shows that the amine tautomer has stronger interactions with neighboring solvent molecules and slightly higher coordination number (g(r) integral) than the imine tautomer. Also the interactions involving the same hydrogen atom in both tautomers (Hd for the imine and Hd1 for the amine, respectively) are stronger for the amine tautomer. In particular, the Hd atom is the only one which forms a very weak intramolecular hydrogen bond, whereas all the other atoms considered in Figure are involved in interactions with solvent molecules.

Further information on the strength and dynamics of the hydrogen bond interactions can be obtained by using the hydrogen bond function (FHB) proposed by Pagliai and co-workers[11,36] and adopted in the present work to determine the amount of charge transfer in both the FQ-CT and FQ-All models. The analysis has been performed separately for the two creatinine tautomers, considering the interactions involving the same atoms for which the radial distribution functions in Figure have been computed.

Figure shows the histogram of the FHB function for the imine tautomer atoms involved in hydrogen bond interactions. The analysis provides additional information with respect to the pair radial distribution functions. In all the histograms, the bar related to FHB = 0 has non-negligible values, indicating that the hydrogen bond does not involve the same molecules but processes of forming and breaking the interaction occur during the simulation. This picture is confirmed by the analysis of FHB as a function of simulation time shown in Figure , where the low average value is essentially due to the Hd atom which is only seldom involved in hydrogen bond interactions with solvent molecules (see Figure ).

Figure 8

Histogram of the FHB function values for the imine tautomer.

Figure 9

Simulation of the imine tautomer: value of the total FHB function summed for all solute atoms and corresponding histogram. This histogram was used to select frames for the subsequent analysis.

The FHB function has been analyzed also for the amine tautomer, showing a stronger overall solute–solvent interaction (see Figure ). This result suggests a higher stability of the amine tautomer in aqueous solution, in agreement with previous results.[35] These conclusions are corroborated by the comparison of (i) the time evolution of the FHB function and (ii) the total histograms of Figures and 11. In particular, this behavior can be related to the interactions involving the Oc and Hd1 atoms, which form stronger hydrogen bonds with solvent molecules.

Figure 10

Histogram of the FHB function values for the amine tautomer.

Figure 11

Simulation of the amine tautomer: value of the total FHB function summed for all solute atoms and corresponding histogram. This histogram was used to select frames for the subsequent analysis.

We used the hydrogen bond information also to select the snapshots to be used in the following spectra calculations. A preliminary principal component analysis (PCA)[37] of the FHB values was run to check for the correlation between the various sites and to evaluate the possibility of limiting the clustering to a 2D problem. The fraction of total variance recovered by using the first two eigenvectors for the imine and amine tautomer was 67.06% and 66.08%, respectively, showing limited correlation between the four atoms involved in hydrogen bonding in both trajectories. Given the limited amount of features, we discarded the idea of using a 2D space, determining, instead, the number of clusters and their centroids by inspection of a specialized plot (known as decision graph), which shows the local density of each point in the data set (ρ) on the x axis and the distance from the nearest high density point (δ) on the y axis. Cluster centroids are determined as points with a high local density and far away from other high density points i. e. as density peaks. These are the outlier points highlighted in the right top of the graph.

The decision graphs for both tautomers shown in Figure clearly indicate that four and three density peaks can be identified in the two trajectories; the sizes of the corresponding clusters (which can be used to estimate their statistical weight) are shown in Table . These clusters correspond to structures of the two tautomers with different values of FHB (see Figures –11).

Figure 12

Decision graph for the imine (left) and amine (right) forms. Cluster centroids are indicated by a yellow hexagon and by their frame label.

Table 1

Centroids and Sizes of Clusters Determined for Imine and Amine Simulations

	cluster number
	1	2	3	4
imine
centroid at frame	1030	1295	2083	2301
number of frames in the cluster	1068	943	1459	1531
amine
centroid at frame	1222	2233	2288
number of frames in the cluster	1891	1443	1667

For the imine tautomer (Figure ), it can be easily observed how the four clusters identify different coordination states for moieties such as the carbonyl oxygen or the Nd-Hd group that can form multiple hydrogen bonds and also the cooperative effect that favors the formation of single bonds for the carbonyl on the Nr-Hr side for all for clusters. An analogous effect can be observed for the amine tautomer in the formation of hydrogen bond networks with the atoms Hd2, Nr, and Ocb (see Figure ).

Figure 13

Cluster centroids for the imine tautomer. The structures corresponds to the simulation frames labeled in the decision graph (Figure , right panel). First neighbor water molecules of hydrogen bonding sites are shown with their donor (or acceptor) hydrogen distance in angstrom. The centroids are shown left to right following the order of Table .

Figure 14

Cluster centroids for the amine tautomer. The structures corresponds to the simulation frames labeled in the decision graph (Figure , right panel). First neighbor water molecules of hydrogen bonding sites are shown with their donor (or acceptor) hydrogen distance in angstrom. The centroids are shown left to right following the order of Table .

Polarization Effects

Figure shows the average charge of the oxygen atom of water molecules as a function of the distance from the creatinine center of mass for the different polarizable models. The analyzed system is made up by 1 QM solute (the imine tautomer of creatinine) and 250 TIP3P-FB water molecules, whose coordinates have been extracted from the molecular dynamics trajectory. Atomic charges have been subsequently obtained a posteriori through QM/MM calculations by modeling water molecules with FQ, FQ-BC, FQ-CT, or FQ-All approaches.

Figure 15

Atomic charges on water oxygen atoms issued from different polarizable models.

While FQ charges have been successfully applied to determine spectroscopic properties of QM solutes in aqueous solution, they have nonphysical behavior near the outer boundary region, as well evidenced by the significant depolarization occurring for distances longer than 10 Å (see Figure ). This behavior does not occur by using the FQ-BC approach, which solves the problem of an incorrect number of hydrogen bond interactions for those water molecules located in the boundary region, thus leading to atomic charges similar to those of the TIP3P-FB model (see Figure ).

Concerning the charge transfer model, Figure shows that the water oxygen charges issuing from the FQ-CT and FQ-All models are quite similar to those provided by the FQ and FQ-BC counterparts, respectively. This is expected for the system under study, where the fragments are relatively weakly bonded (hydrogen bonds only) and each fragment behaves at the same time as donor and acceptor, leading to compensating charge flows in opposite directions.

As a matter of fact, the fragment charges are between −0.05 and 0.05 for more than half the water molecules and between −0.1 and 0.1 for nearly 90% of the water molecules (see Figure ).

Figure 16

Histogram of total charge on water molecules issued from the FQ-All model.

Electronic Absorption Spectra Simulation

The clustering procedure discussed in a previous section permits a strong reduction of the number of snapshots for which electronic spectra must be explicitly computed. Several tests showed that, together with cluster centroids, very few additional snapshots need be considered for obtaining results virtually indistinguishable from the reference results issued from 200 equally spaced snapshots. Although satisfactory results are obtained by adding just a single structure (the most distant from the centroid) for each cluster, we preferred to include about 10 snapshots for each cluster, with the precise number being proportional to its size (given in Table ). Fully converged results are obtained when employing 50 frames for each tautomer, and, in particular 11,9,15,15 frames for the 4 clusters of the imine tautomer and 19,14,17 frames for the 3 clusters of the amine tautomer.

The UV–vis spectra have been computed for all those frames (and then averaged) with QM/MM approaches involving fixed or different flavors of polarizable charges.

From the comparison of panels (a) and (b) of Figure , it is apparent that the simulated UV–vis spectra for the two tautomers show substantial differences. In particular, only the UV–vis spectrum of the amine tautomer is in good agreement with that observed experimentally,[38] confirming that this species is the one present in aqueous solution. A quantitative estimate of the agreement with the experimental data can be obtained by the comparison of the wavelengths experimentally observed[38] for the two bands that characterize the UV–vis spectrum of creatinine in aqueous solution with the counterparts computed with the different polarizable approaches.

Figure 17

UV–vis spectra of creatinine tautomers: (a) imine and (b) amine tautomers. TD-DFT calculations have been carried out at the B3LYP/6-31+G(d) level for the QM part. Solvent molecules have been described with the FQ, FQ-CT, FQ-BC, and FQ-All polarizable models and with TIP3P-FB fixed charges (FX).

The results collected in Table show that the absorption maxima of amine and imine tautomers are similar in the gas phase and also employing the PCM to take into account bulk solvent effects. On the other hand, all the models based on fluctuating charges displace by more than 10 nm the red-side maximum of both tautomers and the blue-side maximum of the amine tautomer but have a negligible effect on the maximum near 190 nm of the imine tautomer.

Table 2

Wavelength (in nm) of Absorption Maxima and Relative Intensity of the Less Intense Band in the UV–Vis Spectra of Imine and Amine Tautomers of Creatininea

	imine tautomer	amine tautomer
gas B3	189, 226 (0.39)	192, 224 (0.64)
gas B2	187, 208	193, 216
PCM	192, 232 (0.31)	190, 227 (0.50)
FX	182, 242 (0.17)	200, 243 (0.39)
FQ	188, 251 (0.14)	213, 255 (0.33)
FQ-BC	187, 249 (0.16)	207, 249 (0.43)
FQ-CT	185, 249 (0.18)	203, 247 (0.52)
FQ-All	185, 250 (0.17)	204, 247 (0.50)
best	183, 232 (0.17)	205, 239 (0.50)
expt[38]		201, 233 (0.53)

The TD-DFT calculations have been carried out at the B3LYP/6-31+G(d) (B3) level for the QM part, while the solvent has been described by the PCM, TIP3P-FB fixed charges (FX), and FQ, FQ-CT, FQ-BC, and FQ-All models. For the isolated molecule (gas), TD-DFT band maxima computed at the B2PLYP/jun-cc-pVTZ level (B2) are also reported and used, together with B3 relative intensities and QM/MM solvent shifts to obtain the best estimate values (best).

Furthermore, the relative intensity of the red-side peak of the imine tautomer is much lower (both in the gas-phase and in aqueous solution) than that of the amine counterpart, with the latter being in remarkable agreement with experiment. Both inclusion of CT and correction for spurious outer boundary effects lead to results intermediate between those issued from fixed charges or bare fluctuating charges, with the role of BC contributions decreasing after introduction of CT. Although the computed solvatochromic shifts of the spectra are not strongly sensitive to small modifications of solvent charges, it is gratifying that the FQ-CT and FQ-All models improve the agreement with experiment, thus suggesting that charge-transfer between solvent molecules could play a non-negligible role in tuning the conventional FQ results, which possibly overcorrect the FX counterparts. The close similarity between the FQ-CT and FQ-All models is related to the marginal effect of water molecules near the outside boundary on the solute UV absorption.

Computation of the spectrum of the isolated molecule at a higher level of theory (B2PLYP in conjunction with the jun-cc-pVTZ basis set[39] at the CCSD(T)-F12/cc-pVDZ-F12 geometries reported in ref (35)) has a negligible effect on the blue-side maximum but significantly improves the agreement between theory and experiment for the red-sideband. At this level, the difference between theory and experiment is comparable for both bands (4 and 6 nm, respectively) and also their relative heights are close to the experimental counterparts. It is noteworthy that the improvement related to more accurate gas-phase computations is comparable to that issued from more refined polarizable models of solvent effects, so that both contributions need to be taken into account for a balanced computation of spectroscopic properties in solution. In this connection, we point out that for broad bands like those characterizing the UV–vis spectrum of creatinine in aqueous solution, the role of vibronic contributions (which could be effectively taken into account by the methodology employed in previous studies[8,9]) is negligible. In summary, a physically sound description of solvent effects rooted into the fluctuating charge model coupled to TD-DFT models employing hybrid and double-hybrid functionals provides remarkably accurate results with nonprohibitive computer times for the UV–vis spectrum of a medium-size chromophore in aqueous solution, allowing also for an unambiguous discrimination between different tautomeric forms.

Conclusions

In this paper, we have extended the fluctuating charge polarizable model to take into account charge transfer and to rectify the unbalanced polarization induced by outer boundary effects. The modifications with respect to the starting model are very simple and can be added straightforwardly to any code implementing the FQ model. A first parametrization for water has shown a remarkable robustness and accuracy for reproducing solvent characteristics. Next, the UV–vis spectra of two creatinine tautomers have been computed at the TD-DFT level by means of different polarizable QM/MM models. In particular, the FQ-All model, which takes into account charge transfer between solvent molecules and enforces correct coordination numbers at the outer boundary, leads to satisfactory results, thus paving the route toward its general use in conjunction with nonperiodic boundary conditions for the study of spectroscopic phenomena in condensed phase. Work is in progress in our laboratory to extend the parametrization of the model to other widely used solvents, e.g., alcohols or CH3CN.