A new method to locate saddle points for reactions in solution by using the free-energy gradient method and the mean field approximation

A new method for calculating saddle points of reactions in solution is presented. The main characteristics of the method are: (1) the solute-solvent system is described by the averaged solvent electrostatic potential/molecular dynamics method (ASEP/MD). This is a quantum mechanics/molecular mechanics method (QM/MM) that makes use of the mean field approximation (MFA) and that permits one to simultaneously optimize the electronic structure and geometry of the solute molecule and the solvent structure around it. (2) The transition state is located by the joint use of the free-energy gradient method and the mean field approximation. An application to the study of the Menshutkin reaction between NH(3) and CH(3)Cl in aqueous solution is discussed. The accuracy and usefulness of the proposed method is checked through comparison with other methods.     

An averaged solvent electrostatic potential from molecular dynamics study of the anomeric equilibrium of D-xylose in aqueous solution

We applied the free-energy perturbation method together with the averaged solvent electrostatic potential from molecular dynamics (ASEP/MD) method to study the anomeric equilibrium of d-xylose in aqueous solution. The level of calculation, 6-311G++(2d,2p) basis set and density functional theory, permits one to explain the main characteristics of the anomeric equilibrium of d-xylopyranose: in vacuo, the anomeric effect predominates and the form is the stabler. In water, solvation leads to the form being the stabler. A comparison between the performances of the ASEP/MD and polarizable continuum models is also presented.Contribution to the Jacopo Tomasi Honorary Issue     

Theoretical study of the relative stability of rotational conformers of alpha and beta-D-glucopyranose in gas phase and aqueous solution

The alpha-beta anomer energy difference and the stability of 10 rotamers of counterclockwise D-glucopyranose were studied in vacuo and in aqueous solution at the B3LYP/6-31+G(d,p) level. To obtain the solute charge distribution and the solvent structure around it, we used the averaged solvent electrostatic potential from molecular dynamics method, ASEP/MD, which alternates molecular dynamics and quantum mechanics calculations in an iterative procedure. The main characteristics of the anomeric equilibrium, both in vacuo and in solution, are well reproduced. The relative stability of the different anomers is related to the availability of the free pairs of electrons in the anomeric oxygen to interact with the water molecules. The influence of solvation in the conformer equilibrium is also analyzed.     

Theoretical study of the 1,3-hydrogen shift of triazene in water

The 1,3-hydrogen shift of triazene in aqueous solution was studied with a combination of QM/MM methods. First, the different species involved were characterized and the activation free-energies calculated with ASEP/MD, a method that makes use of the mean field approximation. Then the reaction dynamics was simulated with a QM/MM/MD method. A very strong influence of the solvent was observed, both specific, with the participation of a water molecule, and from the rest of the solvent. The effect of solvation on the geometry and electron distribution of triazene is important: N-N bond lengths tend to be more similar and the molecule acquires a planar structure. For the transition state structure, a substantial degree of ionic nature was found. Dynamic solvent effects were also analyzed.     

Nitrile and thiocyanate IR probes: molecular dynamics simulation studies

Nitrile- and thiocyanate-derivatized amino acids have been found to be useful IR probes for investigating their local electrostatic environments in proteins. To shed light on the CN stretch frequency shift and spectral lineshape change induced by interactions with hydrogen-bonding solvent molecules, we carried out both classical and quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations for MeCN and MeSCN in water. These QM/MM and conventional force field MD simulation results were found to be inconsistent with the experimental results as well as with the high-level ab initio calculation results of MeCN-water and MeSCN-water potential energies. Thus, a new set of atomic partial charges of MeCN and MeSCN is obtained. By using the MD simulation trajectories and the electrostatic potential model recently developed, the CN and SCN stretching mode frequency trajectories were obtained and used to simulate the IR spectra. The C[Triple Bond]N frequency blueshifts of MeCN and MeSCN in water are estimated to be 9.0 and 1.9 cm(-1), respectively, in comparison with those of gas phase values. These values are found to be in reasonable agreement with the experimentally measured IR spectra of MeCN, MeSCN, beta-cyano-L-alanine, and cyanylated cysteine in water and other polar solvents.     

Solvent effects on internal conversions and intersystem crossings: the radiationless de-excitation of acrolein in water

An extended version of the averaged solvent electrostatic potential from molecular dynamics data (ASEP/MD) method oriented to the study of the solvent effects on internal conversion and intersystem crossing processes is presented. The method allows for the location of crossing points between free energy surfaces both in equilibrium and in frozen solvent conditions. The ground and excited states of the solute molecule are described at the complete active space self-consistent field (CASSCF) level while the solvent structure is obtained from molecular dynamics simulations. As an application, we studied the nonradiative de-excitation of s-trans-acrolein 1(n --> pi*) in aqueous solution. We found that the solvent modifies the relative stability of the different crossing points but not enough as to alter the relative order of stability with respect to the in vacuo situation. The relaxation through an equilibrium path involves a strong solvent reorganization. On the contrary, the nonequilibrium path does not involve solvent motion and the de-excitation could proceed with the same speed as in vacuo.     

Strategies to model the near-solute solvent molecular density/polarization

The solvent molecular distribution significantly affects the behavior of the solute molecules and is thus important in studying many biological phenomena. It can be described by the solvent molecular density distribution, g, and the solvent electric dipole distribution, p. The g and p can be computed directly by counting the number of solvent molecules/dipoles in a microscopic volume centered at r during a simulation or indirectly from the mean force F and electrostatic field E acting on the solvent molecule at r, respectively. However, it is not clear how the g and p derived from simulations depend on the solvent molecular center or the solute charge and if the g(F) and p(E) computed from the mean force and electric field acting on the solvent molecule, respectively, could reproduce the corresponding g and p obtained by direct counting. Hence, we have computed g, p, g(F), and p(E) using different water centers from simulations of a solute atom of varying charge solvated in TIP3P water. The results show that g(F) and p(E) can reproduce the g and p obtained using a given count center. This implies that rather than solving the coordinates of each water molecule by MD simulations, the distribution of water molecules could be indirectly obtained from analytical formulas for the mean force F and electrostatic field E acting on the solvent molecule at r. Furthermore, the dependence of the g and p distributions on the solute charge revealed provides an estimate of the change in g and p surrounding a biomolecule upon a change in its conformation.     

Communication: Quantum polarized fluctuating charge model: a practical method to include ligand polarizability in biomolecular simulations

We present a simple and practical method to include ligand electronic polarization in molecular dynamics (MD) simulation of biomolecular systems. The method involves periodically spawning quantum mechanical (QM) electrostatic potential (ESP) calculations on an extra set of computer processors using molecular coordinate snapshots from a running parallel MD simulation. The QM ESPs are evaluated for the small-molecule ligand in the presence of the electric field induced by the protein, solvent, and ion charges within the MD snapshot. Partial charges on ligand atom centers are fit through the multi-conformer restrained electrostatic potential (RESP) fit method on several successive ESPs. The RESP method was selected since it produces charges consistent with the AMBER/GAFF force-field used in the simulations. The updated charges are introduced back into the running simulation when the next snapshot is saved. The result is a simulation whose ligand partial charges continuously respond in real-time to the short-term mean electrostatic field of the evolving environment without incurring additional wall-clock time. We show that (1) by incorporating the cost of polarization back into the potential energy of the MD simulation, the algorithm conserves energy when run in the microcanonical ensemble and (2) the mean solvation free energies for 15 neutral amino acid side chains calculated with the quantum polarized fluctuating charge method and thermodynamic integration agree better with experiment relative to the Amber fixed charge force-field.     

Molecular simulations of outersphere reorganization energies in polar and quadrupolar solvents. The case of intramolecular electron and hole transfer

Outersphere reorganization energies (lambda) for intramolecular electron and hole transfer are studied in anion- and cation-radical forms of complex organic substrates (p-phenylphenyl-spacer-naphthyl) in polar (water, 1,2-dichloroethane, tetrahydrofuran) and quadrupolar (supercritical CO2) solvents. Structure and charge distributions of solute molecules are obtained at the HF/6-31G(d,p) level. Standard Lennard-Jones parameters for solutes and the nonpolarizable simple site-based models of solvents are used in molecular dynamics (MD) simulations. Calculation of lambda is done by means of the original procedure, which treats electrostatic polarization of a solvent in terms of a usual nonpolarizable MD scheme supplemented by scaling of reorganization energies at the final stage. This approach provides a physically relevant background for separating inertial and inertialless polarization responses by means of a single parameter epsilon(infinity), optical dielectric permittivity of the solvent. Absolute lambda values for hole transfer in 1,2-dichloroethane agree with results of previous computations in terms of the different technique (MD/FRCM, Leontyev, I. V.; et al. Chem. Phys. 2005, 319, 4). Computed lambda values for electron transfer in tetrahydrofuran are larger than the experimental values by ca. 2.5 kcal/mol; for the case of hole transfer in 1,2-dichloroethane the discrepancy is of similar magnitude provided the experimental data are properly corrected. The MD approach gives nonzero lambda values for charge-transfer reaction in supercritical CO2, being able to provide a uniform treatment of nonequilibrium solvation phenomena in both quadrupolar and polar solvents.     

Solvent effect on the absorption spectra of coumarin 120 in water: A combined quantum mechanical and molecular mechanical study

The solvent effect on the absorption spectra of coumarin 120 (C120) in water was studied utilizing the combined quantum mechanical∕molecular mechanical (QM∕MM) method. In molecular dynamics (MD) simulation, a new sampling scheme was introduced to provide enough samples for both solute and solvent molecules to obtain the average physical properties of the molecules in solution. We sampled the structure of the solute and solvent molecules separately. First, we executed a QM∕MM MD simulation, where we sampled the solute molecule in solution. Next, we chose random solute structures from this simulation and performed classical MD simulation for each chosen solute structure with its geometry fixed. This new scheme allowed us to sample the solute molecule quantum mechanically and sample many solvent structures classically. Excitation energy calculations using the selected samples were carried out by the generalized multiconfigurational perturbation theory. We succeeded in constructing the absorption spectra and realizing the red shift of the absorption spectra found in polar solvents. To understand the motion of C120 in water, we carried out principal component analysis and found that the motion of the methyl group made the largest contribution and the motion of the amino group the second largest. The solvent effect on the absorption spectrum was studied by decomposing it in two components: the effect from the distortion of the solute molecule and the field effect from the solvent molecules. The solvent effect from the solvent molecules shows large contribution to the solvent shift of the peak of the absorption spectrum, while the solvent effect from the solute molecule shows no contribution. The solvent effect from the solute molecule mainly contributes to the broadening of the absorption spectrum. In the solvent effect, the variation in C-C bond length has the largest contribution on the absorption spectrum from the solute molecule. For the solvent effect on the absorption spectrum from the solvent molecules, the solvent structure around the amino group of C120 plays the key role.     

Theoretical calculation of the coiled-coil stability in water in the context of its possible use as a molecular rack

The coiled-coil stability and rigidity may be of importance for molecular electronics (electronically bistable molecules). The coiled-coil binding free energy has been calculated using molecular dynamics (MD). The energy has been computed as a difference of the appropriate free energies; derived for the coiled-coil and isolated alpha-helices separately. All MD simulations have been performed using an explicit model of the solvent, whereas the continuum solvent approach has been applied to analyze the MD trajectories. The computed stability of the coiled-coil is of the order of -87 kcal/mol, i.e., about -1.2 kcal/mol per amino acid residue, and arises mainly from the electrostatic interactions and hydrophobic effect. The entropy term has been roughly estimated to be of the order of -22 kcal/mol. This assures that coiled-coil polypeptides may be used as a stable molecular scaffolding.     

Charge separation reaction in clusters of polar molecules: MD simulations

Rate constant of intermolecular electron transfer (ET) in a photoexcited donor-acceptor model system solvated by a cluster of polar molecules has been expressed in terms of the statistical distribution of the electrostatic potential energy difference between the reacting sites. This distribution has been calculated for a particular case of acetonitrile clusters a ≈120 K by MD computer simulation. The MD values of the cluster reorganization energy and the ET rate constant have been compared with the corresponding MD results for the donor-acceptor pair solvated in bulk acetonitrile and with theoretical predictions based on the continuum model.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Effects of solute structure on local solvation and solvent interactions: results from UV/vis spectroscopy and molecular dynamics simulations

Solvation of heterocyclic amines in CO(2)-expanded methanol (MeOH) has been explored with UV/vis spectroscopy and molecular dynamics (MD) simulations. A synergistic study of experiments and simulations allows exploration of solute and solvent effects on solvation and the molecular interactions that affect absorption. MeOH-nitrogen hydrogen bonds hinder the n-pi* transition; however, CO(2) addition causes a blue shift relative to MeOH because of Lewis acid/base interactions with nitrogen. Effects of solute structure are considered, and very different absorption spectra are obtained as nitrogen positions change. MD simulations provide detailed solvent clustering behavior around the solute molecules and show that the local solvent environment and ultimately the spectra are sensitive to the solute structure. This work demonstrates the importance of atomic-level information in determining the structure-property relationships between solute structure, local salvation, and solvatochromism.     

The use of complexation induced proton NMR chemical shifts for structural analysis of host-guest complexes in solution

Proton shielding variations in supramolecular complexes contain a wealth of information on complex geometries in solution that has been until now mostly neglected. We describe herein ways for such analyses with five cyclophane and two cyclodextrin complexes in water, by using a program SHIFT which is based on and parametrized with the analyses of over 300 intramolecular proton shift variations in well defined molecular frameworks such as steroids or cyclophanes. The intermolecular shift changes in the host-guest complexes at 100% complexation (CIS values) are calculated as sum of anisotropy effects deltachi from aromatic ring currents and linear electric field effects LEF, based on force field generated geometries. The conformations with the best agreement between calculated and observed CIS values are at least for noncharged guest compounds close to those obtained from molecular mechanics and/or MD calculations and intermolecular NOEs (where available), noticeably without adjusting the complex geometries to the experimental CIS. Through-space electrostatic field effects LEF, which have been until now often neglected, can be sizeable also for noncharged systems; best agreement between experiment and calculation is observed with Gasteiger atomic charges.     

Fluorescence spectra of organic dyes in solution: a time dependent multilevel approach

Classical all-atom molecular dynamics (MD) simulations and quantum mechanical (QM) time-dependent density functional theory (TD-DFT) calculations are employed to study the conformational and photophysical properties of the first emitter excited state of tetramethyl-rhodamine iso-thiocyanate fluorophore in aqueous solution. For this purpose, a specific and accurate force field has been parameterised from QM data to model the fluorophore's first bright excited state. During the MD simulations, the consequences of the π→π* electronic transition on the structure and microsolvation sphere of the dye has been analysed in some detail and compared to the ground state behaviour. Thereafter, fluorescence has been calculated at the TD-DFT level on configurations sampled from the simulated MD trajectories, allowing us to include time dependent solvent effects in the computed emission spectrum. The latter, when compared with the absorption spectrum, reproduces well the experimental Stokes shift, further validating the proposed multilevel computational procedure.