Free energy perturbation study of water dimer dissociation kinetics

An efficient approach is described for using accurate ab initio calculations to determine the rates of elementary condensation and evaporation processes that lead to nucleation of aqueous aerosols. The feasibility of the method is demonstrated in an application to evaporation rates of water dimer at 230 K. The method, known as ABC-FEP (ab initio/classical free energy perturbation), begins with a calculation of the potential of mean force for the dissociation (evaporation) of small water clusters using a molecular dynamics (MD) simulation with a model potential. The free energy perturbation is used to calculate how changing from the model potential to a potential calculated from ab initio methods would alter the potential of mean force. The difference in free energy is the Boltzmann-weighted average of the difference between the ab initio and classical potential energies, with the average taken over a sample of configurations from the MD simulation. In principle, the method does not require a highly accurate model potential, though more accurate potentials require fewer configurations to achieve a small sampling error in the free energy perturbation step. To test the feasibility of obtaining accurate potentials of mean force from ab initio calculations at a modest number of configurations, the free energy perturbation method has been used to correct the errors when some standard models for bulk water (SPC, TIP4P, and TIP4PFQ) are applied to water dimer. To allow a thorough exploration of sampling issues, a highly accurate fit to results of accurate ab initio calculations, known as SAPT-5s, as been used a proxy for the ab initio calculations. It is shown that accurate values for a point on the potential of mean force can be obtained from any of the water models using ab initio calculations at only 50 configurations. Thus, this method allows accurate simulations of small clusters without the need to develop water models specifically for clusters.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.     

Peptide free‐energy profile is strongly dependent on the force field: Comparison of C96 and AMBER95

The C96 and AMBER95 force fields were compared with small model peptides Ac‐(Ala)_n‐NMe (Ac = CH₃CO, NMe = NHCH₃, n=2 and 3) in vacuo and in TIP3P water by computing the free‐energy profiles using multicanonical molecular dynamics method. The C96 force field is a modified version of the AMBER95 force field, which was adjusted to reproduce the energy difference between extended β‐ and constrained α‐helical energies for the alanine tetrapeptide, obtained by the high level ab initio MO method. The slight modification resulted in a large difference in the free energy profiles. The C96 force field prefers relatively extended conformers, whereas the AMBER95 force field favors turn conformations. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 748–762, 2000     

Temperature and pressure dependence of alanine dipeptide studied by multibaric-multithermal molecular dynamics simulations

We applied the multibaric-multithermal (MUBATH) molecular dynamics (MD) algorithm to an alanine dipeptide in explicit water. The MUBATH MD simulation covered a wide range of conformational space and sampled the states of PII, C5, alphaR, alphaP, alphaL, and C7(ax). On the other hand, the conventional isobaric-isothermal simulation was trapped in local-minimum free-energy states and sampled only a few of them. We calculated the partial molar enthalpy difference DeltaH and partial molar volume difference DeltaV among these states by the MUBATH simulation using the AMBER parm99 and AMBER parm96 force fields and two sets of initial conditions. We compared these results with those from Raman spectroscopy experiments. The Raman spectroscopy data of DeltaH for the C5 state against the PII state agreed with both MUBATH data with the AMBER parm96 and parm99 force fields. The partial molar enthalpy difference DeltaH for the alphaR state and the partial molar volume difference DeltaV for the C5 state by the Raman spectroscopy agreed with those for the AMBER parm96 force field. On the other hand, DeltaV for the alphaR state by the Raman spectroscopy was consistent with our AMBER-parm99 force-field result. All the experimental results fall between those of simulations using AMBER parm96 and parm99 force fields, suggesting that the ideal force-field parameters lie between those of AMBER parm96 and parm99.     

Controlling the shape and flexibility of arylamides: a combined ab initio, ab initio molecular dynamics, and classical molecular dynamics study

Using quantum chemistry plus ab initio molecular dynamics and classical molecular dynamics methods, we address the relationship between molecular conformation and the biomedical function of arylamide polymers. Specifically, we have developed new torsional parameters for a class of these polymers and applied them in a study of the interaction between a representative arylamide and one of its biomedical targets, the anticoagulant drug heparin. Our main finding is that the torsional barrier of a C(aromatic)-C(carbonyl) bond increases significantly upon addition of an o-OCH2CH2NH3+ substituent on the benzene ring. Our molecular dynamics studies that are based on the original general AMBER force field (GAFF) and GAFF modified to include our newly developed torsional parameters show that the binding mechanism between the arylamide and heparin is very sensitive to the choice of torsional potentials. Ab initio molecular dynamics simulation of the arylamide independently confirms the degree of flexibility we obtain by classical molecular dynamics when newly developed torsional potentials are used.     

Ab initio calculations of the melting temperatures of refractory bcc metals

We present ab initio calculations of the melting temperatures for bcc metals Nb, Ta and W. The calculations combine phase coexistence molecular dynamics (MD) simulations using classical embedded-atom method potentials and ab initio density functional theory free energy corrections. The calculated melting temperatures for Nb, Ta and W are, respectively, within 3%, 4%, and 7% of the experimental values. We compare the melting temperatures to those obtained from direct ab initio molecular dynamics simulations and see if they are in excellent agreement with each other. The small remaining discrepancies with experiment are thus likely due to inherent limitations associated with exchange-correlation energy approximations within density-functional theory.     

An ab initio force field for the cofactors of bacterial photosynthesis

This article presents a new ab initio force field for the cofactors of bacterial photosynthesis, namely quinones and bacteriochlorophylls. The parameters has been designed to be suitable for molecular dynamics simulations of photosynthetic proteins by being compatible with the AMBER force field. To our knowledge, this is the first force field for photosynthetic cofactors based on a reliable set of ab initio density functional reference data for methyl bacteriochlorophyll a, methyl bacteriopheophytin a, and of a derivative of ubiquinone. Indeed, the new molecular mechanics force field is able to reproduce very well not only the experimental and ab initio structural properties and the vibrational spectra of the molecules, but also the eigenvectors of the molecular normal modes. For this reason it might also be helpful to understand vibrational spectroscopy results obtained on reaction center proteins.     

Potential-energy and free-energy surfaces of glycyl-phenylalanyl-alanine (GFA) tripeptide: experiment and theory

The free-energy surface (FES) of glycyl-phenylalanyl-alanine (GFA) tripeptide was explored by molecular dynamics (MD) simulations in combination with high-level correlated ab initio quantum chemical calculations and metadynamics. Both the MD and metadynamics employed the tight-binding DFT-D method instead of the AMBER force field, which yielded inaccurate results. We classified the minima localised in the FESs as follows: a) the backbone-conformational arrangement; and b) the existence of a COOH...OC intramolecular H-bond (families CO(2)H(free) and CO(2)H(bonded)). Comparison with experimental results showed that the most stable minima in the FES correspond to the experimentally observed structures. Remarkably, however, we did not observe experimentally the CO(2)H(bonded) family (also predicted by metadynamics), although its stability is comparable to that of the CO(2)H(free) structures. This fact was explained by the former's short excited-state lifetime. We also carried out ab initio calculations using DFT-D and the M06-2X functional. The importance of the dispersion energy in stabilising peptide conformers is well reflected by our pioneer analysis using the DFT-SAPT method to explore the nature of the backbone/side-chain interactions.     

Derivation of class II force fields. VIII. Derivation of a general quantum mechanical force field for organic compounds

A class II valence force field covering a broad range of organic molecules has been derived employing ab initio quantum mechanical "observables." The procedure includes selecting representative molecules and molecular structures, and systematically sampling their energy surfaces as described by energies and energy first and second derivatives with respect to molecular deformations. In this article the procedure for fitting the force field parameters to these energies and energy derivatives is briefly reviewed. The application of the methodology to the derivation of a class II quantum mechanical force field (QMFF) for 32 organic functional groups is then described. A training set of 400 molecules spanning the 32 functional groups was used to parameterize the force field. The molecular families comprising the functional groups and, within each family, the torsional angles used to sample different conformers, are described. The number of stationary points (equilibria and transition states) for these molecules is given for each functional group. This set contains 1324 stationary structures, with 718 minimum energy structures and 606 transition states. The quality of the fit to the quantum data is gauged based on the deviations between the ab initio and force field energies and energy derivatives. The accuracy with which the QMFF reproduces the ab initio molecular bond lengths, bond angles, torsional angles, vibrational frequencies, and conformational energies is then given for each functional group. Consistently good accuracy is found for these computed properties for the various types of molecules. This demonstrates that the methodology is broadly applicable for the derivation of force field parameters across widely differing types of molecular structures. Copyright 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1782-1800, 2001     

Direct assessment of quantum nuclear effects on hydrogen bond strength by constrained-centroid ab initio path integral molecular dynamics

The impact of quantum nuclear effects on hydrogen (H-) bond strength has been inferred in earlier work from bond lengths obtained from path integral molecular dynamics (PIMD) simulations. To obtain a direct quantitative assessment of such effects, we use constrained-centroid PIMD simulations to calculate the free energy changes upon breaking the H-bonds in dimers of HF and water. Comparing ab initio simulations performed using PIMD and classical nucleus molecular dynamics (MD), we find smaller dissociation free energies with the PIMD method. Specifically, at 50 K, the H-bond in (HF)(2) is about 30% weaker when quantum nuclear effects are included, while that in (H(2)O)(2) is about 15% weaker. In a complementary set of simulations, we compare unconstrained PIMD and classical nucleus MD simulations to assess the influence of quantum nuclei on the structures of these systems. We find increased heavy atom distances, indicating weakening of the H-bond consistent with that observed by direct calculation of the free energies of dissociation.     

Effects of all-atom force fields on amyloid oligomerization: replica exchange molecular dynamics simulations of the Aβ(16-22) dimer and trimer

The aim of this work is to investigate the effects of molecular mechanics force fields on amyloid peptide assembly. To this end, we performed extensive replica exchange molecular dynamics (REMD) simulations on the monomer, dimer and trimer of the seven-residue fragment of the Alzheimer's amyloid-β peptide, Aβ(16-22), using the AMBER99, GROMOS96 and OPLS force fields. We compared the force fields by analysing the resulting global and local structures as well as the free energy landscapes at 300 K. We show that AMBER99 strongly favors helical structures for the monomer and does not predict any β-sheet structure for the dimer and trimer. In contrast, the dimer and trimer modeled by GROMOS96 form antiparallel β-sheet structures, while OPLS predicts diverse structures. Overall, the free energy landscapes obtained by three force fields are very different, and we also note a weak structural dependence of our results on temperature. The implications of this computational study on amyloid oligomerization, fibril growth and inhibition are also discussed.     

Overcoming entropic barrier with coupled sampling at dual resolutions

An enhanced sampling method is proposed for ab initio protein folding simulations. The new method couples a high-resolution model for accuracy and a low-resolution model for efficiency. It aims to overcome the entropic barrier found in the exponentially large protein conformational space when a high-resolution model, such as an all-atom molecular mechanics force field, is used. The proposed method is designed to satisfy the detailed balance condition so that the Boltzmann distribution can be generated in all sampling trajectories in both high and low resolutions. The method was tested on model analytical energy functions and ab initio folding simulations of a beta-hairpin peptide. It was found to be more efficient than replica-exchange method that is used as its building block. Analysis with the analytical energy functions shows that the number of energy calculations required to find global minima and to converge mean potential energies is much fewer with the new method. Ergodic measure shows that the new method explores the conformational space more rapidly. We also studied imperfect low-resolution energy models and found that the introduction of errors in low-resolution models does decrease its sampling efficiency. However, a reasonable increase in efficiency is still observed when the global minima of the low-resolution models are in the vicinity of the global minimum basin of the high-resolution model. Finally, our ab initio folding simulation of the tested peptide shows that the new method is able to fold the peptide in a very short simulation time. The structural distribution generated by the new method at the equilibrium portion of the trajectory resembles that in the equilibrium simulation starting from the crystal structure.     

Frozen density functional free energy simulations of redox proteins: computational studies of the reduction potential of plastocyanin and rusticyanin

The evaluation of reduction potentials of proteins by ab initio approaches presents a major challenge for computational chemistry. This is addressed in the present investigation by reporting detailed calculations of the reduction potentials of the blue copper proteins plastocyanin and rusticyanin using the QM/MM all-atom frozen density functional theory, FDFT, method. The relevant ab initio free energies are evaluated by using a classical reference potential. This approach appears to provide a general consistent and effective way for reproducing the configurational ensemble needed for consistent ab initio free energy calculations. The FDFT formulation allows us to treat a large part of the protein quantum mechanically by a consistently coupled QM/QM/MM embedding method while still retaining a proper configurational sampling. To establish the importance of proper configurational sampling and the need for a complete representation of the protein+solvent environment, we also consider several classical approaches. These include the semi-macroscopic PDLD/S-LRA method and classical all-atom simulations with and without a polarizable force field. The difference between the reduction potentials of the two blue copper proteins is reproduced in a reasonable way, and its origin is deduced from the different calculations. It is found that the protein permanent dipole tunes down the reduction potential for plastocyanin compared to the active site in regular water solvent, whereas in rusticyanin it is instead tuned up. This electrostatic environment, which is the major effect determining the reduction potential, is a property of the entire protein and solvent system and cannot be ascribed to any particular single interaction.     

Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes

A hybrid quantum/classical path integral Monte Carlo (QC-PIMC) method for calculating the quantum free energy barrier for hydrogen transfer reactions in condensed phases is presented. In this approach, the classical potential of mean force along a collective reaction coordinate is calculated using umbrella sampling techniques in conjunction with molecular dynamics trajectories propagated according to a mapping potential. The quantum contribution is determined for each configuration along the classical trajectory with path integral Monte Carlo calculations in which the beads move according to an effective mapping potential. This type of path integral calculation does not utilize the centroid constraint and can lead to more efficient sampling of the relevant region of conformational space than free-particle path integral sampling. The QC-PIMC method is computationally practical for large systems because the path integral sampling for the quantum nuclei is performed separately from the classical molecular dynamics sampling of the entire system. The utility of the QC-PIMC method is illustrated by an application to hydride transfer in the enzyme dihydrofolate reductase. A comparison of this method to the quantized classical path and grid-based methods for this system is presented.     

Ab initio molecular dynamics study of formate ion hydration

We perform ab initio molecular dynamics simulations of the aqueous formate ion. The mean number of water molecules in the first solvation shell, or the hydration number, of each formate oxygen is found to be consistent with recent experiments. Our ab initio pair correlation functions, however, differ significantly from many classical force field results and hybrid quantum mechanics/molecular mechanics predictions. They yield roughly one less hydrogen bond between each formate oxygen and water than force field or hybrid methods predict. Both the BLYP and PW91 exchange correlation functionals give qualitatively similar results. The time dependence of the hydration numbers are examined, and Wannier function techniques are used to analyze electronic configurations along the molecular dynamics trajectory.     

Insights from first principles molecular dynamics studies toward infrared multiple-photon and single-photon action spectroscopy: case study of the proton-bound dimethyl ether dimer

We have used finite temperature ab initio molecular dynamics simulations in conjunction with computation of critical quantum nuclear effects to probe the differences between single-photon argon tagged action spectral results and infrared multiple-photon dissociation experiments for a proton bound molecular ion system. We find that the principal difference between the results in these experimental techniques is essentially that of cluster temperature. The multiple-photon dissociation experiments conducted using room temperature ions reflect a larger degree of conformational freedom compared to the colder single-photon argon tagged action spectral results. Our ab initio molecular dynamics simulation techniques accurately capture the effects of conformational sampling, adequately reproduce both spectra, and can be utilized to assign the dynamically averaged finite temperature spectra.