Reliable molecular simulations of solute-solvent systems with a minimum number of solvent shells

Computational Spectroscopy in Solution by Integration of Variational and Perturbative Approaches on Top of Clusterized Molecular Dynamics

Multiscale QM/MM approaches have become the most suitable and effective methods for the investigation of spectroscopic properties of medium- or large-size chromophores in condensed phases. On these grounds, we are developing a novel workflow aimed at improving the generality, reliability, and ease of use of the available tools. In the present paper, we report the latest developments of such an approach with specific reference to a general workplan starting with the addition of acetonitrile to the panel of solvents already available in the General Liquid Optimized Boundary (GLOB) model enforcing nonperiodic boundary conditions (NPBC). Next, the solvatochromic shifts induced by acetonitrile on both rigid (uracil and thymine) and flexible (thyrosine) chromophores have been studied introducing in our software a number of new features ranging from rigid-geometry NPBC molecular dynamics based on the quaternion formalism to a full integration of variational (ONIOM) and perturbative (perturbed matrix method (PMM)) approaches for describing different solute–solvent topologies and local fluctuations, respectively. Finally, thymine and uracil have been studied also in methanol to point out the generality of the computational strategy. While further developments are surely needed, the strengths of our integrated approach even in its present version are demonstrated by the accuracy of the results obtained by an unsupervised approach and coupled to a computational cost strongly reduced with respect to that of conventional QM/MM models without any appreciable accuracy deterioration.

Multiscale Simulation of Protein Hydration Using the SWINGER Dynamical Clustering Algorithm

To perform computationally efficient concurrent multiscale simulations of biological macromolecules in solution, where the all-atom (AT) models are coupled to supramolecular coarse-grained (SCG) solvent models, previous studies resorted to modified AT water models, such as the bundled-simple point charge (SPC) models, that use semiharmonic springs to restrict the relative movement of water molecules within a cluster. Those models can have a significant impact on the simulated biomolecules and can lead, for example, to a partial unfolding of a protein. In this work, we employ the recently developed alternative approach with a dynamical clustering algorithm, SWINGER, which enables a direct coupling of original unmodified AT and SCG water models. We perform an adaptive resolution molecular dynamics simulation of a Trp-Cage miniprotein in multiscale water, where the standard SPC water model is interfaced with the widely used MARTINI SCG model, and demonstrate that, compared to the corresponding full-blown AT simulations, the structural and dynamic properties of the solvated protein and surrounding solvent are well reproduced by our approach.

Polarizable molecular dynamics in a polarizable continuum solvent

We present, for the first time, scalable polarizable molecular dynamics (MD) simulations within a polarizable continuum solvent with molecular shape cavities and exact solution of the mutual polarization. The key ingredients are a very efficient algorithm for solving the equations associated with the polarizable continuum, in particular, the domain decomposition Conductor-like Screening Model (ddCOSMO), which involves a rigorous coupling of the continuum with the polarizable force field achieved through a robust variational formulation and an effective strategy to solve the coupled equations. The coupling of ddCOSMO with nonvariational force fields, including AMOEBA, is also addressed. The MD simulations are feasible, for real-life systems, on standard cluster nodes; a scalable parallel implementation allows for further acceleration in the context of a newly developed module in Tinker, named Tinker-HP. NVE simulations are stable, and long-term energy conservation can be achieved. This paper is focused on the methodological developments, the analysis of the algorithm, and the stability of the simulations; a proof-of-concept application is also presented to attest to the possibilities of this newly developed technique.

Elastic Barrier Dynamical Freezing in Free Energy Calculations: A Way To Speed Up Nonequilibrium Molecular Dynamics Simulations by Orders of Magnitude

An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called "elastic barrier dynamical freezing", according to which a stiff potential-energy barrier is enforced at the boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speedup exceeding a few orders of magnitude.

Boundary condition effects on the dynamic and electric properties of hydration layers

Water solvation has a central role in several biochemical processes ranging from protein folding to biomolecular recognition and enzyme catalysis. Because of its importance, the structure and dynamics of hydration layers around biological macromolecules have been the targets of a great number of experimental and computational studies. In the present contribution, we have investigated the effects of periodic boundary conditions (PBCs), as used in conjunction with molecular dynamics (MD) simulations, on the dynamic and electric properties of water layers. In particular, we have systematically performed MD simulations of neat water and biomolecules in aqueous solutions by imposing a different external dielectric constant, a generally overlooked parameter in PBC simulations. The effect of the system size has also been addressed. Overall, our results consistently indicate that the dipole moment properties of water layers, and specifically the dipole moment fluctuations and the reorientational correlation functions, can be sensitive to the choice of the external boundary conditions, whereas other molecular properties, such as the self-diffusion coefficient and the reorientational relaxation times, are not affected. We think that our investigation may help to assess appropriate simulation conditions for modeling the aqueous environment of relevant biochemical systems and processes.

Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II.Towards Massively Parallel Computations using Smooth Particle Mesh Ewald

In this paper, we present a scalable and efficient implementation of point dipole-based polarizable force fields for molecular dynamics (MD) simulations with periodic boundary conditions (PBC). The Smooth Particle-Mesh Ewald technique is combined with two optimal iterative strategies, namely, a preconditioned conjugate gradient solver and a Jacobi solver in conjunction with the Direct Inversion in the Iterative Subspace for convergence acceleration, to solve the polarization equations. We show that both solvers exhibit very good parallel performances and overall very competitive timings in an energy-force computation needed to perform a MD step. Various tests on large systems are provided in the context of the polarizable AMOEBA force field as implemented in the newly developed Tinker-HP package which is the first implementation for a polarizable model making large scale experiments for massively parallel PBC point dipole models possible. We show that using a large number of cores offers a significant acceleration of the overall process involving the iterative methods within the context of spme and a noticeable improvement of the memory management giving access to very large systems (hundreds of thousands of atoms) as the algorithm naturally distributes the data on different cores. Coupled with advanced MD techniques, gains ranging from 2 to 3 orders of magnitude in time are now possible compared to non-optimized, sequential implementations giving new directions for polarizable molecular dynamics in periodic boundary conditions using massively parallel implementations.

Combining the Fluctuating Charge Method, Non-Periodic Boundary Conditions and Meta-Dynamics: Aqua Ions as case studies

We present the current status of development of our code for performing Molecular Dynamics (MD) simulations exploiting a polarizable force field based on the Fluctuating Charge (FQ) method and non-Periodic Boundary Conditions (NPBC). Continuing on the path set in a previous work, we increased the capabilities of the code by implementing a number of new features, including: a non-iterative algorithm for rigid trigonal molecule simulations; two additional temperature coupling schemes; a meta-dynamics based approach for effective free energy evaluations. Although these are well known algorithms, each present in one or more widely used MD packages, they have now been tested, for the first time, in the context of the FQ model coupled with NPBC. As case studies, we considered three aqueous ions of increasing charge, namely Na⁺, Ca²⁺ and La³⁺, at infinite dilution. In particular, by exploiting a computational approach recently proposed by our group and based on the metadynamics technique, we focused on the important role played by solvent polarization on ionic hydration structures, also investigating the free energy landscapes of ion coordination and the water exchange rates. Such an approach, previously tested with standard non-polarizable models, was applied here to evaluate the effects of explicit polarization on water exchange barriers between different solvent coordination structures. Moreover, we have analyzed and discussed in some detail non-linear electrostatic effects arising from solvent polarization while going from a mono- to a di- and trivalent ion.

QM/MM investigations of organic chemistry oriented questions

About 35 years after its first suggestion, QM/MM became the standard theoretical approach to investigate enzymatic structures and processes. The success is due to the ability of QM/MM to provide an accurate atomistic picture of enzymes and related processes. This picture can even be turned into a movie if nuclei-dynamics is taken into account to describe enzymatic processes. In the field of organic chemistry, QM/MM methods are used to a much lesser extent although almost all relevant processes happen in condensed matter or are influenced by complicated interactions between substrate and catalyst. There is less importance for theoretical organic chemistry since the influence of nonpolar solvents is rather weak and the effect of polar solvents can often be accurately described by continuum approaches. Catalytic processes (homogeneous and heterogeneous) can often be reduced to truncated model systems, which are so small that pure quantum-mechanical approaches can be employed. However, since QM/MM becomes more and more efficient due to the success in software and hardware developments, it is more and more used in theoretical organic chemistry to study effects which result from the molecular nature of the environment. It is shown by many examples discussed in this review that the influence can be tremendous, even for nonpolar reactions. The importance of environmental effects in theoretical spectroscopy was already known. Due to its benefits, QM/MM can be expected to experience ongoing growth for the next decade.In the present chapter we give an overview of QM/MM developments and their importance in theoretical organic chemistry, and review applications which give impressions of the possibilities and the importance of the relevant effects. Since there is already a bunch of excellent reviews dealing with QM/MM, we will discuss fundamental ingredients and developments of QM/MM very briefly with a focus on very recent progress. For the applications we follow a similar strategy.

Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system or box size) and treatment of electrostatic interactions (Coulombic, lattice-sum, or cutoff-based) used during these simulations. However, as shown by Kastenholz and Hünenberger [J. Chem. Phys. 124, 224501 (2006)], correction terms can be derived for the effects of: (A) an incorrect solvent polarization around the ion and an incomplete or/and inexact interaction of the ion with the polarized solvent due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite-size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site, and the possible presence of a polarized air-liquid interface or of a constraint of vanishing average electrostatic potential in the simulated system; and (D) an inaccurate dielectric permittivity of the employed solvent model. Comparison with standard experimental data also requires the inclusion of appropriate cavity-formation and standard-state correction terms. In the present study, this correction scheme is extended by: (i) providing simple approximate analytical expressions (empirically-fitted) for the correction terms that were evaluated numerically in the above scheme (continuum-electrostatics calculations); (ii) providing correction terms for derivative thermodynamic single-ion solvation properties (and corresponding partial molar variables in solution), namely, the enthalpy, entropy, isobaric heat capacity, volume, isothermal compressibility, and isobaric expansivity (including appropriate standard-state correction terms). The ability of the correction scheme to produce methodology-independent single-ion solvation free energies based on atomistic simulations is tested in the case of Na(+) hydration, and the nature and magnitude of the correction terms for derivative thermodynamic properties is assessed numerically.

Molecular dynamics simulations of a reversibly folding beta-heptapeptide in methanol: influence of the treatment of long-range electrostatic interactions

Eight 100-ns molecular dynamics simulations of a beta-heptapeptide in methanol at 340 K (within cubic periodic computational boxes of about 6-nm edge) are reported and compared. These simulations were performed with three different charge-state combinations at the peptide termini, one of them with or without a neutralizing chloride counterion, and using either the lattice-sum (LS) or reaction-field (RF) scheme to handle electrostatic interactions. The choice of the electrostatic scheme has essentially no influence on the folding-unfolding equilibrium when the peptide termini are uncharged and only a small influence when the peptide is positively charged at its N-terminus (with or without inclusion of a neutralizing chloride counterion). However, when the peptide is zwitterionic, the LS scheme leads to preferential sampling of the high-dipole folded helical state, whereas the RF scheme leads to preferential sampling of a low-dipole unfolded salt-bridged state. A continuum electrostatics analysis based on the sampled configurations (zwitterionic case) suggests that the LS scheme stabilizes the helical state through artificial periodicity, but that the magnitude of this perturbation is essentially negligible (compared to the thermal energy) for the large box size and relatively polar solvent considered. The results thus provide clear evidence (continuum electrostatics analysis) for the absence of LS artifacts and some indications (still not definitive because of the limited sampling of the folding-unfolding transition) for the presence of RF artifacts in this specific system.

Interpreting the Coulomb-field approximation for generalized-Born electrostatics using boundary-integral equation theory

The importance of molecular electrostatic interactions in aqueous solution has motivated extensive research into physical models and numerical methods for their estimation. The computational costs associated with simulations that include many explicit water molecules have driven the development of implicit-solvent models, with generalized-Born (GB) models among the most popular of these. In this paper, we analyze a boundary-integral equation interpretation for the Coulomb-field approximation (CFA), which plays a central role in most GB models. This interpretation offers new insights into the nature of the CFA, which traditionally has been assessed using only a single point charge in the solute. The boundary-integral interpretation of the CFA allows the use of multiple point charges, or even continuous charge distributions, leading naturally to methods that eliminate the interpolation inaccuracies associated with the Still equation. This approach, which we call boundary-integral-based electrostatic estimation by the CFA (BIBEE/CFA), is most accurate when the molecular charge distribution generates a smooth normal displacement field at the solute-solvent boundary, and CFA-based GB methods perform similarly. Conversely, both methods are least accurate for charge distributions that give rise to rapidly varying or highly localized normal displacement fields. Supporting this analysis are comparisons of the reaction-potential matrices calculated using GB methods and boundary-element-method (BEM) simulations. An approximation similar to BIBEE/CFA exhibits complementary behavior, with superior accuracy for charge distributions that generate rapidly varying normal fields and poorer accuracy for distributions that produce smooth fields. This approximation, BIBEE by preconditioning (BIBEE/P), essentially generates initial guesses for preconditioned Krylov-subspace iterative BEMs. Thus, iterative refinement of the BIBEE/P results recovers the BEM solution; excellent agreement is obtained in only a few iterations. The boundary-integral-equation framework may also provide a means to derive rigorous results explaining how the empirical correction terms in many modern GB models significantly improve accuracy despite their simple analytical forms.

Accurate density functional calculations of near-edge x-ray and optical absorption spectra of liquid water using nonperiodic boundary conditions: the role of self-interaction and long-range effects

An integrated theoretical approach based on a discrete-continuum solvent model combined with molecular dynamics and time-dependent density functional theory leads to near-edge x-ray and optical absorption spectra of liquid water in excellent agreement with experiments. For quantitative results both self-interaction and long-range effects should be properly corrected. While scaling with respect to a monomer effectively compensates the self-interaction error, long-range effects require inclusion of exact exchange. Further, new physical insights on the not fully elucidated nature of the first optical band are provided.

A quantum mechanical/molecular dynamics/mean field study of acrolein in aqueous solution: analysis of H bonding and bulk effects on spectroscopic properties

A novel molecular dynamics methodology recently proposed by our group [Rega et al., Chem. Phys. Lett. 422, 367 (2006)], which is based on an integrated hybrid potential rooted in high level quantum mechanical methods using localized basis functions and nonperiodic boundary conditions, has been applied to study acrolein in aqueous solution. The solute structural rearrangement and its hydrogen-bonding pattern due to the interactions with water have been analyzed in some detail. Moreover, the solvent effects on the UV n-->pi* vertical transition and on the NMR 13C and 17O shielding constants of acrolein have been investigated theoretically by performing a posteriori quantum mechanical calculations on a statistically significant number of snapshots extracted from both gas-phase and aqueous solution simulations. Results show that such effective computational strategy can be successfully used to improve our understanding, at atomic level, of important spectroscopic observables.

Interplay of Intrinsic, Environmental, and Dynamic Effects in Tuning the EPR Parameters of Nitroxides: Further Insights from an Integrated Computational Approach

The role of stereoelectronic, environmental, and short-time dynamic effects in tuning the hyperfine and gyromagnetic tensors of a prototypical nitroxide spin probe has been investigated by an integrated computational approach based on extended Lagrangian molecular dynamics and discrete-continuum solvent models. Trajectories were generated in two protic solvents as well as in the gas phase for reference; structural analysis of the dynamics, and comparison with optimized solute-solvent clusters, allowed for the identification of the prevailing solute-solvent hydrogen-bonding patterns and helped to define the strategy for the computation of magnetic parameters. This was performed in a separate step, on a large number of frames, by a high-level DFT approach coupling the PBE0 hybrid functional with a tailored basis set and with proper account of specific and bulk solvent effects. Remarkable changes in solvation networks are found on going from aqueous to methanol solution, thus providing a rationalization of indirect experimentally available evidence. The computed magnetic parameters are in satisfactory agreement with the available measured values and allow for an unbiased evaluation of the role of different effects in tuning the overall EPR observables. Apart from their intrinsic interest, our results pave the route toward the development of tunable detection protocols based on specific spectroscopic signatures.

Effective potentials for 1:1 electrolyte solutions incorporating dielectric saturation and repulsive hydration

Implicit water potentials are developed for the study of thermodynamic and structural properties of solutions of NaCl, LiCl, and KCl. The interaction potential between cations and anions is parametrized from the ionic crystal potential. Two short-range corrections were added to the system to account for the water solvent. The first is due to dielectric saturation which reduces the dielectric permittivity in the vicinity of an ion. The second is a repulsive Gaussian potential which represents the first hydration shell around the ions. Grand canonical Monte Carlo simulations were performed to calculate the mean ionic activity coefficients. Molecular dynamics simulations were performed to calculate the radial distribution functions of 1.0 molal solutions at 298 K which were used to compare the structure of the explicit and implicit water simulations. The implementation of dielectric saturation and a repulsive hydration potential results in an excellent description of the mean activity coefficient and is able to capture structural features of contact ion pairs and solvent separated ions.

Change in a protein's electronic structure induced by an explicit solvent: Anab initio fragment molecular orbital study of ubiquitin

The effect of solvation on the electronic structure of the ubiquitin protein was analyzed using the ab initio fragment molecular orbital (FMO) method. FMO calculations were performed for the protein in vacuo, and the protein was immersed in an explicit solvent shell as thick as 12 A at the HF or MP2 level by using the 6-31G* basis set. The protein's physical properties examined were the net charge, the dipole moment, the internal energy, and the solvent interaction energy. Comparison of the computational results revealed the following changes in the protein upon solvation. First, the positively charged amino acid residues on the protein surface drew electrons from the solvent, while the negatively charged ones transfer electrons to the solvent. Second, the dipole moment of the protein was enhanced as a result of the polarization. Third, the internal energy of the protein was destabilized, but the destabilization was more than compensated for by the generation of a favorable protein-solvent interaction. Finally, the energetic changes were elicited both by the electron correlation effect of the first solvent shell and by the electrostatic effect of more distant solvent molecules. These findings were consistent with the picture of the solvated protein being a polarizable molecule dissolved in a dielectric media.