An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

Computing the (1)H NMR spectrum of a bulk ionic liquid from snapshots of car-parrinello molecular dynamics simulations.

We have investigated the performance of several computational protocols in predicting the NMR spectrum of a molecular ion in a complex liquid phase such as an ionic liquid. To do this, we computed the proton NMR chemical shifts of the 1-ethyl-3-methylimidazolium cation [emim](+) in [emim][Cl]. Environmental effects on the imidazolium ring proton chemical shifts are quite significant and must be taken into account explicitly. Calculations performed on the isolated imidazolium cation as well as on the [emim][Cl] ion pair grossly fail to reproduce the correct spacing between proton signals. In contrast, calculations performed on clusters extracted from the trajectory of a Car-Parrinello molecular dynamics simulation yield very good results.     

Density functional theory for efficient ab initio molecular dynamics simulations in solution

We present a density functional for first-principles molecular dynamics simulations that includes the electrostatic effects of a continuous dielectric medium. It allows for numerical simulations of molecules in solution in a model polar solvent. We propose a smooth dielectric model function to model solvation into water and demonstrate its good numerical properties for total energy calculations and constant energy molecular dynamics.     

A systematic study of minima in alanine dipeptide

The alanine dipeptide is a standard system to model dihedral angles in proteins. It is shown that obtaining the Ramachandran plot accurately is a hard problem because of many local minima; depending on the details of geometry optimizations, different Ramachandran plots can be obtained. To locate all energy minima, starting from geometries from MD simulations, 250,000 geometry optimizations were performed at the level of RHF/6-31G*, followed by re-optimizations of the located 827 minima at the level of MP2/6–311++G**, yielding 30 unique minima, most of which were not previously reported in literature. Both in vacuo and solvated structures are discussed. The minima are systematically categorized based on four backbone dihedral angles. The Gibbs energies are evaluated and the structural factors determining the relative stabilities of conformers are discussed. © 2018 Wiley Periodicals, Inc.     

Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide

Methodology is discussed for mixed ab initio quantum mechanics/molecular mechanics modeling of systems where the quantum mechanics (QM) and molecular mechanics (MM) regions are within the same molecule. The ab initio QM calculations are at the restricted Hartree–Fock level using the pseudospectral method of the Jaguar program while the MM part is treated with the OPLS force fields implemented in the IMPACT program. The interface between the QM and MM regions, in particular, is elaborated upon, as it is dealt with by “breaking” bonds at the boundaries and using Boys-localized orbitals found from model molecules in place of the bonds. These orbitals are kept frozen during QM calculations. Results from tests of the method to find relative conformational energies and geometries of alanine dipeptides and alanine tetrapeptides are presented along with comparisons to pure QM and pure MM calculations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1468–1494, 1999     

The influence of hydrophobic amino acid side groups on the acidity of the aromatic imidazole ring of histidine: A theoretical study

In this study we have calculated the acidity constant (pKa) of imidazole ring in Histidine‐Hydrophobic amino acid dipeptides using the quantum chemistry and continuum solvation methods. Density functional theory calculations with the large basis sets are used to determine the Gibbs free energy of deprotonate in the gas and liquid phases. Based on our results ΔG_S values are located between ?69.38 and ?18.82 kcal mol^?1 which are related to His⁺–Gly and His forms, respectively. pKa of the dipeptides in the aqueous phase was obtained from the calculated gas‐phase and solvation free energies through a thermodynamic cycle and the solvation model chemistry of Martin Karplus et al. Solvation effects are treated using a self‐consistent reaction field formalism involving polarized continuum models. According to our calculations pKa values are between 5.50 and 8.19 that are belong to His⁺–ILe and His⁺–Ala forms, respectively. Natural bond orbital analysis of dipeptides reveals that the electron delocalization in imidazole ring is the most effective factor in determination of acidity order for these compounds. Structural analysis confirmed that the orientation of carbonyl group with respect to imidazole ring is an effective factor in imidazole ring stability. © 2011 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Ionization and protonation of MgC3: A theoretical study

A theoretical study of the MgC₃⁺ and MgC₃H⁺ species has been carried out. Predictions for their geometries and vibrational frequencies have been made at both second‐order Møller–Plesset (MP2) and B3LYP levels, whereas electronic energies have been computed at G2 and coupled cluster single and double excitation model augmented with a noniterative triple excitation conection (CCSD(T)) levels. The predicted global minimum for MgC₃⁺ is a rhombic structure (²A₁ electronic state), whereas a T‐shaped structure and an open‐chain isomer lie about 10 and 12 kcal/mol, respectively, higher in energy. In the case of MgC₃H⁺ the predicted global minimum is also a four‐membered ring obtained upon protonation of the most stable neutral isomer. Low ionization potentials and high proton affinities are generally obtained, especially for the most stable MgC₃⁺ isomer. The estimated values at the CCSD(T) level for the predicted global minimum are 7.20 eV [ionization potential (IP)] and 256.5 kcal/mol [proton affinities (PA)]. Therefore, if present in the interstellar medium, MgC₃ should be easily ionized and would react quite easily to give the protonated species. © 2002 John Wiley & Sons, Inc. Int J Quantum Chem, 2001     

Hydrogen bonding, solvation, and hydrolysis of cisplatin: a theoretical study

Density functional calculations on a range of hydrogen bonded clusters of cisplatin are reported. A systematic search of 1:1 cisplatin:water complexes reveals only three stable minima, which contain a number of common, recurring interactions, such as an N-H...O-H...Cl bridging mode. Expanding these clusters by adding water molecules leads to a model of the first solvation shell of cisplatin, which contains the above motifs along with several strong water-water interactions. The strengths of such interactions are rationalized on the basis of electrostatic potentials, and quantified by use of Atoms in Molecules properties. This analysis also allows us to estimate cisplatin's position on Abraham's hydrogen bond acidity and basicity scales, indicating that cisplatin is a strong donor and acceptor of hydrogen bonds due to the dominance of hard, electrostatic interactions. The effects of this explicit solvation on the barrier to hydrolysis, and hence activation, of cisplatin are explored, indicating a slightly higher barrier than in the gas phase, leading to better agreement with experiment than either gas phase or continuum solvation calculations.     

Protonation Equilibrium in the Active Site of the Photoactive Yellow Protein

The role and existence of low-barrier hydrogen bonds (LBHBs) in enzymatic and protein activity has been largely debated. An interesting case is that of the photoactive yellow protein (PYP). In this protein, two short HBs adjacent to the chromophore, p-coumaric acid (pCA), have been identified by X-ray and neutron diffraction experiments. However, there is a lack of agreement on the chemical nature of these H-bond interactions. Additionally, no consensus has been reached on the presence of LBHBs in the active site of the protein, despite various experimental and theoretical studies having been carried out to investigate this issue. In this work, we perform a computational study that combines classical and density functional theory (DFT)-based quantum mechanical/molecular mechanical (QM/MM) simulations to shed light onto this controversy. Furthermore, we aim to deepen our understanding of the chemical nature and dynamics of the protons involved in the two short hydrogen bonds that, in the dark state of PYP, connect pCA with the two binding pocket residues (E46 and Y42). Our results support the existence of a strong LBHB between pCA and E46, with the H fully delocalized and shared between both the carboxylic oxygen of E46 and the phenolic oxygen of pCA. Additionally, our findings suggest that the pCA interaction with Y42 can be suitably described as a typical short ionic H-bond of moderate strength that is fully localized on the phenolic oxygen of Y42.     

A comparative study of O2, CO,and NO binding to iron–porphyrin

Minimum-energy structures of O₂, CO, and NO iron–porphyrin (FeP) complexes, computed with the Car–Parrinello molecular dynamics, agree well with the available experimental data for synthetic heme models. The diatomic molecule induces a 0.3–0.4 Å displacement of the Fe atom out of the porphyrin nitrogen (N_p) plane and a doming of the overall porphyrin ring. The energy of the iron–diatomic bond increases in the order Fe(SINGLE BOND)O₂ (9 kcal/mol) < Fe(SINGLE BOND)CO (26 kcal/mol) < Fe(SINGLE BOND)NO (35 kcal/mol). The presence of an imidazole axial ligand increases the strength of the Fe(SINGLE BOND)O₂ and Fe(SINGLE BOND)CO bonds (15 and 35 kcal/mol, respectively), with few structural changes with respect to the FeP(CO) and FeP(O₂) complexes. In contrast, the imidazole ligand does not affect the energy of the Fe(SINGLE BOND)NO bond, but induces significant structural changes with respect to the FeP(NO) complex. Similar variations in the iron–imidazole bond with respect to the addition of CO, O₂, and NO are also discussed. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 31–35, 1998     

Ionic dissociations of chlorosulfonic acid in microsolvated clusters: A density functional theory and ab initio MO study

Ionic dissociation of chlorosulfonic acid (HSO3Cl) in the molecular clusters HSO3Cl-(H2O)n (n = 1-4) and HSO3Cl-NH3-(H2O)n (n = 0-3) was investigated by density functional theory and ab initio molecular orbital theory. The equilibrium structures, binding energies, and thermodynamic properties, such as relative enthalpy and relative Gibbs free energy, and were calculated using the hybrid density func- tional (B3LYP) method and the second order M?ller-Plesset approximation (MP2) method with the 6-311 G** basis set. Chlorosulfonic acid was found to require a minimum of three water molecules for ionization to occur and at least one water molecule to protonate ammonia. The corresponding clusters with fewer water molecules were found to be strongly hydrogen-bonded. The related properties and acid strength of chlorosulfonic acid were discussed and compared to the acid strengths of perchloric acid and sulfuric acid in the context of clusters with ammonia and water. The relative stabilities of these clusters were also investigated.     

Effect of solvation on the ionization of guanine nucleotide: A hybrid QM/EFP study

Ionization of nucleobases is affected by their biological environment, which includes both the effect of adjacent nucleotides as well as the presence of water around it. Guanine and its nucleotide have the lowest ionization potentials among the various DNA bases. Therefore, the threshold of ionization is dependent on that of guanine and its characterization is crucial to the prediction of interaction of light with DNA. We investigate the effect of solvation on the vertical ionization energies (VIEs) of guanine and its nucleotide. In this work, we have used hybrid quantum mechanics/molecular mechanics (QM/MM) approach with effective fragment potential as the MM method of choice and equation‐of‐motion coupled‐cluster for ionization potential with singles and doubles (EOM‐IP‐CCSD) as the QM method. The performance of the hybrid scheme with respect to the full QM method shows an accuracy of 0.02–0.04 eV. The lowest few ionizations of the nucleotide are found to be from different parts of the moiety, that is, the nucleic acid base, phosphate, or sugar, and these ionization energies are very closely spaced giving rise to a very complicated spectrum. Furthermore, microsolvation has large effects on these ionizations and can lead to red or blue shift depending on the position of the water molecule. Even a single water molecule can change the order of ionized states in the nucleotide. The VIEs of the bulk solvated chromophores are predicted and compared to existing experimental spectra. The predominant role of polarization in the solvatochromic shift is noticed. © 2017 Wiley Periodicals, Inc.     

Determining the strengths of hydrogen bonds in solid-state ammonia and urea: insight from periodic DFT calculations

Plane-wave density functional theory has been applied to determine the strengths of hydrogen bonds in the phase I crystal structures of ammonia and urea. For ammonia, each component of the trifurcated hydrogen bond has been found to be almost as strong as a standard N-H.N interaction, and for urea the strengths of the two different N-H.O interactions have been determined by a quantum mechanical technique for the first time.     

van der Waals interactions are critical in Car–Parrinello molecular dynamics simulations of porphyrin–fullerene dyads

The interplay between electrostatic and van der Waals (vdW) interactions in porphyrin‐C₆₀ dyads is still under debate despite its importance in influencing the structural characteristics of such complexes considered for various applications in molecular photovoltaics. In this article, we sample the conformational space of a porphyrin‐C₆₀ dyad using Car–Parrinello molecular dynamics simulations with and without empirical vdW corrections. Long‐range vdW interactions, which are poorly described by the commonly used density functional theory functionals, prove to be essential for a proper dynamics of the dyad moieties. Inclusion of vdW corrections brings porphyrin and C₆₀ close together in an orientation that is in agreement with experimental observations. The structural differences arising from the vdW corrections are shown to be significant for several properties and potentially less important for others. Additionally, our Mulliken population analysis reveals that contrary to the common belief, porphyrin is not the primary electron donating moiety for C₆₀. In the considered dyad, fullerene's affinity for electrons is primarily satisfied by charge transfer from the amide group of the linker. However, we show that in the absence of another suitable bound donor, C₆₀ can withdraw electrons from porphyrin if it is sufficiently close. © 2015 Wiley Periodicals, Inc.     

Conformational preference and mechanism of decarboxylation of levodopa. A quantum dynamics/quantum mechanics study

The present study addresses the conformational preferences and the mechanism of decarboxylation of levodopa (LD). LD is used to increase dopamine concentrations in the treatment of Parkinson's disease. LD crosses the protective blood–brain barrier, where it is converted into dopamine by the process of decarboxylation. Molecular dynamics simulation has been carried out at the DFT/6‐31++G level of theory to identify the global minimum structure of LD. Conformational preferences of the amino acid side chain of LD has been investigated at the B3LYP/6‐311++G** level of theory. Fourier transform analysis has been performed to identify the origin of the rotational barriers. Electrostatic dipole moment and bond interactions underlie the observed potential energy barriers for rotation of the amino acid side chain of LD. The vital biological process of decarboxylation of LD has been examined in the gas phase and in aqueous solution. Without the presence of water, there is only one possible route for the decarboxylation of LD. In this concerted mechanism, a proton transfer and breakage of the C₁₀—C₁₈ bond, take place simultaneously (ΔE^# = 73.2 kcal/mol). In solution, however, two possible decarboxylation routes are available for LD. The first involve the formation of a zwitterionic intermediate (ΔE^# = 72.4 kcal/mol). The zwitterionic form of LD have been localized using explicitly bound water molecules to model short‐range solvent effects and self‐consistent reaction field polarized continuum model to estimate long‐range solvent interactions. The second route involve the formation of a cyclic structure in which a water molecule acts as a bridge linking the anticarboxylic hydrogen and α‐position carbon atom (ΔE^# = 59.8 kcal/mol). Natural bond orbital (NBO) analysis reveals that the conformational and overall stability of the amino acid side chain is facilitated by the antiperiplanar interactions between the phenyl moiety C—H and C—C bonds and C—X bonds of the amino acid side chain. However, much of the major donor–acceptor interactions is of the lone pair type and is localized within the amino acid side chain itself. Results of the present work reveal that NBO data reflect nicely and identify clearly reaction coordinates at the transition species. © 2013 Wiley Periodicals, Inc.