Assessing Molecular Dynamics Simulations with Solvatochromism

Article pubs.acs.org/JPCB

Assessing Molecular Dynamics Simulations with Solvatochromism Modeling Tobias Schwabe* Center for Bioinformatics and Physical Chemistry Institute, University of Hamburg, Bundesstraße 43, D-20146 Hamburg, Germany

Downloaded by SUNY UPSTATE MEDICAL UNIV on September 9, 2015 | http://pubs.acs.org Publication Date (Web): August 6, 2015 | doi: 10.1021/acs.jpcb.5b05206

S Supporting Information *

ABSTRACT: For the modeling of solvatochromism with an explicit representation of the solvent molecules, the quality of preceding molecular dynamics simulations is crucial. Therefore, the possibility to apply force fields which are derived with as little empiricism as possible seems desirable. Such an approach is tested here by exploiting the sensitive solvatochromism of p-nitroaniline, and the use of reliable excitation energies based on approximate second-order coupled cluster results within a polarizable embedding scheme. The quality of the various MD settings for four different solvents, water, methanol, ethanol, and dichloromethane, is assessed. In general, good agreement with the experiment is observed when polarizable force fields and special treatment of hydrogen bonding are applied.

■

methods.19 These embedding approaches can range from very elaborate and computational demanding subsystem DFT20−22 and wave function-in-frozen density approaches23−25 to simple point charge embedding.26,27 The latter is also referred to as quantum mechanics/molecular mechanics approach (QM/MM) and is very efficient with almost no additional computational cost in comparison to an equivalent vacuum computation but lack the electronic response of the environment (solvent) to an electronic excitation in the solute and the mutual coupling. Because this static picture does not always hold, the coupling can be introduced via the computation of induced dipoles or fluctuating charges in the classical (solvent) region. Several such approaches have been combined with quantum mechanical methods for the study of the surroundings on response properties28−42 and have been used before with semiempirical methods.43,44 Actually, already the very first paper on QM/MM made use of induced dipoles.45 Note that such advanced schemes can also be combined with continuum solvation models to reduce the size of the region with explicit solvent molecules.46−48 One of the most elaborate schemes is the polarizable embedding (PE) approach which has been implemented for various quantum mechanical methods.34,36,40−42 It takes the coupling to the environment also for the electronic response into account which has been shown to be crucial for an accurate description of environmental effects on electronic excitations. 41,49 It has been applied successfully to model solvatochromism34,36,50,51 with a special focus on two-photon

INTRODUCTION Solvatochromism has been of constant interest for spectroscopic studies because on the one hand it allows for the analysis of the interactions between solute and solvent and on the other hand reveals the character of an excited state. Therefore, it has also been of great interest to theoretical chemistry because of a fruitful interplay: theory aids in interpreting experiments, and experiments on solvatochromism allow one to assess theoretical models of solutions because of the relative ease of comparing computed and measured results directly.1,2 For modeling solvatochromism, the complete solute/solvent system (or at least a large part of the solute’s solvation shell) should be treated quantum mechanically. Due to computational resource constraints, only for a few studies, this approach was chosen then it was often necessary to resort to semiempirical quantum mechanical methods.3−6 In most cases, two less expensive schemes are applied: only the solute is described on an adequate quantum mechanical level and electronic excitations are locally confined to this region while solvent molecules are approximated by either an implicit or an explicit but lower-level representation. Both schemes are still subject to current developments and improvements. For the former, the molecules are replaced by an interaction with a dielectric continuum. Early examples which are still in use are PCM7 and COSMO8 and a more recent development in the field is SMSSP.9 Especially PCM and its further improvements in combination with time-dependent density functional theory (DFT) are very common for the study of solvent effects on electronic excitations,10−13 but also other quantum mechanical methods in combination with continuum solvation models have regained attention lately.14−18 Schemes with an explicit representation of solvent molecules are sometimes subsumed under the label embedding © 2015 American Chemical Society

Received: June 1, 2015 Revised: July 14, 2015 Published: July 28, 2015 10693

DOI: 10.1021/acs.jpcb.5b05206 J. Phys. Chem. B 2015, 119, 10693−10700


The Journal of Physical Chemistry B

■

absorption52−55 but also for other properties in solution56 or to describe effects of the protein surroundings.57−62 Beside the very reliable results for the actual PE computation, another appealing feature is the application of classical potentials which are derived from quantum chemical computations and are applied without further parametrization. These potentials can also be used for the molecular dynamics (MD) simulation which is required to obtain a set of independent solute/solvent configurations which resemble the real distribution in liquid phase. Computing these configurations is a key step for all schemes which make use of an explicit representation of solvent molecules. In principle, this can be done with ab initio methods like Car−Parrinello MD63 or semiempirical approaches likes density functional tight binding,64 but in most cases, the computational costs are prohibitive and therefore QM/MM schemes or pure molecular mechanics simulations are used. In the latter cases, the quality of the configurations depends on the quality of the applied force field to large extent and any deficiencies here can ruin the subsequent computations to study the solvatochromic effect. Modern force fields are highly parametrized but normally with other focuses than chromophore contacts in solution. Therefore, there are no guarantees that such solvent force fields lead to reliable results for the study of solvatochromism. Resorting to interacting potentials with low empiricism might relieve the burden of choosing the right one, and indeed, good results with the PE potentials could be obtained in previous studies.34,50 Especially, the PE potentials allow for a straight derivation of atomic polarizabilities so that the more advanced interaction scheme of polarizable force fields can be used. Polarizable force fields, at least for highly polarizable solvents like water, are deemed to be necessary for reliable MD simulations.65,66 For subsequent studies in our lab with the same setup, it turned out that the treatment of hydrogen bonding in protic solvents is a crucial factor. Especially, it was found that removing the polarizability from H atoms was beneficial in some cases. But such a treatment cannot be justified based on the computations of local polarizabilities in vacuum. More empiric force fields treat the polarizability at the H donor differently but normally without any deeper rationalization.65−70 Normally, main arguments for removing it are related to the efficiency of the approach. Another example is the quantum mechanically derived force field (QMDFF),71 which includes an explicit hydrogen bond correction72 without further justification. A systematic investigation of the H-bond treatment in the context of solvatochromism modeling is an open question. In this paper, this question is addressed by computing solvent shifts of the electronic excitation energy for the first π−π*excitation in p-nitroaniline (PNA) in various solvents and based on different MD settings. For rating the observed differences, other parameters than the H-bond treatment are also varied. PNA is known to be quite sensitive to solvent effects, which are exploited to study Kamlet−Taft solvent parameters1 or preferential solution in solvent mixture73 and which also still attracts interest from the theoretical side for development testing.37,49,74−79 Here, this sensitivity of PNA and reliable excitation energies including solvent effects based on approximate second-order coupled cluster results with polarizable embedding (PERI-CC2)40 are used to assess the quality of the MD settings in the preceding simulations.

Article

COMPUTATIONAL DETAILS

Potential Determination. The multipole moments and polarizabilities for water and for methanol are taken from a previous study.50 Missing data for PNA in ethanol (representative for all protic solvents) and in dichloromethane, for dichloromethane (DCM) itself and ethanol, are computed by the same protocol as before. First, the molecules are optimized at the B3LYP80,81/aug-cc-pVTZ82,83/PCM84,85 level of theory where we apply the standard PCM parameters as provided in GAUSSIAN 09.86 Next, multipole moments up to quadrupoles and anisotropic polarizabilities are obtained with the LoProp87 approach as implemented in MOLCAS.88 The properties are calculated at the B3LYP/a-aug-cc-pVTZ82,83 level for the isolated molecules in vacuum using the geometry obtained from the PCM optimization. The a-aug-cc-pVTZ basis set is a recontraction of aug-cc-pVTZ to an ANO-type basis, which is necessary for the LoProp procedure. Finally, we also obtain atomic point charges with the CHelpG procedure89 fitted against the electrostatic potential from a B3LYP/aug-cc-pVTZ calculation again at the PCM optimized geometries. As an additional constraint, the reproduction of the molecular dipole moment is imposed on the fitting. These point charges implicitly contain higher-order multipole effects and therefore should not be combined with higher multipole moments. For comparison, solvent potentials have been recomputed at the HF/ANO-S level of theory for the LoProp approach and HF/cc-pVDZ for the CHelpG parameter and are used to replace either the polarizabilities or the point charges or both. All potential data which have been obtained for this study are given in the Supporting Information. Molecular Dynamics Simulations. Periodic boundary conditions are applied for all molecular simulations and are based on a cubic box of 31.925 Å edge length. In each box, one PNA molecule was placed, and the box was filled with remaining solvent molecules to obtain the experimental density of the solvent. This corresponds to a PNA concentration of 0.05 mol l−1. An equilibration run of 0.4 ns in time steps of 2 fs followed by a production run of 1.2 ns is carried out. Configurations are sampled every 10 ps, which should ensure statistically uncorrelated configurations and results in 120 configurations from each simulation. All molecules are treated as rigid bodies and thus only intermolecular potentials are considered. Van-der-Waals interactions are based on a 12-6-Lennard−Jones potential. The Lennard−Jones parameters are taken from the OPLS-AA force field90 as shipped with TINKER.91 Further, atomic point charges and induced dipoles from atomic isotropic polarizabilities are taken into account. On the basis of these common settings, different simulation setups have been chosen. As reference, an NVT ensemble at T = 298 K is used. All interactions up to a cutoff radius rcut = 15 Å are considered beyond which a reaction field method based on experimental dielectric permittivities is applied to avoid boundary effects for the induced dipoles. The dielectric permittivities are the same as in the PCM computations and are taken from GAUSSIAN. Atomic point charges and polarizabilities are taken from the B3LYP computations with the larger triple-ζ basis set (see previous section). For H atoms at H-bond donors, Lennard−Jones parameters and polarizabilities are set to zero. To conserve the molecular 10694


Article


Table 1. Comparison of Electronic Excitation Energies and Solvent Shifts of PNA in Different Solvents and Based on Various MD Setupsa aq


a

MeOH

EtOH

DCM

MD setup

ΔE(exc)

shift

ΔE(exc)

shift

ΔE(exc)

shift

ΔE(exc)

shift

exp NVT NPT rcut =10 Å pol(HF) all(HF) M*P0 pot M*P0(HF) pot pol H-bond no H-bond

3.29 3.30 3.32 3.31 3.31 3.25 3.40 3.42 3.32 3.57

−1.00 −0.99 −0.97 −0.98 −0.98 −1.04 −0.89 −0.87 −0.97 −0.72

3.37 3.46 3.51 3.49 3.50 3.54 3.58 3.60 3.56 3.71

−0.92 −0.83 −0.78 −0.80 −0.79 −0.75 −0.71 −0.69 −0.73 −0.58

3.36 3.51 3.42 3.45 3.48 3.52 3.57 3.60 3.59 3.72

−0.93 −0.78 −0.87 −0.84 −0.81 −0.77 −0.72 −0.69 −0.70 −0.57

3.58 3.63 3.69 3.65 3.64 3.65 3.63 3.61 3.63 3.64

−0.71 −0.66 −0.60 −0.64 −0.65 −0.64 −0.66 −0.68 −0.66 −0.65

All computed data are obtained by averaging over 120 configurations. Experimental reference data are taken from Ref 73. All values in eV.

■

RESULTS AND DISCUSSION The electronic excitation energies and solvatochromic shifts of the bright π−π*-state of PNA in aqueous, methanol, ethanol, and dichloromethane solution based on the various MD simulation settings are compiled in Table 1. The transition of interest is not the lowest excitation in gas phase but is subject to a larger solvent shift than the lower excitations. Therefore, always the lowest excitation energy found for a given configuration was chosen. As a check, oscillator strength was also inspected and found to be consistent with the π−π*excitation. The NVT setting has usually been applied in previous studies wherefore it is used as a reference setting here. Comparing these results to experimental data, very good agreement for the shifts but also for absolute excitation energies can be found for water and dichloromethane. For methanol, the agreement is still satisfying and the results deviate less than 0.1 eV from the experimental ones. For ethanol, the results are a bit worse and deviate by 0.15 eV, which is slightly higher than what has been found before for PERI-CC2 computations. Next, the influence of the varying simulation settings is considered. For dichloromethane, results are almost insensitive to the underlying simulation and are not changing by more than 0.08 eV. Even computing the force field parameter at the lower HF level does not influcence the results significantly. For most of the aqueous solution simulations, results are also quite stable. But when induced dipoles are not considered, the solvatochromic shift decreases by 0.10 eV (or even 0.12 eV when HF point charges are used). For the configurations obtained from these simulations, the excited state is less stabilized than what is found in the experiment. This gets worse when Lennard−Jones parameters are placed on the water H atoms. Now, the solvatochromic shift is only −0.72 eV instead of −0.99 eV. It should be emphasized that this is only a consequence of the underlying configurations used for the PERI-CC2 computations. The latter are always carried out with the same computational setup. Of course, the chosen parameters have not been trained for being used in H bonds. Instead, they have been taken from H atom data for different chemical environments. Further, the parameter for the oxygen atom in water has also not been reoptimized. This inexpedient choice has been made deliberately to study how sensitive PE results are to such more drastic changes in the underlying simulations.

polarizability, the polarizability from the H atom was equally distributed over all other atoms. This reference setting has been altered in the following ways: (i) the NPT ensemble is applied, using experimental compressibilities (χH2O = 0.51 GPa−1,92 χMeOH = 1.2 GPa−1,93 χEtOH = 1.1 GPa−1,93 χDCM = 1.0 GPa−194), (ii) rcut is reduced to 10 Å, (iii) polarizabilities are taken from a LoProp computation with Hartree−Fock and a double-ζ basis set, (iv) atomic point charges and polarizabilities are computed with Hartree−Fock and small basis sets of double-ζ quality, (v) only point charges but no induced dipoles are used (dubbed M*P0 here for consistency with previous studies), (vi) only point charges from the HF computations but no induced dipoles are used, (vii) no special care of polarizabilities at H atoms on H-bond donors, and (viii) no special treatment for H atoms on H-bond donors at all. All simulations have been carried out with the MOLSIM program package.95 Polarizable Embedding Computations for Electronic Excitations. The quantum mechanical treatment following the molecular simulation step has been the same regardless of the simulation setup. The input for each PE computation is generated with the WHIRLPOOL program.96 For each configuration, the solute geometry is extracted, and a cutoff radius of 12 Å from the center of mass of the solute is applied. All solvent molecules within the cutoff are included in the electrostatic embedding potential which is built from multipole moments up to quadrupoles and anisotropic polarizabilities. Because it has been shown recently that density functional approximations are not capable of treating solvatochromism in PNA right,77 PE calculations are performed using PERI-CC240/ (aug-)cc-pVDZ,82,83 that is, diffuse functions are only added on nonhydrogen atoms. For the resolution-of-the-identity (RI) approximations, the corresponding basis sets have been used.97,98 Hartree−Fock reference states are computed with a direct SCF approach,99 which takes the polarizable embedding already at this level into account. The final excitation energy in solution is obtained by averaging over all 120 individual results. To determine solvatochromic shifts, excitation energies for PNA in vacuum based on a B3LYP/aug-cc-pVTZ geometry are computed also at the RI-CC2100,101/(aug-)cc-pVDZ level. All quantum chemical computations in this step are done with a development version of the TURBOMOLE program.102 10695


Article



due to the solute interacting with the polarizable embedding are missing. Nevertheless, the results should reveal general differences in the electric field felt by the solute. Perusing the obtained data, most variation of the electric field strength from different simulation settings is found at the nitro end of PNA (site 3), except for simulations with dicholormethane. This agrees well with the finding that the solvatochromic shift from the latter simulations is quite stable. Comparing to Table 1, a large field strength at site 3 can be correlated to a large solvatochromic shift. Conversely, a lower field at this site leads to less agreement with the experiment. The π−π*excitation has a considerable charge transfer character, which results in additional electron density at the nitro group terminus. A larger field strength at this site could support this process, leading to a reduction in excitation energy. The postulated correlation is also in agreement with the findings for ethanol: the largest electric field strength at site 3 is found for the simulation with the NPT ensemble, and it is this simulation for which the smallest deviation to experiment concerning the solvatochromic shift is obtained. As an additional test, radial distribution functions (RDFs) at the two ends of PNA as well as for the solvent−solvent distributions are plotted in Figure 2. For clarity, only the RDFs from the reference, the point charge only, and the two different H-bond settings are given. Consistently with the previous results, for dichloromethane, all RDFs look the same. For the protic solvents, this is different. The distributions from simulations without special treatment of H bonds show distinct different shapes. As could be assumed, the first solvation shell around both ends of PNA form at too large distances. This effect is more pronounced for the alcohols. But not only the distribution around the solute is affected. Likewise, the whole solvent structure changes. Interestingly, this leads to more structure in aqueous solutions, but for the alcohols, the first and second solvation shell can no longer be distinguished. Results from the point charge only force fields also show some differences in the RDFs, but less pronounced. The plots show slight deviations for the probability to find a solvent molecule next to the functional groups in PNA. This could explain the differences in the electric field strength. On the other hand, the RDFs show also that the overall solvent structure could be related to differences in the description of the solvatochromic shift. In the end, it might as well be a coupled effect. Strictly, from the given data one cannot deduce if the change for the interaction with the solute leads to a change in the solvent structure or the other way around, but it

For methanol, the effect of different simulation settings parallels the one for water except that only adding induced dipoles at hydrogen atoms in H bonds deteriorates the result much more, worsening the solvatochromic shift by 0.10 eV. Again, not taking special care of H bonds is even worse, yielding a solvatochromic shift of only −0.58 eV instead of −0.83 eV from the reference setting. Finally, ethanol is the only case where significant improvement compared to experiment can be observed by altering the simulation settings. The best result is obtained for the NPT ensemble, improving the solvatochromic shift by 0.09 eV. In addition, an improvement by 0.06 eV can be found when electrostatic interactions are only fully treated within a cutoff radius of 10 Å. Otherwise, the findings are similar to methanol and water: not considering induced dipoles in the simulations leads to worse agreement, and applying no special treatment for H bonds is even less advantageous. For further understanding of what exactly is responsible for the varying predictions of the solvatochromic shift, the local field within the solvent cavity is tested. For this, the solute atoms are removed and three test sites are defined for which the electric field including all many-body contributions from the induced dipoles in the surroundings is computed. The test sites are placed in the center of the phenyl ring, at the nitrogen atom of the amine group, and at the nitrogen of the nitro group (see also Figure 1). The electric field strength is determined for

Figure 1. Definition of test sites for the electric field probing.

each configuration individually, and the results are averaged for comparison. These results are given in Table 2. Note that in these computations, the additional many-body contributions

Table 2. Average Electric Field Strength and Its Standard Deviation Found at Three Test Sites within PNA (see also Fig 1) in Different Solvents and Based on Various Simulation Settings, Given in 10−3 a.u.

MD setup NVT NPT rcut =10 Å pol(HF) all(HF) M*P0 pot M*P0(HF) pot pol H-bond no H-bond

1 7.9 7.8 8.0 8.5 8.4 6.9 6.1 8.3 7.3

± ± ± ± ± ± ± ± ±

aq

EtOH

DCM

site

site

site

2 3.7 3.5 3.8 3.8 4.0 3.0 2.5 4.1 3.1

9.2 8.1 9.4 9.6 9.1 9.5 8.4 9.3 9.0

± ± ± ± ± ± ± ± ±

3 4.1 3.5 4.4 3.8 4.5 4.1 3.8 4.1 3.9

16.1 15.3 15.6 15.1 16.9 12.6 12.0 14.9 11.4

± ± ± ± ± ± ± ± ±

1 4.6 4.7 4.5 4.8 5.6 3.9 3.2 4.9 4.1

5.5 6.3 5.9 5.6 5.9 4.9 4.7 5.0 4.1

± ± ± ± ± ± ± ± ±

2 2.0 1.6 2.6 1.8 1.9 1.5 1.5 1.6 1.7 10696

8.8 9.3 8.7 8.2 8.8 7.4 8.1 10.3 5.5

± ± ± ± ± ± ± ± ±

3 3.1 3.0 3.0 3.0 3.0 2.6 2.6 3.0 2.3

9.8 10.7 10.0 9.5 9.4 7.5 7.1 6.7 6.7

± ± ± ± ± ± ± ± ±

1 3.5 4.1 3.7 4.0 4.6 3.0 2.9 4.1 2.6

6.3 5.6 6.1 6.3 6.1 5.9 6.1 6.0 6.7

± ± ± ± ± ± ± ± ±

2 2.5 2.2 2.5 2.5 2.2 2.3 2.4 2.6 2.9

6.2 5.6 7.3 6.6 6.4 6.7 6.9 7.4 6.7

± ± ± ± ± ± ± ± ±

3 2.8 2.4 2.9 3.4 2.9 3.1 3.0 2.8 2.7

9.2 8.8 8.8 8.7 9.0 8.4 9.6 9.0 9.5

± ± ± ± ± ± ± ± ±

2.9 3.1 3.1 3.1 3.4 3.3 3.3 3.1 3.2


Article



fields) perform significantly worse for the protic solvents used in the present study. Whether this is common for all protic solvents or which other properties of the solvent lead to this effect cannot be concluded for now and has to be investigated in a following study. The question of how to treat polarizablities at H atoms attached to an H-bond donor is still not fully answered. From this study, it seems preferable to remove induced dipoles here. The exact reason for this needs further investigation. It is also interesting to note that removing polarizabilities at these H atoms for the subsequent PE computations is not beneficial (data not shown). However, that hydrogen bonds have such a particular influence, it is not surprising. Even today, the physical interpretation of this special bond situation is still a matter of debate.103,104 And although the models might differ to varying degrees, they all agree that the nature of the hydrogen bond can only be revealed when it has already formed (i.e., it cannot be deduced from the isolated hydrogen bond donor and the acceptor). This might change, once the general model is derived but is certainly true for the interaction model based on (classical) electrostatics and polarization. Therefore, even quantum mechanically derived force fields need testing and fine-tuning. On the other hand, it is quite encouraging that results obtained from force field parameters at the HF level do not show significantly different results for the study of the solvatochromic shift. It has been observed that the modeling of ethanol solutions gain accuracy when the NPT ensemble or when a smaller cutoff radius is applied. The finding for the NPT ensemble is not transferable to the other solvent molecules. And the small cutoff radius is below what is normally recommended. Therefore, the better results are most likely just a case of error compensation. Where these errors are coming from in detail cannot be answered here. Possible sources are the use of rigid bodies and the therewith involved limitation of degrees of freedom as well as the neglect of hyperpolarization, which should damp the induced dipole interactions. In the end, the deviation from experiment is a little bit larger for the alcohols than for water and dichloromethane. Further testing and analysis of all steps in the sequential MD plus QM/ MM approach is needed to fully understand this. Certainly, a critical step would be going beyond the rigid-body approximation. Now, approaches like the QMDFF allow one to address these questions without the need for tedious force field tuning.71 Work in this direction is being persued at the moment. Anyhow, the results, which are computed with many parameters obtained from ab initio data, are already in good agreement with experiment. The more we learn about our solvent models and the complex interactions in solution, the easier it will be to derive force field parameters with less empiricism. For this purpose, solvatochromism modeling is an interesting alternative to the simulation of macroscopic properties of solvent and solutions because it is more sensitive to contacts and distributions at an atomic level. With the ability to reliably describe environmental effects on electronic excitations with advanced models like the PE method, further progress in that direction is possible.

Figure 2. RDFs from MDs of PNA in (A) water, (B) methanol, (C) ethanol, and (D) dichloromethane for the distance of (i) O(PNA)X(solvent), (ii) NH(PNA)-X(solvent), and (iii) solvent−solvent, where X is the oxygen atom for the ROH solvents and the carbon atom for DCM.

is very likely that different RDFs will also be found when pure solvent simulations based on the different setting are carried out. Again, some increase of this effect might be present due to coupling to a perturbing effect of the solute.

■

CONCLUSION The effect of varying settings for preceding MD simulations on subsequent solvatochromism computing of PNA has been studied carefully with a focus on interaction potentials derived from ab inito computations without further parametrization. In general, very good agreement with the experiment can be obtained, especially for water and dicholormethane as the solvent. For these solvents, the PERI-CC2 results for PNA deviate less than 0.1 eV to the experimental solvatochromic shift as well as to the absolute excitation energy. But it is also found that the performance depends significantly on a reliable MD setup. Obviously, to obtain solvatochromic shifts in agreement with experimental measurements, the interactions of the solvent with the solute have to be modeled correctly. Severe shortcomings like an unbalanced treatment of H bonds therefore are reflected in a less accurate description of solvatochromism. Note also that even models based on point charges only for the electrostatic interactions (which is the standard in most force 10697


Article


■

(17) Caricato, M. A corrected-linear response formalism for the calculation of electronic excitation energies of solvated molecules with the CCSD-PCM method. Comput. Theor. Chem. 2014, 1040−1041, 99−105. (18) Eriksen, J. J.; Solanko, L. M.; NÅbo, L. J.; Wüstner, D.; Sauer, S. P.; Kongsted, J. The second-order polarization propagator approximation (SOPPA) method coupled to the polarizable continuum model. Comput. Theor. Chem. 2014, 1040−1041, 54−60. (19) Severo Pereira Gomes, A.; Jacob, C. R. Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2012, 108, 222−277. (20) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesolowski, T. A. The merits of the frozen-density embedding scheme to model solvatochromic shifts. J. Chem. Phys. 2005, 122, 094115. (21) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. WIREs Comput. Mol. Sci. 2014, 4, 325−362. (22) Wesolowski, T. A. Embedding potentials for excited states of embedded species. J. Chem. Phys. 2014, 140, 18A530. (23) Gomes, A. S. P.; Jacob, C. R.; Visscher, L. Calculation of local excitations in large systems by embedding wave-function theory in density-functional theory. Phys. Chem. Chem. Phys. 2008, 10, 5353− 5362. (24) Höfener, S.; Severo Pereira Gomes, A.; Visscher, L. Molecular properties via a subsystem density functional theory formulation: A common framework for electronic embedding. J. Chem. Phys. 2012, 136, 044104. (25) Höfener, S.; Visscher, L. Calculation of electronic excitations using wave-function in wave-function frozen-density embedding. J. Chem. Phys. 2012, 137, 204120. (26) Parac, M.; Doerr, M.; Marian, C. M.; Thiel, W. QM/MM calculation of solvent effects on absorption spectra of guanine. J. Comput. Chem. 2010, 31, 90−106. (27) Mennucci, B. Modeling environment effects on spectroscopies through QM/classical models. Phys. Chem. Chem. Phys. 2013, 15, 6583−6594. (28) Kongsted, J.; Osted, A.; Christiansen, O.; Mikkelsen, K. V. The QM/MM approach for wavefunctions, energies and response functions within self-consistent field and coupled cluster theories. Mol. Phys. 2002, 100, 1813−1828. (29) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Linear response functions for coupled cluster/molecular mechanics including polarization interactions. J. Chem. Phys. 2003, 118, 1620− 1633. (30) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. A discrete solvent reaction field model within density functional theory. J. Chem. Phys. 2003, 118, 514−521. (31) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. A discrete solvent reaction field model for calculating molecular linear response properties in solution. J. Chem. Phys. 2003, 119, 3800−3809. (32) Nielsen, C. B.; Christiansen, O.; Mikkelsen, K. V.; Kongsted, J. Density functional self-consistent quantum mechanics/molecular mechanics theory for linear and nonlinear molecular properties: Applications to solvated water and formaldehyde. J. Chem. Phys. 2007, 126, 154112. (33) Söderhjelm, p.; Husberg, C.; Strambi, A.; Olivucci, M.; Ryde, U. Protein influence on electronic spectra modeled by multipoles and polarizabilities. J. Chem. Theory Comput. 2009, 5, 649−658. (34) Olsen, J. M.; Aidas, K.; Kongsted, J. Excited states in solution through polarizable embedding. J. Chem. Theory Comput. 2010, 6, 3721−3734. (35) Slipchenko, L. V. Solvation of the excited states of chromophores in polarizable environment: Orbital relaxation versus polarization. J. Phys. Chem. A 2010, 114, 8824−8830. (36) Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. The polarizable embedding coupled cluster method. J. Chem. Phys. 2011, 134, 104108.

ASSOCIATED CONTENT

S Supporting Information *

This material is available free of charge via the Internet at http://pubs.acs.org/. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b05206. A compilation of the newly computed force field potential data is provided (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes Downloaded by SUNY UPSTATE MEDICAL UNIV on September 9, 2015 | http://pubs.acs.org Publication Date (Web): August 6, 2015 | doi: 10.1021/acs.jpcb.5b05206

The authors declare no competing financial interest.

■

REFERENCES

(1) Reichardt, C., Welton, T. Solvents and Solvent Effects in Organic Chemistry; Wiley-VCH: Weinheim, Germany, 2011. (2) Solvation Effects on Molecules and Biomolecules; Canuto, S., Ed.; Springer: Netherlands, 2008. (3) Coutinho, K.; Canuto, S. Solvent effects in emission spectroscopy: A Monte Carlo quantum mechanics study of the n→π* shift of formaldehyde in water. J. Chem. Phys. 2000, 113, 9132−9139. (4) Murugan, N. A.; Rinkevicius, Z.; Ågren, H. Modeling solvatochromism of Nile red in water. Int. J. Quantum Chem. 2011, 111, 1521−1530. (5) Kongsted, J.; Mennucci, B.; Coutinho, K.; Canuto, S. Solvent effects on the electronic absorption spectrum of camphor using continuum, discrete or explicit approaches. Chem. Phys. Lett. 2010, 484, 185−191. (6) Eilmes, A. Solvatochromic probe in molecular solvents: implicit versus explicit solvent model. Theor. Chem. Acc. 2014, 133, 1538− 1550. (7) Cossi, M.; Barone, V. Time-dependent density functional theory for molecules in liquid solutions. J. Chem. Phys. 2001, 115, 4708−4717. (8) Klamt, A. Calculation of UV/Vis spectra in solution. J. Phys. Chem. 1996, 100, 3349−3353. (9) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Uniform treatment of solute-solvent dispersion in the ground and excited electronic states of the solute based on a solvation model with state-specific polarizability. J. Chem. Theory Comput. 2013, 9, 3649−3659. (10) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3094. (11) Mennucci, B.; Cappelli, C.; Guido, C. A.; Cammi, R.; Tomasi, J. Structures and properties of electronically excited chromophores in solution from the polarizable continuum model coupled to the timedependent density functional theory. J. Phys. Chem. A 2009, 113, 3009−3022. (12) Pedone, A.; Biczysko, M.; Barone, V. Environmental effects in computational spectroscopy: Accuracy and interpretation. ChemPhysChem 2010, 11, 1812−1832. (13) Mennucci, B. Polarizable continuum model. WIREs Comput. Mol. Sci. 2012, 2, 386−404. (14) Cammi, R.; Fukuda, R.; Ehara, M.; Nakatsuji, H. Symmetryadapted cluster and symmetry-adapted cluster-configuration interaction method in the polarizable continuum model: Theory of the solvent effect on the electronic excitation of molecules in solution. J. Chem. Phys. 2010, 133, 024104. (15) Cammi, R. Coupled-cluster theory for the polarizable continuum model. III. A response theory for molecules in solution. Int. J. Quantum Chem. 2012, 112, 2547−2560. (16) Lunkenheimer, B.; Köhn, A. Solvent effects on electronically excited states using the conductor-like screening model and the second-order correlated method ADC(2). J. Chem. Theory Comput. 2013, 9, 977−994. 10698


Article


The Journal of Physical Chemistry B (37) Kosenkov, D.; Slipchenko, L. V. Solvent effects on the electronic transitions of p-nitroaniline: A QM/EFP study. J. Phys. Chem. A 2011, 115, 392−401. (38) DeFusco, A.; Ivanic, J.; Schmidt, M. W.; Gordon, M. S. Solventinduced shifts in electronic spectra of uracil. J. Phys. Chem. A 2011, 115, 4574−4582. (39) DeFusco, A.; Minezawa, N.; Slipchenko, L. V.; Zahariev, F.; Gordon, M. S. Modeling solvent effects on electronic excited states. J. Phys. Chem. Lett. 2011, 2, 2184−2192. (40) Schwabe, T.; Sneskov, K.; Olsen, J. M. H.; Kongsted, J.; Christiansen, O.; Hättig, C. PERI-CC2: A polarizable embedded RICC2 method. J. Chem. Theory Comput. 2012, 8, 3274−3283. (41) Eriksen, J. J.; Sauer, S. P. A.; Mikkelsen, K. V.; Jensen, H. J. A.; Kongsted, J. On the importance of excited state dynamic response electron correlation in polarizable embedding methods. J. Comput. Chem. 2012, 33, 2012−2022. (42) Hedegård, E. D.; List, N. H.; Jensen, H. J. A.; Kongsted, J. The multi-configuration self-consistent field method within a polarizable embedded framework. J. Chem. Phys. 2013, 139, 044101. (43) Thompson, M. A. QM/MMpol: A consistent model for solute/ solvent polarization. application to the aqueous solvation and spectroscopy of formaldehyde, acetaldehyde, and acetone. J. Phys. Chem. 1996, 100, 14492−14507. (44) Gao, J.; Byun, K. Solvent effects on the n→ π* transition of pyrimidine in aqueous solution. Theor. Chem. Acc. 1997, 96, 151−156. (45) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: dielectric, electostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227−249. (46) Li, H. Quantum mechanical/molecular mechanical/continuum style solvation model: Linear response theory, variational treatment, and nuclear gradients. J. Chem. Phys. 2009, 131, 184103. (47) Steindal, A. H.; Ruud, K.; Frediani, L.; Aidas, K.; Kongsted, J. Excitation energies in solution: The fully polarizable QM/MM/PCM method. J. Phys. Chem. B 2011, 115, 3027−3037. (48) Lipparini, F.; Cappelli, C.; Barone, V. Linear response theory and electronic transition energies for a fully polarizable QM/classical Hamiltonian. J. Chem. Theory Comput. 2012, 8, 4153−4165. (49) Sneskov, K.; Schwabe, T.; Christiansen, O.; Kongsted, J. Scrutinizing the effects of polarization in QM/MM excited state calculations. Phys. Chem. Chem. Phys. 2011, 13, 18551−18560. (50) Schwabe, T.; Olsen, J. M. H.; Sneskov, K.; Kongsted, J.; Christiansen, O. Solvation effects on electronic transitions: exploring the performance of advanced solvent potentials in polarizable embedding calculations. J. Chem. Theory Comput. 2011, 7, 2209−2217. (51) Murugan, N. A.; Kongsted, J.; Rinkevicius, Z.; Ågren, H. Demystifying the solvatochromic reversal in Brookers merocyanine dye. Phys. Chem. Chem. Phys. 2011, 13, 1290−1292. (52) Arul Murugan, N.; Kongsted, J.; Rinkevicius, Z.; Aidas, K.; Mikkelsen, K. V.; Ågren, H. Hybrid density functional theory/ molecular mechanics calculations of two-photon absorption of dimethylamino nitro stilbene in solution. Phys. Chem. Chem. Phys. 2011, 13, 12506−12516. (53) Silva, D. L.; Murugan, N. A.; Kongsted, J.; Rinkevicius, Z.; Canuto, S.; Ågren, H. The role of molecular conformation and polarizable embedding for one- and two-photon absorption of disperse orange 3 in solution. J. Phys. Chem. B 2012, 116, 8169−8181. (54) Olesiak-Banska, J.; Matczyszyn, K.; Zaleśny, R.; Murugan, N. A.; Kongsted, J.; Ågren, H.; Bartkowiak, W.; Samoc, M. Revealing spectral features in two-photon absorption spectrum of Hoechst 33342: A combined experimental and quantum-chemical study. J. Phys. Chem. B 2013, 117, 12013−12019. (55) Wielgus, M.; Zaleśny, R.; Murugan, N. A.; Kongsted, J.; Ågren, H.; Samoc, M.; Bartkowiak, W. Two-photon solvatochromism II: Experimental and theoretical study of solvent effects on the twophoton absorption spectrum of Reichardtÿ dye. ChemPhysChem 2013, 14, 3731−3739. (56) Olsen, J. M. H.; Kongsted, J. In Advances in Quantum Chemistry; Sabin, J. R., Brändas, E., Eds.; Academic Press: Waltham, MA, 2011; Vol. 61, pp 107−143.

(57) Steindal, A. H.; Olsen, J. M. H.; Ruud, K.; Frediani, L.; Kongsted, J. A combined quantum mechanics/molecular mechanics study of the one- and two-photon absorption in the green fluorescent protein. Phys. Chem. Chem. Phys. 2012, 14, 5440−5451. (58) Sneskov, K.; Olsen, J. M. H.; Schwabe, T.; Hättig, C.; Christiansen, O.; Kongsted, J. Computational screening of one- and two-photon spectrally tuned channelrhodopsin mutants. Phys. Chem. Chem. Phys. 2013, 15, 7567−7576. (59) List, N. H.; Pimenta, F. M.; Holmegaard, L.; Jensen, R. L.; Etzerodt, M.; Schwabe, T.; Kongsted, J.; Ogilby, P. R.; Christiansen, O. Effect of chromophore encapsulation on linear and nonlinear optical properties: the case of “miniSOG”, a protein-encased flavin. Phys. Chem. Chem. Phys. 2014, 16, 9950−9959. (60) Beerepoot, M. T. P.; Steindal, A. H.; Kongsted, J.; Brandsdal, B. O.; Frediani, L.; Ruud, K.; Olsen, J. M. H. A polarizable embedding DFT study of one-photon absorption in fluorescent proteins. Phys. Chem. Chem. Phys. 2013, 15, 4735−4743. (61) Beerepoot, M. T. P.; Steindal, A. H.; Ruud, K.; Olsen, J. M. H.; Kongsted, J. Convergence of environment polarization effects in multiscale modeling of excitation energies. Comput. Theor. Chem. 2014, 1040-1041, 304−311. (62) Schwabe, T.; Beerepoot, M. T. P.; Olsen, J. M. H.; Kongsted, J. Analysis of computational models for an accurate study of electronic excitations in GFP. Phys. Chem. Chem. Phys. 2015, 17, 2582−2588. (63) Car, R.; Parrinello, M. Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (64) Seifert, G.; Joswig, J.-O. Density-functional tight binding−an approximate density-functional theory method. WIREs Comput. Mol. Sci. 2012, 2, 456−465. (65) Halgren, T. A.; Damm, W. Polarizable force fields. Curr. Opin. Struct. Biol. 2001, 11, 236−242. (66) Cisneros, G. A.; Karttunen, M.; Ren, P.; Sagui, C. Classical electrostatics for biomlecular simulations. Chem. Rev. 2014, 114, 779− 814. (67) Ren, P.; Ponder, J. W. Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B 2003, 107, 5933− 5947. (68) Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. A simple polarizable model of water based on classical Drude oscillators. J. Chem. Phys. 2003, 119, 5185−5198. (69) Yu, H.; Geerke, D. P.; Liu, H.; van Gunsteren, W. F. Molecular dynamics simulations of liquid methanol and methanol-water mixtures with polarizable models. J. Comput. Chem. 2006, 27, 1494−1504. (70) Bachmann, S. J.; van Gunsteren, W. F. An improved simple polarisable water model for use in biomolecular simulation. J. Chem. Phys. 2014, 141, 22D515. (71) Grimme, S. A general quantum mechanically derived force field (QMDFF) for molecules and condensed phase simulations. J. Chem. Theory Comput. 2014, 10, 4497−4514. (72) Korth, M. Empirical hydrogen-bond potential functions−an old hat reconditioned. ChemPhysChem 2011, 12, 3131−3142. (73) Patel, S.; Gorai, S.; Malik, P. K. Preferential solvation through selective functional group recognition in p-nitroaniline. J. Photochem. Photobiol., A 2011, 219, 76−83. (74) Sok, S.; Willow, S. Y.; Zahariev, D.; Gordon, M. S. Solventinduced shift of the lowest singlet π → π* charge-transfer excited state of p-nitroaniline in water: An application of the TDDFT/EFP1 method. J. Phys. Chem. A 2011, 115, 9801−9809. (75) Liang, W.; Chapman, C. T.; Ding, F.; Li, X. Modeling ultrafast solvated electronic dynamics using time-dependent density functional theory and polarizable continuum model. J. Phys. Chem. A 2012, 116, 1884−1890. (76) Frutos-Puerto, S.; Aguilar, M. A.; Galván, I. F. Theoretical study of the preferential solvation effect on the solvatochromic shifts of paranitroaniline. J. Phys. Chem. B 2013, 117, 2466−2474. (77) Eriksen, J. J.; Sauer, S. P. A.; Mikkelsen, K. V.; Christiansen, O.; Jensen, H. J. A.; Kongsted, J. Failures of TDDFT in describing the lowest intramolecular charge-transfer excitation in para-nitroaniline. Mol. Phys. 2013, 111, 1235−1248. 10699


Article


The Journal of Physical Chemistry B (78) Hidalgo, M.; Rivelino, R.; Canuto, S. Origin of the red shift for the lowest singlet π → π* charge-transfer absorption of p-aitroaniline in supercritical CO2. J. Chem. Theory Comput. 2014, 10, 1554−1562. (79) Piacente, G.; Aschi, M.; Cerichelli, G.; Chiarini, M.; Amadei, A.; D’Aiuto, G. P. A. V. Inclusion of cybotactic effect in the theoretical modeling of absorption spectra of liquid-state systems with perturbed matrix method and molecular dynamics simulations: the UV-Vis absorption spectrum of para-nitroaniline as a case study. Theor. Chem. Acc. 2014, 133, 1478−1487. (80) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (81) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density dunctional force fields. J. Phys. Chem. 1994, 98, 11623−11627. (82) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (83) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796−6806. (84) Cancés, E.; Mennucci, B.; Tomasi, J. A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics. J. Chem. Phys. 1997, 107, 3032−3031. (85) Mennucci, B.; Cancés, E.; Tomasi, J. Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: Theoretical bases, computational implementation, and numerical applications. J. Phys. Chem. B 1997, 101, 10506−10517. (86) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, revision A.02, Gaussian, Inc.: Wallingford CT, 2009. (87) Gagliardi, L.; Lindh, R.; Karlström, G. Local properties of quantum chemical systems: The LoProp approach. J. Chem. Phys. 2004, 121, 4494−4500. (88) Aquilante, F.; de Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P.A.; Neogrády, P.; Pedersen, T. B.; Pitoňaḱ , M.; Reiher, M.; Roos, B. O.; et al. MOLCAS 7: The next generation. J. Comput. Chem. 2010, 31, 224−247. (89) Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11, 361−373. (90) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (91) Ponder, J. W. TINKER: Software tools for molecular design, 5.1; Washington University School of Medicine: Saint Louis, MO, 2010. (92) Fine, R. A.; Millero, F. J. Compressibility of water as a function of temperature and pressure. J. Chem. Phys. 1973, 59, 5529−5537. (93) Diaz Peña, M.; Tardajos, G. Isothermal compressibilities of n-1alcohols from methanol to 1-dodecanol at 298.15, 308.15, and 333.15 K. J. Chem. Thermodyn. 1979, 11, 441−445. (94) Nath, J.; Dixit, P. Ultrasonic velocities in, and adiabatic compressibilities and excess volumes for, binary liquid mixtures of acetone with tetrachloroethylene, trichloroethylene, methylene chloride, 1,2-dichloroethane, and cyclohexane. J. Chem. Eng. Data 1984, 29, 313−316. (95) Linse, P. MOLSIM 5.2; Lund University: Lund, Sweden, 2011. (96) Aidas, K. Whirlpool, a QM/MM analysis program, 1.0, 2010. (97) Weigend, F.; Häser, M. RI-MP2: First derivatives and global consistency. Theor. Chem. Acc. 1997, 97, 331−340. (98) Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116, 3175−3183.

(99) Häser, M.; Ahlrichs, R. Improvements on the direct SCF method. J. Comput. Chem. 1989, 10, 104−111. (100) Christiansen, O.; Koch, H.; Jørgensen, P. The second-order approximate coupled cluster singles and doubles model CC2. Chem. Phys. Lett. 1995, 243, 409−418. (101) Hättig, C.; Weigend, F. CC2 excitation energy calculations on large molecules using the resolution of the identity approximation. J. Chem. Phys. 2000, 113, 5154−5161. (102) TURBOMOLE V6.6, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007; TURBOMOLE GmbH: Karlsruhe, Germany, 2014. (103) Hydrogen Bonding: New Insights; Grabowski, S. J., Ed.; Springer: Netherlands, 2006. (104) Weinhold, F.; Klein, R. A. What is a hydrogen bond? Mutually consistent theoretical and experimental criteria for characterizing Hbonding interactions. Mol. Phys. 2012, 110, 565−579.

10700


Assessing Molecular Dynamics Simulations with Solvatochromism

Recommend Documents