Electronic absorption spectra of pyridine and nicotine in aqueous solution with a combined molecular dynamics and polarizable QM/MM approach.

The electronic absorption spectra of pyridine and nicotine in aqueous solution have been computed using a multistep approach. The computational protocol consists in studying the solute solvation with accurate molecular dynamics simulations, characterizing the hydrogen bond interactions, and calculating electronic transitions for a series of configurations extracted from the molecular dynamics trajectories with a polarizable QM/MM scheme based on the fluctuating charge model. Molecular dynamics simulations and electronic transition calculations have been performed on both pyridine and nicotine. Furthermore, the contributions of solute vibrational effect on electronic absorption spectra have been taken into account in the so called vertical gradient approximation.

Introduction

The determination of spectroscopic properties of molecules in aqueous solution requires a detailed understanding of solute/solvent interactions, particularly when complex effects, such as hydrogen bonding, play a significant role. To obtain useful information, accurate models describing solute/solvent interactions are needed. Although ab initio molecular dynamics (MD)1 or quantum‐mechanics/molecular mechanics (QM/MM) approaches2, 3 have been revealed to be effective in the description of a series of hydrogen bonded systems, their application has been limited by computational cost in considering systems characterized by an elevated number of atoms and/or by relative dynamics that occur with relatively long time scale (typically, larger than 10 ps). In such cases, classical methods, in particular MD simulations, can provide the required amount of phase space sampling. The intra‐ and inter‐molecular interactions in classical methods are modeled through force fields like, for example, AMBER,4, 5, 6, 7, 8, 9 CHARMM,10, 11, 12, 13 GROMOS,14, 15, 16, 17, 18, 19, 20, 21 and OPLS,22, 23, 24, 25, 26 which are continuously developed to improve the accuracy and to extend the application to different systems. One of the limitations of classical methods is represented by the impossibility to have access to spectroscopic properties, which require the knowledge of the electronic structure of the studied systems. To overcome this problem, it is possible to adopt a number of computational strategies. For example, the calculations of the sought properties can be obtained performing QM calculations on the solute and a limited number of hydrogen bonded solvent molecules and eventually considering the average solvent effects with a polarizable continuum model (PCM)27, 28, 29 or similar methods.30, 31 This approach allows to avoid problems related to the description of the intermolecular interactions, which perturb the electronic structure and consequently the spectroscopic properties of the solute, but it fails, for example, whenever the solute can populate different conformations or it can form stable hydrogen bond interactions with a different number of solvent molecules. A different approach relies on the adoption of QM/MM methods,2, 3 which can suffer however from some polarization problem when the solvent molecules are treated with fixed charges (FX).32 In this respect, many different methods have been developed in the last years to include polarization effects in the total energy of QM/MM schemes,33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44 allowing a mutual polarization between the QM charge density and the MM charge distribution. Among these, we have selected a recent QM/MM implementation including both fluctuating (FQ) and fixed point‐charge distributions44 and based on the electronegativity equalization principle.45, 46, 47

The UV‐vis spectra of nicotine in water have been calculated with this QM/MM approach for a series of configurations extracted from a MD trajectory. Nicotine molecule represents a challenging system, which has been adopted as a test case in gas phase48 to verify the accuracy of a series of computational method to determine electronic excitations. The results show that TD‐DFT calculations provide results similar to multireference and coupled cluster approaches and allow to obtain first insights on the solvent effects on the electronic absorption spectrum, simulating nicotine in aqueous solution with PCM and QM/MM methods.49

To perform accurate simulations including proper sampling of the conformational freedom of nicotine related to the dihedral angle (see Fig. 1), a new force field specific for this system has been developed. The approach used to parametrize the electrostatic part and, in particular, the anisotropy in the interactions due to the lone pair on the nitrogen atom, has been validated by MD simulations of the rigid aromatic moiety of nicotine (i.e., pyridine) in water. The simulations on the heterocycle have been performed using a force field developed with the same approach adopted for nicotine, explicitly describing the lone pair of the nitrogen atom using a virtual site (VS). Pyridine has been also adopted as test case to calculate electronic spectra in aqueous solution through a polarizable QM/MM approach,44 considering the vibrational effects in the harmonic approximation with the vertical gradient (VG) method.50, 51 The computational approach has been subsequently applied to the study of nicotine.

Figure 1

Molecular structure of pyridine (left) and the nicotine A (middle) and B (right) conformers. The spheres in cyan represent the position of the virtual sites. The distance is indicated as 0.38 or 0.39 Å depending on the algorithm adopted. The most probable dihedral angle value ( ) is reported for the A and B conformers, respectively. [Color figure can be viewed at wileyonlinelibrary.com]

Methodology

Nicotine force field

The molecular structure of the nicotine molecule, shown in Figure 1, consists of two building blocks, namely a pyridine rigid ring and a methylpyrrolidine five‐membered ring. The dihedral angle around the bond connecting the two aforementioned moieties ( ) gives rise to a number of different conformations, which are expected to be populated at room temperature.

The force field used in the MD simulations has been developed performing a series of DFT calculations at B3LYP/6‐31 + G(d) level of theory considering implicitly the interactions with water through the Conductor‐like Polarizable Continuum Model (C‐PCM)52, 53 to accurately reproduce the intramolecular features, with particular regard to the distribution of the dihedral angle.

It has been shown49, 54 that nicotine presents two different conformers (labeled A and B) characterized by two different values of the dihedral angle , as reported in Figure 1. The potential energy curve (PES) along the dihedral angle for a series of optimized molecular structures allows to identify the most stable conformers of nicotine, verifying that the conformer A is more stable than B. Consequently the conformer A has been adopted as reference in the development of the force field.

To properly describe Coulomb interactions, a population analysis based on Charge Model 5 (CM5)55 has been carried out. It has been recently shown by Jorgensen and coworkers,56, 57 that to take into account polarization effects the atomic charges have to be increased with an appropriate scaling factor to obtain a better agreement with experiments. To alleviate the solute polarization problem, a different approach has been adopted, which consists in the determination of CM5 charges performing DFT calculations at B3LYP/6‐31 + G(d) level of theory, describing the solvent effects with the C‐PCM approach52, 53 imposing the value of the scaling factor for the sphere radius ( ) to 1.05. All DFT calculations have been performed with the Gaussian suite of programs.58 To further improve the description of the interactions between the heteroatom of the six‐membered ring with the environment, a VS has been introduced in the model to properly describe the lone‐pair directional character on the sp nitrogen atom. The VS position has been obtained through a localization of the molecular orbitals with the Boys method.59, 60 The centroid of the molecular orbital describing the lone pair of nitrogen atom has been adopted as the VS position, as shown in Figure 1. The centroids of localized orbitals have been calculated with the approach described by Vidossich et al.61 This procedure has been successfully applied in the study of structural and dynamic properties of formamide and its two N‐methyl derivatives, N‐methylformamide, and N,N‐dimethylformamide in liquid phase.62

The procedure adopted to assign the charge to the VS can be easily explained in the case of the pyridine molecule and it consists in the following step:

the charge on nitrogen atom is moved on the VS;

the charge on the VS is reduced to avoid overpolarization effects imposing the constraint that the dipole moment of the C‐VS‐C and C‐N‐C groups has to be equal.

This correction allows to obtain a very similar dipole moment for pyridine with and without VS. It has been verified that adopting the localization algorithm proposed by Pipek and Mezey,63 the position of the centroid related to the nitrogen atom lone‐pair is only slightly different, leaving the description of the hydrogen bond interactions unchanged.

The atomic charges, dipole moment and of pyridine molecule are summarized in Table 1.

Table 1

Atomic charges, dipole moment, and nitrogen VS distance for the models without, NO_VS, and with virtual site, VS, respectively.

Atom	NO_VS	VS
N	−0.424	0.000
C1	0.042	−0.033
C2	−0.095	−0.095
C3	−0.075	−0.075
H1	0.121	0.121
H2	0.120	0.120
H3	0.124	0.124
VS	0.000	−0.275
	NO_VS	VS_B
μ/D	2.812	2.812
r_N⋯VS/Å	–	0.38

The atom labeling proposed by Jorgensen et al.64 has been used

The nitrogen atom on the five‐membered ring does not require a VS as a consequence of the tetrahedral local structure due to the sp hybridization.

A set of 36 relaxed nicotine structures, with spanning the [− , ] interval (one configuration every ), and their relative energies have been adopted to determine the torsional potential related to the variation of the dihedral angle through a fitting procedure using the Joyce program.54 As it can be observed in Figure 2, the PES along the coordinate is characterized, as expected, by the presence of two different minima separated by two 25 kJ/mol energy barriers.

Figure 2

Solid black line: Torsional potential obtained by the fitting procedure with Joyce program.54 Red open circles: Potential energy differences computed for each relaxed structure of nicotine varying the value of the dihedral angle . The potential energy of conformer A has been taken as reference in the fitting procedure. [Color figure can be viewed at wileyonlinelibrary.com]

The parameters defining the force field of nicotine and of the adopted water model are listed in Table 2.

Table 2

Force field specifications, following the atom numbers and names in inset figure.

MD simulations

MD simulations of nicotine in aqueous solution have been performed with the GROMACS software,65 using the new developed force field. The procedure adopted to set up the electrostatic part of the force field has been validated performing a series of MD simulations on pyridine. In fact, due to the rigidity of the aromatic ring, only intermolecular parameters are critical for pyridine. OPLS‐AA parameters developed by Jorgensen et al.64 have been adopted to model pyridine, with the exception of the atomic charges. The water molecules have been described through the three‐site, fixed‐geometry model proposed by Wang et al.66 (TIP3P‐FB), which improves the determination of structural and dynamic properties with respect to the original TIP3P model.67 MD simulations have been carried out in the NVT ensemble, adopting the algorithm proposed by Bussi et al.68 The simulation box, using periodic boundary conditions, is made up of one solute molecule together with 512 (cubic box with side length 24.423 Å) or 3196 (orthorhombic box with Å, Å, Å) solvent molecules for pyridine and nicotine, respectively. The time‐step has been set to 0.2 fs, for a total accumulation run of 5 and 100 ns for the simulations of pyridine and nicotine, respectively. The particle mesh Ewald69, 70 method has been used to compute long‐range interactions.

Hydrogen bond characterization

The hydrogen bond interactions with solvent during the simulations of both nicotine and pyridine have been analyzed by means of the function introduced by Pagliai et al.71 The function for water molecule interacting with the nitrogen atoms of pyridine and nicotine is defined as: with and given by:

The values of the parameters are obtained from the histograms of the H‐bond distance, (r), and angle, ( ). represent the position of the first peak in (r) and ( ), respectively, while are the half widths at half maximum (HWHM) in (r) and ( ), respectively. are the instantaneous distance and angle involved in the interaction between the water molecule and the nitrogen atoms of pyridine and nicotine. The function assumes values in the range between 1 and 0, depending of the deviation of from the reference values .

The TD‐DFT/FQ/FX approach for UV‐vis spectra

The TD‐DFT/FQ/FX approach for the calculation of vertical excitation energies has been recently presented by Carnimeo et al.44 and is briefly summarized in the following. The total QM/MM ground state energy for a molecular system partitioned into a QM subsystem, a MM subsystem endowed with FQ charges, and a MM subsystem endowed with FX charges, can be written as where is the total energy of the QM subsystem, which depends on the particular level of theory used in the calculations, is the total energy of the MM force field, including—among the other terms—also the electrostatic Coulomb interaction between the FQ and FX charge distributions, is a parametrization of the dispersion‐repulsion interactions between the QM and the MM subsystems, while ( ) is the electrostatic interaction energy between the QM and the FQ (FX) charge distributions.

For a MM system composed by polarizable atoms and FX atoms, and when the QM charge density distribution is expanded in atomic orbital basis set reads where the ( ) index runs over the FQ (QM) atoms, is the atomic number of the atom , is the element of the AO ground state density matrix, is the array containing the polarizable charges, and the integrals are defined as

The and refer to the coordinates of the atoms and , respectively.

The energy term corresponding to the FX charge distribution is analogous to eq. (3), where the is replaced by the array containing the FX, and the proper atomic coordinates are used.

The FQ charges are obtained at each step of the SCF procedure by solving the linear system42 where is the charge‐charge interaction matrix, based on a semi‐empirical kernel72 using the atomic hardnesses as parameters and also including the proper terms for the charge conservation constraints; is the atomic electronegativity array, is the vector of the electrostatic potential generated by the QM and FX charge distributions and evaluated at the coordinates of the FQ atoms, and is the array containing the total charge on each FQ fragment, which defines the charge conservation constraints.

Then, at each step of the SCF, the FQ charges are polarized by the QM and FX charge densities through the external potential , as well as the QM charge density is polarized by the FQ charges by adding to the total Fock matrix the contribution related to eq. (3) ( )

Once the SCF is at convergence, the TD‐DFT equations are solved using standard numerical approaches,73, 74, 75, 76, 77 the contribution from eq. (3) ( ) being explicitely included in the coupling matrix elements

Starting from the MD trajectories of both pyridine and nicotine aqueous solutions, 100 snapshots equispaced in time have been extracted and used as input configurations for the subsequent excited state calculations. Although the MD simulations were carried out with periodic boundary conditions using orthorhombic boxes, only the water molecules inside a sphere with origin at the center of mass of the solute and a radius of 11 Å and 17 Å (in the case of pyridine and nicotine, respectively), were retained for the subsequent TD‐DFT/FQ/FX calculations. The FQ parameters used in this study (summarized in Table 3) have been taken from a previous work,44 where a new parametrization was proposed to reproduce the changes in the polarization of the water molecules for different environments, ranging from the isolated molecules in the gas phase to the bulk phase of condensed systems. Regarding the excited state calculations of hydrogen bonded systems, it was found that the vertical excitation energies were better reproduced by the TD‐DFT/FQ/FX model by including a proper screening in the electronegativity of the hydrogen atom directly involved in an hydrogen bond with a QM donor. Thus, for every water molecule in the FQ subsystem, the magnitude of the hydrogen bond with the QM solute has been estimated through the function [eq. (1)].71 The molecules showing values lower than a fixed threshold ( ) were considered as bulk water molecules, and the hydrogen and oxygen parameters in Table 3 were used; whereas, the molecules showing larger values where considered as directly interacting with the QM charge distribution, and the hydrogen/QM parameters were used for the hydrogen atom accepting the hydrogen bond.

Table 3

Parameters (in atomic units) for water molecules in the QM/MM calculations.

	Oxygen	Hydrogen	Hydrogen/QM
χ	0.189194	0.012767	0.042000
η	0.623700	0.637512	0.637512
q	−0.659	0.329	0.329

Electronegativity, hardness, and atomic charges are labeled with , and , respectively. The parameters of the FQ hydrogen atom directly bonded to the QM nitrogen are labeled as hydrogen/QM.

In the same work,44 it was proposed that a good strategy for a reliable representation of the MM charge density is to endow the outer layer of water molecules with fixed (FX)—rather than polarizable—charges as the electrostatic potential acting on such molecules is frustrated by the truncated coordination, and leads to nonphysical polarization effects. The solvation spheres end with a 2 Å layer made up by fixed‐charge TIP3P‐FB66 water molecules as shown in Figure 3 for both pyridine and nicotine systems.

Figure 3

One snapshot taken from the MD simulation of pyridine in water. The oxygen atom of the water molecule directly bonded to the QM is colored in red. The pink color refers to the FQ water molecules, whereas the green color has been adopted for the external layer of FX66 molecules. [Color figure can be viewed at wileyonlinelibrary.com]

The TD‐DFT/FQ/FX results were compared with TD‐DFT calculations in vacuo and TD‐DFT/PCM calculations using the C‐PCM52, 53 scheme, for the two molecules at optimized geometry, and in all cases the B3LYP exchange and correlation functional and the 6‐31 + G(d) basis set were used, including 12 and 24 states for pyridine and nicotine, respectively.

Eventually, for each snapshot and each transition , energies and oscillator strengths were convoluted in the energy domain with broadening functions , which in principle depend on the particular transition, and with a custom HWHM . The final electronic spectrum was then obtained by averaging the signals originated by each snapshot as and is reported in the wavelength domain, for a better comparison with experimental data expressed in nm.

A more refined approach is based on the fact that the sampling rate adopted in the MD simulations is much larger than the typical time step of the vibrational motions, so that the vibronic broadening effects can be evaluated separately from the solvent effects. First, the vibronic stick spectra were obtained within the VG approach (see Refs. [ 78, 79] and references therein.) by evaluating excited state potential energy surfaces in the harmonic approximation from ground state vibrational frequencies, normal modes and excited state forces. This choice was supported by the fact that a reliable approach for the calculation of the polarizable QM/MM excited state gradients has been recently developed and validated,44 so that the VG method is currently more appealing than other approaches80, 81 based on Hessian calculations (e.g., adiabatic Hessian) requiring further steps of parametrization and validation.

Vibronic calculations have been performed on pyridine or nicotine, interacting with one or two explicit water molecules, on the geometry optimized with the PCM. Then, the stick vibronic spectra of the relevant transitions were normalized and convoluted with Gaussian functions with a very small HWHM (20 cm ), and such spectral line‐shapes were used to define the functions. The solvent effects were then added using eq. (8) to convolute the vibronic line‐shapes with the excitation energies and oscillator strengths obtained along the MD trajectories.

Results

Pyridine in aqueous solution interacts with water essentially by forming hydrogen bonds, which involve the lone pair on the nitrogen atom of the heterocyclic compound with solvent.82, 83 To validate the proposed model, in which the lone pair on nitrogen atom is explicitly described through the VS, the H‐bond interactions have been compared with a model without VS (hereafter NO_VS) and previous experimental82, 84, 85, 86, 87, 88, 89, 90, 91 and computational92, 93, 94, 95, 96 results. The first information on the hydrogen bond structure has been obtained by determining the radial pair distribution function, g(r), for the N H and N O contacts, reported in Figure 4.

Figure 4

Radial and angular distribution functions for NO_VS (red line) and VS (black line). a) Radial distribution function (solid line) and running integration number (dashed line) for the N H intermolecular interactions. b) Radial distribution functions for the N O contacts. C: Angular distribution functions of the HO N angle, limited to the interactions with N H and N O distances within the first minimum position in the respective g(r). [Color figure can be viewed at wileyonlinelibrary.com]

The position of the first g(r) maximum obtained by both models is in agreement with the most recent experimental measurements by Kameda et al.97 (1.9 Å), with the exception of the coordination number, which is definitively lower. In fact the coordination number obtained by a fitting procedure of experimental data provides a value of 2.5, while the coordination number extracted from g(r) is 1.3 and 1.8 for VS and NO_VS, respectively. However, Bakó et al.83 have observed that an accurate estimation of both interaction distances and coordination number for the N O contacts is very difficult as a consequence of the concomitant presence of the inter‐molecular H‐bond (1.8–1.9 Å) between solvent molecules. Lower values for the number of hydrogen bonds between the heterocyclic nitrogen atom and the hydrogen atom of solvent (in the range 1–2) are obtained by both experiments82, 84, 85, 86, 87, 88, 89, 90, 91 and models to interpret experimental findings,92, 93, 94, 95, 96 with a preponderance of 1:1 complexes.82, 95

It is expected that the inclusion of the VS allows to obtain a more directional character in the H‐bond interaction. This behavior can be observed by calculating the angular distribution functions related to the N angle. Figure 3 shows the angular distribution functions limited to the N H and N O interactions with distances shorter than those of the corresponding first minimum. The function obtained with VS goes faster to zero than in the case of simulations with pyridine without VS, confirming the hypothesis on the VS effect. In all cases, the angular distribution functions reach a zero value for angle higher than 30 usually adopted as threshold (in conjunction with other parameters, such as for example the potential energy93, 94) to ascertain the presence of H‐bond interactions.

Further details about the number of hydrogen bonds between pyridine and water can be accessed determining the value of during the simulations. On the basis of experimental results and on the initial analysis of g(r) for the N H contacts, it is expected that the instantaneous value of is close to one. The analysis reported in Figure 5 confirms the expected value, especially in the case of the model with VS. The results are summarized in the histogram of the distribution of hydrogen bonds and are compared with those obtained adopting geometrical criteria, based on the position of first minimum in g(r) and on the zero value in g( ) shown in Figure 4.

Figure 5

Hydrogen bond analysis for the NO_VS (a and b panels) and VS (c and d panels) models. a and c panels: Evolution of . b and d panels: Blue line represents the probability distribution of ; Red bars represent the coordination numbers valuated with geometrical criteria. [Color figure can be viewed at wileyonlinelibrary.com]

The most probable coordination number computed during the simulations of pyridine with the VS model is 1, independently of the criterium adopted to calculate this quantity and is in agreement with the values extrapolated from the g(r). The main complex observed in aqueous solution is represented by the solute molecule bound to one water molecule as hypothesized to interpret a series of spectroscopic results.82, 95, 98 The NO_VS model presents instead a higher coordination number with a non‐negligible number of configurations with pyridine involved in two hydrogen bonds.

The correct representation of the H‐bond interactions with the VS model permits to adopt the present approach to model also the pyridine ring in nicotine moiety.

MD simulations on nicotine system have been carried out with the developed force field, analyzing initially the variation of the dihedral angle during the 100 ns accumulation run and the relative distribution. As observed in previous MD simulations on nicotine both in water49 and methanol,54 the most populated configurations are essentially related to the A and B conformers in agreement with the DFT computed PES, as Figure 6 shows. Results of experimental measurements have been correctly interpreted on the basis of models with a similar relative population of the two conformers.99, 100

Figure 6

Left panel: Evolution of the dihedral angle during the molecular dynamics simulations. Right panel: Histogram of the probability distribution of obtained by the analysis of the molecular dynamics trajectory.

The study of nicotine solvation in water requires a particular attention due to the different conformations accessible and to the presence of a nitrogen atom on both the pyridine and methylpyrrolidine rings, labeled N6 and N5 in the following. Preliminary information regarding the interaction with solvent has been obtained by analyzing the radial and angular distribution functions involving each nitrogen atom. Figure 7 shows that the position of the first peak maximum in the g(r) for the N H and N O interactions is nearly the same for the two different nitrogen atoms, whereas the stability of the hydrogen bond appears to be higher for N5 than for N6, as confirmed by the depth of the first minimum in the g(r), which indicates a lower rate of exchange between molecules of the first and second shell. This behavior can be explained on the basis of the higher basicity of N5 with respect to N6.49, 101, 102 The coordination number is closer to one for N5 than for N6.

Figure 7

a) Radial distribution functions for the N H contacts (solid lines) and running integration number (dashed lines). b) Radial distribution functions for the N O contacts (solid lines) and running integration number (dashed lines) c) Angular distribution functions of the HO N angle, only for N H and N O distances lower than the position of the first minimum in the respective . Blue and red lines refer to nitrogen atoms N5 and N6, respectively. [Color figure can be viewed at wileyonlinelibrary.com]

To characterize the hydrogen bond interactions the angular distribution functions have been computed for all the N H and N O contacts with distances shorter than that of the first minimum. The results are reported for the two nitrogen atoms in Figure 7, showing a slightly faster decay for the H‐bond with N5. To confirm the hypothesis that the main difference of the H‐bond interaction between water and the two nitrogen atoms on nicotine is essentially due to the rate of exchange between molecules of first and second solvation shells, the number of times that the hydrogen bonded molecule differs in two consecutive steps has been valuated. The exchange rate is 4 times higher for N6 than for N5.

Electronic absorption spectra of pyridine and nicotine in water

In Table 4, the vertical excitation wavelengths and oscillator strengths for the first six transitions of the pyridine molecule at B3LYP/6‐31 + G(d) level of theory (in the range 190–300 nm) are reported, as computed with TD‐DFT calculations on the optimized molecule in vacuo (VAC) and in PCM (PCM) and TD‐DFT/PCM calculations of the pyridine optimized with one explicit QM water molecule directly interacting with the nitrogen (PCM + H O). The results from TD‐DFT/FQ calculations on the QM pyridine directly interacting with one FQ water molecule—FQ(H O)—taken from a previous work44 have been also reported. The electronic absorption spectrum of pyridine has been studied with several computational methods,103, 104, 105, 106, 107, 108, 109, 110 which have been shown difficulties in the correct determination of electronic transitions energies using different levels of theory. To verify the effects of exchange and correlation functional or basis set, the calculations reported in Table S1 and Table S2, respectively, of the Supporting Information have been performed. The observed variations on the excitation energies and oscillator strengths are very small and do not change the overall description of the computed electronic spectrum of pyridine obtained using the B3LYP exchange and correlation functional in conjunction with the 6‐31 + G(d) basis set.

Table 4

Wavelength (nm) and oscillator strengths (f) for the first six excited state of pyridine in the gas phase (VAC) and in aqueous solution.

			Pyridine (from optimized geometry)
			VAC		PCM		PCM+H _2O		FQ(H _2O)
Label	Symm. (C _2v)	Assign.	nm	f	nm	f	nm	f	nm	f
1	B _1	n→π^*	256	0.0043	247	0.0058	239	0.0045	241	0.0047
D	A _2	n→π^* (Dark)	242	0.0000	231	0.0000	221	0.0000	226	0.0000
B	B _2	π→π^* (Bright)	225	0.0388	227	0.0666	227	0.0721	225	0.0433
4	A _1	π→π^*	197	0.0162	199	0.0199	198	0.0150	197	0.0127
5th	–	–	191	0.0027	184	0.0000	194	0.0000	184	0.0001
6th	–	–	188	0.0000	183	0.0739	183	0.0001	177	0.0002

Symmetry labels have been taken from gas phase calculations of pyridine constrained at the C point group symmetry. PCM and PCM + H O refer to TD‐DFT calculations on the optimized pyridine in PCM with and without one explicit hydrogen bonded QM water molecule. FQ(H O) refers to TD‐DFT/FQ calculations on the optimized pyridine treated at the QM level, and one explicit water molecule at the FQ level (from Ref. 44). M1 and M2 refer to the Mixed states found in the MD simulations. FQ(MD) and PCM(MD) refer to wavelengths and oscillator strengths averaged over the snapshots of the MD simulation, using the TD‐DFT/FQ/FX and the TD‐DFT/PCM models, respectively. The averages including the Bright/Dark states analysis are reported as FQ (MD) and PCM (MD) (see text). Standard deviations are given as subscripts.

Excited state calculations on optimized geometries

The features of the calculated UV‐vis spectrum of pyridine in vacuo are due to one weak transition at 256 nm well separated from all the others, one dark transition at 242 nm, and two bright transitions at 225 and 197 nm. When the PCM is used for modeling bulk solvation effects, with both the PCM and PCM + H O schemes, the wavelength of the bright transitions at 225 nm is quite unaffected by the solvent, being red shifted to 227 nm, while the dark transition at 242 nm is strongly influenced by the solvent, being shifted to 231 nm and 221 nm with the PCM and PCM + H O schemes, respectively. In the latter model, the larger solvatochromic shift results in an inversion of the ordering between the dark and the bright transitions. For this reason, we adopted the notation of indicating the first weak transition and the fourth bright transition as 1 and 4, respectively, being well separated from the others, while the second and the third are labeled as D (dark) and bright (B), to avoid possible confusion due to the inverse ordering.

Such a progressive red shift found in the wavelengths of the electronic absorption bands reported in Table 4 when moving from the isolated molecules to the solvated systems is in agreement with previous works.44, 52, 53, 111, 112 Marenich et al.111, 112 found that the presence of the directional interaction induces a red shift in the first excited state transition, showing that the average position of the first excited state band goes from 263.1 nm in the isolated pyridine to 254.0 nm with the C‐PCM52, 53 scheme, finally reaching the value of 248.5 when the hydrogen bonded water molecule is explicitly included. A similar result (249.5 nm) was obtained by means of QM/MM calculations,44 where the water molecules surrounding pyridine were explicitly described through the multilayer model based on the FQ approach.

If we compare the transitions reported in Table 4 with the experimental UV‐vis vertical spectrum of Figure 8, we notice that significant discrepancies between the wavelengths of the two bright transitions and the experimental maxima are observed. Such discrepancies, are partially due to the fact that vertical excitation energies (and oscillator strengths) of pyridine in gas phase, predicted by computational approaches are quite dependent on the specific method adopted.103, 104, 105, 106, 107, 108, 109, 110 Hence, the results discussed in the following should be considered as a best estimate of the real spectrum within the B3LYP/6‐31 + G(d) level of theory.

Figure 8

UV‐vis spectrum of pyridine in aqueous solution. The black line is the experimental spectrum113; the blue dashed line is the spectrum obtained by convolution with Gaussian functions (HWHM = 200 cm ) of the vertical spectra taken from 100 snapshots of the MD simulation (1 configuration every 50 ps). The red line is the spectrum obtained using Gaussian functions with HWHM = 500 cm in the convolution. [Color figure can be viewed at wileyonlinelibrary.com]

In particular, the maximum of the experimental first absorption band is observed at about 257 nm, and this value could be either compared with the B absorption band at 227 nm, and/or with the transition 1, predicted at 247 and 239 nm by the PCM and PCM + H O models, respectively. In the first case, the discrepancy is due to a wrong solvatochromic shift predicted in the wavelength of transition B, while in the second case, the error on the wavelength of transition 1 would be smaller but the intensity should be much more enhanced by the solvent. Furthermore, the experimental band shape appears to be very complex, with the presence of many minor absorption peaks at 246, 251, 262 nm, suggesting that beyond the solvation effects modeled by the PCM and PCM + H O schemes, other effects are likely to affect the overall spectral line‐shape.

The inclusion of vibronic features and broadening effects

To improve the agreement between simulated and experimental spectra, we will mainly focus on two main effects in the following: the solvent broadening effects coming from the intrinsic disorder of the molecule in solvent (i), and the effects coming from the vibronic structure of the electronic transitions (ii).

The former effect (i) can be included by the MD. In Table 4 averages and standard deviations of wavelengths and oscillator strengths of the first six states computed with the TD‐DFT/FQ/FX model on the 100 snapshots extracted from the MD simulation of pyridine in water, are labeled as FQ(MD). We observe that transition 1 and 4, falling at 247 and 202 nm, respectively, are well separated from the others, and the intensity of transition 1 appears increased with respect to the PCM and PCM + H O models, suggesting that the intrinsic disorder of the geometry of the snapshots leads to an average increase of the oscillator strength. On the other side, the bright and dark character of the second and third transition cannot be recognized anymore, and they appear as two mixed bright transitions with lower oscillator strengths than the ones predicted by the PCM and PCM + H O models. This could be explained by recalling that while the wavelength of the bright transition is poorly affected by the solvent, the dark transition shows a large solvatochromic shift, possibly leading to the inversion of the levels; thus, there still could exist a Bright and Dark character for the second and third transitions, but as the two states are randomly inverted in the different snapshots of the MD, they appear mixed in the averages. As a first level of approximation we simply neglected such an issue adopting a very simple and intuitive approach of calculating the vertical UV‐vis spectrum on each snapshot by convoluting transition energies and oscillator strengths of each transition with the same Gaussian function with fixed HWHM = 200 cm , and for comparison with HWHM = 500 cm , independently of the assignment. Then, the final the UV‐vis spectrum has been obtained as the sum of the contributions from each snapshot, with normalized intensities. This has been done using the VMS program114, 115, 116 recently developed by our group and some handmade programs written on purpose, and the final spectrum is shown in dashed line in Figure 8, together with the experimental one.

From such an approach we observe that there is quite a broad distribution of peaks in the low energy zone of the spectrum, possibly due to the fact that transition 1 in the MD simulation has higher average oscillator strengths with respect to the PCM and PCM + H O models and a broad distribution of wavelengths, as suggested by the value of standard deviation of 11 nm (Table 4). In this respect, the small peaks at about 240, 250, 260, and 270 nm of the FQ(MD) spectrum can be assigned to such a transition. The strong absorption band between 225 and 235 nm is due to the contributions from the mixed transitions M1 and M2, being worth of note that the solvent broadening leads to a red‐shift of the wavelength with respect to the value of 227 nm found with the PCM and PCM + H O model (Table 4), approaching the experimental maxima. Comparing the simulated and experimental band‐shape, it seems that the simulated spectrum generally reproduces relative intensities and wavelengths differences of the transitions, and just suffers from a rigid blue‐shift of few tens of nanometers. There are many reasons which could possibly lead to such a shift. While some sources of error (such as the particular choice of the DFT functional and basis set, the parametrizations of the FQ model, MM force field and function, the lack of dispersion effects in the description of the interactions between the orbitals of pyridine and the water molecules) cannot be removed being intrinsically related with the fact that the medium/large size of the molecular system under study requires the use of approximate methodologies, other effects such as the vibronic substructure underlying the electronic transitions can be analyzed in more detail.

Indeed, further insights in the reproduction of the band‐shapes can be added by the vibronic effects (ii). We applied the VG approach50, 51 to simulate the vibrational structure of the transitions 1, B and 4 in the PCM + H O model, and in Figure 9 the corresponding spectra are shown.

Figure 9

Vertical gradient spectrum of pyridine in aqueous solution according to the PCM + H O model. The red line is the spectrum obtained using Gaussian functions with HWHM = 200 cm−1 in the convolution. [Color figure can be viewed at wileyonlinelibrary.com]

Transition 1 has an absorption maximum at about 240 nm, transition B has two main maxima at about 230 and 235 nm, and transition 4 has two maxima at about 200 and 205 nm. Comparing such values with the vertical absorption wavelengths of the PCM + H O model reported in Table 4, the general macroscopic effect coming from the vibronic structure seems to be an overall red‐shift of the absorption maxima with respect to the vertical values of about 5–7 nm. As such a shift is very similar for all the three observed transitions, it could be reasonable to assume that the vibrational features contribute to rigidly shift the overall simulated spectrum toward higher wavelength zone, at least in this region. Furthermore, we observe that the vibronic line‐shapes of transitions B and 4 are far from being similar to a Gaussian line‐shape, suggesting that using Gaussian functions with fixed HWHM for the convolution of the vertical transitions of the snapshots is a coarse approximation.

As in the real system the vibronic effects and the solvent broadening should be both present at the same time, we tried to combine these two effect in a simple and empirical way for the first absorption band, by convoluting the vertical excitation energies obtained along the MD simulation with the vibronic line‐shapes of the transitions 1 and B, in place of using simple Gaussian functions.

To do that, we have to be sure that we are convoluting the line‐shape of the right transition with the right excitation wavelength for each snapshot. From one side, the FQ(MD) and FQ(H O) (Ref. 44) schemes predict the weak transition 1 at 241 and 247 nm, respectively, and the bright transition 4 at 197 and 202 nm, respectively, in agreement with the PCM + H O model. Furthermore, the latter appear to be well separated from the 5th transition, occurring at 184 and 183 nm, respectively. On the other side, in many snapshots of the FQ(MD) model the dark/bright assignment of the second and third transition is lost, and they appear mixed.

We decided to proceed in an approximate way, using the ratio between the oscillator strengths of the different transitions as a criterion to determine whether a given snapshot shows an electronic structure similar to the PCM + H O one, or not. In particular, we retained only those snapshots showing a ratio between oscillator strengths of the second and third transitions higher than an arbitrary threshold value (3.0). In this way, we were sure that the contributions from the bright and dark transitions were not mixed, being similar to the PCM + H O one. In Figure 10, the distribution of the second and third transition energies and oscillator strengths is plotted.

Figure 10

Assignment of the second and third excited states from the FQ(MD) simulation, in Bright/Dark states (71 snapshots) and Unknown states (29 snapshots). [Color figure can be viewed at wileyonlinelibrary.com]

By applying the aforementioned criterion, we found that in 71 snapshots the second and third states presented very different oscillator strengths, and could be assigned as bright and dark transitions. Otherwise in 29 snapshots, those state presented similar oscillator strengths, so that nothing could be concluded about the assignment (we labeled those states as ″Mixed” or ″Unknown” states, as no comparison could be made with the reference PCM + H O spectral features). The values of wavelengths and oscillator strengths averaged over the 71 selected snapshots are reported in Table 4 as FQ (MD), together with the standard deviation. The average value of wavelength and oscillator strengths of transition B are 229 nm and 0.0390, respectively, in agreement with PCM and PCM + H O predictions. On the other side, transition D, which is very sensitive to the solvation environment is found at about 230 nm, in close agreement with the PCM model. It is worth noting that the Mixed transitions are found at 231 and 224 nm, very close to the PCM + H O values, and this could suggest that the explicit water directly interacting with the QM nitrogen induces some mixing in the oscillator strengths of the bright and dark transitions. The average wavelength for transition 1 is of 248 nm, with a standard deviation of 11 nm. As the next transition occurs at 230 m, the risk of level inversion should be very small. Furthermore, such an analysis confirms the enhancement of the oscillator strength of transition 1 (0.0112) with respect to the PCM and PCM + H O models, as was also found with the FQ(MD) approach. The average wavelength of transition 4 is 202 nm, with a standard deviation of 3 nm, and the other closest transitions occur at 224 and 183 nm, being well separated.

To be sure that such results are not dependent on some particular computational choice, such as the parametrization of the FQ model for TD‐DFT/FQ/FX calculations or some unphysical arrangement of the solvent molecules in the MD simulations, we removed from the snapshots all the water molecules except the one which was found to be directly bonded with the QM nitrogen (from the function analysis), and computed the vertical absorption spectrum with TD‐DFT/PCM calculations. The results are referred in Table 4 as PCM(MD) and PCM (MD), and the values of wavelengths and oscillator strengths are in very close agreement with the corresponding FQ(MD) and FQ (MD) calculations.

Then, if a given snapshot showed an electronic structure similar to the PCM + H O one, we projected the vibronic line‐shape on the transitions 1 and B, assuming that the spectral structure of such a snapshot should be reasonably similar to the PCM + H O one; in the other cases we simply rejected that snapshot without considering its contribution to the final overall spectrum.

In Figure 11, the first absorption band obtained by convoluting the vibronic features of the 1 and B transitions of the PCM + H O system with the vertical energies obtained from the 71 selected snapshots is shown.

Figure 11

First absorption band of pyridine in aqueous solution. Upper panel: The vibronic features of transition 1 and B of the PCM + H2O model have been projected onto the vertical excitation energies coming from the MD simulation. Both bands normalized to one. Lower panel: The black line is the experimental spectrum, the dashed line is the FQ(MD) spectrum, the straight red and blue lines are the spectra of the 1 and B transitions showed in the upper panel weighted with a ratio of 1:3, according to the oscillator strengths of the FQ model shown in Table 4. [Color figure can be viewed at wileyonlinelibrary.com]

In the upper panel of Figure 11, vibronic stick spectra and convoluted line‐shapes of the two transitions are shown normalized on the same scale, to evidence the comparison. In the lower panel of Figure 11, the convoluted spectra of the two transitions have been plotted together with the experimental and the FQ(MD) line‐shapes. Furthermore, the intensities have been scaled to keep into account the enhancement of the oscillator strength of transition 1 found in the MD simulation, so that the intensity of transition 1 has been set to be 29% of transition B, according to the average values of the oscillator strengths of 0.0112 and 0.0390, respectively (Table 4).

In this way, the red shift coming from the vibronic features and the intensity enhancement coming from the solvation effects can be visualized together, showing a better agreement with the experimental spectrum with respect to the simple Gaussian functions used in the FQ(MD) model, both in terms of wavelengths and line‐shapes. It is worth noting that although the effect of the Mixed states has been neglected in this approach, as they were found to be important for only the 30% of the snapshots, they still should give a contribution to the final spectrum to some extent.

Eventually, from our calculations we can conclude that the red shift of the UV‐vis spectrum of pyridine, moving from the gas phase to the aqueous solution should be mainly due to a combination of effects: a moderate red shift of the spectrum to the vibronic structure, an enhancement of the oscillator strength of the transition due to the solvent broadening, and the mixing of the oscillator strengths of the B and D transitions possibly due to the fact that in the solvent environment the A and B symmetries of the dark and bright transitions are lost.

UV‐vis spectrum of nicotine in water

The experimental electronic absorption spectrum of nicotine at basic pH,101, 102 resembles the one of pyridine shown in Figure 8, in agreement with the detailed description reported by Clayton et al.101, 102 suggesting that the UV spectrum of nicotine in water is dominated by the transition related to the pyridine ring. Thus, the same computational protocol adopted for pyridine and described in the previous paragraph has been used for nicotine.

In Table 5, the vertical excitation energies of the first seven excited states are reported, as computed using the PCM + 2H O, FQ(MD) and FQ (MD) models. In the high wavelength zone of the spectrum of the nicotine A conformer, two weak transitions related to the pyrrolidine ring of nicotine are found, and they have been labeled as 1′ and 2′. Then, moving to lower wavelengths the electronic structure of pyridine is recovered, with the two transitions at 241 and 221 nm (labeled as 1 and D according to Table 4), and the two transitions at 232, 207 nm (B and 4). The transitions of nicotine B are very similar to those of nicotine A, with minor changes in wavelengths and oscillator strengths, with the only noticeable remark that an inversion occur between the transition 4 and the 5th.

Table 5

Wavelength (nm) and oscillator strengths (f) for the first eight excited states of nicotine in aqueous solution.

		PCM+2H _2O Nicotine A		Nicotine B		FQ(MD)			FQ _B/D(MD)
Label	Assign.	nm	f	nm	f	〈nm〉	〈f〉	N	〈nm〉	〈f〉	N
1′	n^′→π^⋆	267	0.0033	264	0.0069	〈287〉_23	〈0.0059〉_0.0061	100	〈288〉_24	〈0.0056〉_0.0045	68
2′	n^′→π^⋆	242	0.0026	249	0.0011	〈264〉_15	〈0.0087〉_0.0100	100	〈264〉_16	〈0.0078〉_0.0087	68
1	n→π^⋆	241	0.0046	237	0.0058	〈246〉_8	〈0.0176〉_0.0178	100	〈246〉_8	〈0.0155〉_0.0148	68
B	π→π^⋆(Bright)	234	0.0942	232	0.0769				〈233〉_6	〈0.0539〉_0.0220	68
D	n→π^⋆(Dark)	221	0.0019	223	0.0004				〈232〉_9	〈0.0071〉_0.0058	68
M1	(Mixed)					〈237〉_6	〈0.0317〉_0.0242	100	〈236〉_6	〈0.0277〉_0.0126	32
M2	(Mixed)					〈228〉_7	〈0.0276〉_0.0254	100	〈227〉_8	〈0.0281〉_0.0167	32
4	π→π^⋆	207	0.0349	205	0.0459	〈217〉_7	〈0.0396〉_0.0401	100	〈217〉_7	〈0.0398〉_0.0398	68
5th	–	204	0.0423	208	0.0613	〈208〉_7	〈0.0343〉_0.0349	100	〈209〉_7	〈0.0329〉_0.0352	68
6th	–	195	0.0000	194	0.1215	〈202〉_7	〈0.0346〉_0.0413	100	〈202〉_7	〈0.0366〉_0.0456	68

PCM + 2H O refer to TD‐DFT calculations on the optimized pyridine in PCM with two explicit hydrogen bonded QM water molecules. FQ(MD) refers to wavelengths and oscillator strengths averaged over the snapshots of the MD simulation, using the TD‐DFT/FQ/FX model. The averages including the Bright/Dark states analysis are reported as FQ (MD) (see text). M1 and M2 refer to the Mixed states found in the MD simulations. Standard deviations are given as subscripts.

Moving to the FQ(MD) averages, the transitions located on the pyrrolidine ring seem to be quite well separated from the others, being the average wavelengths at 287 and 264 nm, and the values of oscillator strengths relatively close to the ones found with the PCM + 2H O model. Then we found the first transition at 246 nm shows an increased intensity with respect to the PCM + 2H O calculations, and the second and third transitions, again with values of oscillator strengths intermediate between the values found for the B and D transitions with the PCM + 2H O model. In Figure 12, the spectra computed by convoluting the values of wavelengths and oscillator strengths in the 100 snapshots with Gaussian functions at fixed HWHM = 200 cm and HWHM = 500 cm are shown. Besides the issues related to the specific choice of the computational methods used for the calculations,48 the discrepancies between the experimental and simulated spectra can be related either to an underestimation of the wavelength of transition B, or to an underestimation of the intensity of transitions 1′, 2′, and 1. Regarding the latter point, we recall that the oscillator strength of transition 1 in the FQ(MD) approach is larger than the corresponding value computed with the PCM + 2H O model, and in Figure 12 the peak related to such a transition appears at about 245–250 nm.

Figure 12

UV‐vis spectrum of nicotine in aqueous solution. The black line is the experimental spectrum; the dashed blue line is the spectrum obtained by the MD simulation configurations (1 snapshot every 1 ns of the MD simulation) convolved with Gaussian functions at fixed HWHM = 200 cm . The red line is the spectrum obtained using Gaussian functions with HWHM = 500 cm in the convolution. [Color figure can be viewed at wileyonlinelibrary.com]

To split the effect of the inversion of the transitions due to similar wavelengths, from the physical effect of mixing of oscillator strengths, we applied the same criterion used for the pyridine and presented in the previous section, and the results are presented in the FQ (MD) column. Also in this case about 70% of the snapshots showed an electronic structure similar to the PCM + 2H O model, with one bright transition at 233 nm, and the dark one at 232 nm, while in the remaining 30% of the snapshots the mixing of the oscillator strengths occurred due to the breaking of the symmetry of the transitions. The other transition seems to be quite unaffected by the level inversion as the values of averages and standard deviations for both the wavelength and the oscillator strengths are quite similar between the FQ(MD) and the FQ (MD) models. Furthermore, the oscillator strength of transition 1 in this model is about 28% of the intensity of transition B, in quantitative agreement with the results for pyridine.

Moving to the vibronic effects, we found a major difference with respect to the pyridine spectrum. In Figure 13, the VG spectra of the B transition for both the nicotine A and B conformers are shown, and the vibronic line‐shape resembles the one of a Gaussian function, suggesting that the vibronic shift should be smaller in this case than for the pyridine. To be sure that such a Gaussian shape is really due to the vibronic structure, we used a very small HWHM of just 20 cm , evidencing the subtle features of the spectra.

Figure 13

Vertical gradient spectrum of the transitions 1 and B of nicotine A (a) and nicotine B (b) with the PCM + 2H O model. [Color figure can be viewed at wileyonlinelibrary.com]

Finally, in Figure 14, we plotted the convolution of the vibronic features related to transitions 1 and B, with the wavelengths and oscillator strengths of the 68 snapshots selected in the FQ approach. While in the upper panel the two transition have been plotted together and both scaled to unity, to evidence the two separate contributions, in the lower panel the convoluted line‐shapes of transition 1 has been set to be the 28% of the one of transition B, to show the correct contributions, in agreement with the oscillator strengths of Table 5. In this case, we note that the vibronic features do not shift the maximum of the absorption band, but provide only broadening effects.

Figure 14

Upper panel: Convolution of the vibronic features related to transitions 1 and B of nicotine in PCM + H2O, with the vertical excitation energies coming from the shapshots of the MD simulation. Lower panel: The black line represent the experimental spectrum of nicotine in aqueous solution. The red line is the spectrum obtained by convolution of the vibronic structure of the state B, whereas the cyan line is the spectrum obtained by convolution of the vibronic structure of the state 1. The blue dashed line is the spectrum obtained without vibronic resolution. [Color figure can be viewed at wileyonlinelibrary.com]

Conclusions

The electronic absorption spectra of molecules in aqueous solution have been determined with an integrated computational approach made up by MD simulations of the target molecule in water followed by TD‐DFT calculation in a QM/MM framework with the solvent molecules described through fluctuating or FX models depending on the distance from the solute.44 Both the intramolecular potential and atomic charges of the solute molecule in MD simulations have been accurately parametrized to improve the description of the intermolecular interactions, with particular regard to the hydrogen bonds. The electronic states have been computed on a series of configurations extracted by MD trajectories. Because it has been observed44 that the MM water molecules directly interacting with the QM part require parameters which differ from those adopted in the bulk phase, these are locate using the function developed by Pagliai et al.71 To avoid polarization problems of the most external fluctuating charges water molecules, these have been replaced by the FX characterizing the TIP3P‐FB model. This approach has allowed an accurate description of the structural variation of nicotine molecule in aqueous solution, observing the presence of two different conformers in agreement with experimental data.99, 100 As deduced by different measurements,49, 101, 102 the N5 atom forms more stable hydrogen bonds than N6, even if a VS has been needed to describe the intermolecular directional character due to the ibridization of the latter. The correct interpretation of the solute/solvent interactions has been a valuable help in the definition of a model for the electronic absorption spectra, taking into account the vibrational effects through vibronic calculations. Together with the specific interest of the studied molecules, our results pave the route toward systematic use of the proposed tool as a companion to experimental spectroscopic studies.

Supporting Information

Click here for additional data file.