Showing papers in "Journal of Chemical Theory and Computation in 2007"
TL;DR: There is no one perfect "one size fits all" algorithm for clustering MD trajectories and that the results strongly depend on the choice of atoms for the pairwise comparison, so the best performance was observed with the average-linkage, means, and SOM algorithms.
Abstract: Molecular dynamics simulation methods produce trajectories of atomic positions (and optionally velocities and energies) as a function of time and provide a representation of the sampling of a given molecule's energetically accessible conformational ensemble. As simulations on the 10-100 ns time scale become routine, with sampled configurations stored on the picosecond time scale, such trajectories contain large amounts of data. Data-mining techniques, like clustering, provide one means to group and make sense of the information in the trajectory. In this work, several clustering algorithms were implemented, compared, and utilized to understand MD trajectory data. The development of the algorithms into a freely available C code library, and their application to a simple test example of random (or systematically placed) points in a 2D plane (where the pairwise metric is the distance between points) provide a means to understand the relative performance. Eleven different clustering algorithms were developed, ranging from top-down splitting (hierarchical) and bottom-up aggregating (including single-linkage edge joining, centroid-linkage, average-linkage, complete-linkage, centripetal, and centripetal- complete) to various refinement (means, Bayesian, and self-organizing maps) and tree (COBWEB) algorithms. Systematic testing in the context of MD simulation of various DNA systems (including DNA single strands and the interaction of a minor groove binding drug DB226 with a DNA hairpin) allows a more direct assessment of the relative merits of the distinct clustering algorithms. Additionally, means to assess the relative performance and differences between the algorithms, to dynamically select the initial cluster count, and to achieve faster data mining by "sieved clustering" were evaluated. Overall, it was found that there is no one perfect "one size fits all" algorithm for clustering MD trajectories and that the results strongly depend on the choice of atoms for the pairwise comparison. Some algorithms tend to produce homogeneously sized clusters, whereas others have a tendency to produce singleton clusters. Issues related to the choice of a pairwise metric, clustering metrics, which atom selection is used for the comparison, and about the relative performance are discussed. Overall, the best performance was observed with the average-linkage, means, and SOM algorithms. If the cluster count is not known in advance, the hierarchical or average-linkage clustering algorithms are recommended. Although these algorithms perform well, it is important to be aware of the limitations or weaknesses of each algorithm, specifically the high sensitivity to outliers with hierarchical, the tendency to generate homogenously sized clusters with means, and the tendency to produce small or singleton clusters with average-linkage.
716 citations
TL;DR: The comparison shows that two newly developed density functional theory methods, PWB6K and M05-2X, give the best performance for this benchmark database of 22 noncovalent complexes, including both hydrogen-bonding and dispersion-dominated complexes.
Abstract: Forty density functionals and one wavefunction method are assessed against a recently published database of accurate noncovalent interaction energies of biological importance. The comparison shows that two newly developed density functional theory (DFT) methods, PWB6K and M05-2X, give the best performance for this benchmark database of 22 noncovalent complexes, including both hydrogen-bonding and dispersion-dominated complexes. In contrast, the more popular B3LYP and PBEh functionals fail to describe the interactions in the dispersion-dominated complexes. The local spin density approximation and BHandH functionals give good performance for dispersion-dominated interactions at the expense of a large error for hydrogen bonding. PWB6K and M05-2X constitute a new generation of DFT methods based on simultaneously optimized exchange and correlation functionals that include kinetic energy density in both the exchange and correlation functional, and the present study confirms that they have greatly improved performance for noncovalent interactions as compared to previous DFT methods. We interpret this as being due to an improved treatment of medium-range correlation effects by the exchange-correlation functional. We recommend the PWB6K and M05-2X methods for investigating large biological systems and soft materials.
568 citations
TL;DR: Variants of WHAM, STWHAM and PTWHAM are presented, derived with the same set of assumptions, that can be directly applied to several generalized ensemble algorithms, including simulated tempering, parallel tempering (better known as replica-exchange among temperatures), and replica- Exchanges simulated Tempering.
Abstract: The growing adoption of generalized-ensemble algorithms for biomolecular simulation has resulted in a resurgence in the use of the weighted histogram analysis method (WHAM) to make use of all data generated by these simulations. Unfortunately, the original presentation of WHAM by Kumar et al. is not directly applicable to data generated by these methods. WHAM was originally formulated to combine data from independent samplings of the canonical ensemble, whereas many generalized-ensemble algorithms sample from mixtures of canonical ensembles at different temperatures. Sorting configurations generated from a parallel tempering simulation by temperature obscures the temporal correlation in the data and results in an improper treatment of the statistical uncertainties used in constructing the estimate of the density of states. Here we present variants of WHAM, STWHAM and PTWHAM, derived with the same set of assumptions, that can be directly applied to several generalized ensemble algorithms, including simulated tempering, parallel tempering (better known as replica-exchange among temperatures), and replica-exchange simulated tempering. We present methods that explicitly capture the considerable temporal correlation in sequentially generated configurations using autocorrelation analysis. This allows estimation of the statistical uncertainty in WHAM estimates of expectations for the canonical ensemble. We test the method with a one-dimensional model system and then apply it to the estimation of potentials of mean force from parallel tempering simulations of the alanine dipeptide in both implicit and explicit solvent.
457 citations
TL;DR: A data set of 19 second-row transition-metal complexes has been collated from sufficiently precise gas-phase electron-diffraction experiments and used for evaluating errors in DFT optimized geometries and the TPSSh hybrid meta-GGA is slightly inferior to the best hybrid GGAs.
Abstract: A set of 41 metal-ligand bond distances in 25 third-row transition-metal complexes, for which precise structural data are known in the gas phase, is used to assess optimized and zero-point averaged geometries obtained from DFT computations with various exchange-correlation functionals and basis sets. For a given functional (except LSDA) Stuttgart-type quasi-relativistic effective core potentials and an all-electron scalar relativistic approach (ZORA) tend to produce very similar geometries. In contrast to the lighter congeners, LSDA affords reasonably accurate geometries of 5d-metal complexes, as it is among the functionals with the lowest mean and standard deviations from experiment. For this set the ranking of some other popular density functionals, ordered according to decreasing standard deviation, is BLYP > VSXC > BP86 ≈ BPW91 ≈ TPSS ≈ B3LYP ≈ PBE > TPSSh > B3PW91 ≈ B3P86 ≈ PBE hybrid. In this case hybrid functionals are superior to their nonhybrid variants. In addition, we have reinvestigated the previous test sets for 3d- (Buhl M.; Kabrede, H. J. Chem. Theory Comput. 2006, 2, 1282-1290) and 4d- (Waller, M. P.; Buhl, M. J. Comput. Chem. 2007, 28, 1531-1537) transition-metal complexes using all-electron scalar relativistic DFT calculations in addition to the published nonrelativistic and ECP results. For this combined test set comprising first-, second-, and third-row metal complexes, B3P86 and PBE hybrid are indicated to perform best. A remarkably consistent standard deviation of around 2 pm in metal-ligand bond distances is achieved over the entire set of d-block elements.
454 citations
TL;DR: A new universal continuum solvation model (where "universal" denotes applicable to all solvents), called SM8, is presented, and it improves on earlier SMx universal solvation models by including free energies of solvation of ions in nonaqueous media in the parametrization.
Abstract: A new universal continuum solvation model (where “universal” denotes applicable to all solvents), called SM8, is presented. It is an implicit solvation model, also called a continuum solvation model, and it improves on earlier SMx universal solvation models by including free energies of solvation of ions in nonaqueous media in the parametrization. SM8 is applicable to any charged or uncharged solute composed of H, C, N, O, F, Si, P, S, Cl, and/or Br in any solvent or liquid medium for which a few key descriptors are known, in particular dielectric constant, refractive index, bulk surface tension, and acidity and basicity parameters. It does not require the user to assign molecular-mechanics types to an atom or group; all parameters are unique and continuous functions of geometry. It may be used with any level of electronic structure theory as long as accurate partial charges can be computed for that level of theory; we recommend using it with self-consistently polarized Charge Model 4 or other self-consis...
419 citations
TL;DR: An earlier pairwiseGB model is extended by a simple analytic correction term that largely alleviates the problem by correctly describing the solvent-excluded volume of each pair of atoms and preserves the efficiency of the pairwise GB models while making them a better approximation to reality.
Abstract: Generalized Born (GB) models provide a computationally efficient means of representing the electrostatic effects of solvent and are widely used, especially in molecular dynamics (MD). A class of particularly fast GB models is based on integration over an interior volume approximated as a pairwise union of atom spheres-effectively, the interior is defined by a van der Waals rather than Lee-Richards molecular surface. The approximation is computationally efficient, but if uncorrected, allows for high dielectric (water) regions smaller than a water molecule between atoms, leading to decreased accuracy. Here, an earlier pairwise GB model is extended by a simple analytic correction term that largely alleviates the problem by correctly describing the solvent-excluded volume of each pair of atoms. The correction term introduces a free energy barrier to the separation of non-bonded atoms. This free energy barrier is seen in explicit solvent and Lee-Richards molecular surface implicit solvent calculations, but has been absent from earlier pairwise GB models. When used in MD, the correction term yields protein hydrogen bond length distributions and polypeptide conformational ensembles that are in better agreement with explicit solvent results than earlier pairwise models. The robustness and simplicity of the correction preserves the efficiency of the pairwise GB models while making them a better approximation to reality.
334 citations
TL;DR: An overview of the SIBFA polarizable molecular mechanics procedure, which is formulated and calibrated on the basis of quantum chemistry (QC), and the development of a novel methodology, the Gaussian electrostatic model (GEM), which relies on ab initio-derived fragment electron densities to compute the components of the total interaction energy.
Abstract: We present an overview of the SIBFA polarizable molecular mechanics procedure, which is formulated and calibrated on the basis of quantum chemistry (QC). It embodies nonclassical effects such as electrostatic penetration, exchange-polarization, and charge transfer. We address the issues of anisotropy, nonadditivity, and transferability by performing parallel QC computations on multimolecular complexes. These encompass multiply H-bonded complexes and polycoordinated complexes of divalent cations. Recent applications to the docking of inhibitors to Zn-metalloproteins are presented next, namely metallo-β-lactamase, phosphomannoisomerase, and the nucleocapsid of the HIV-1 retrovirus. Finally, toward third-generation intermolecular potentials based on density fitting, we present the development of a novel methodology, the Gaussian electrostatic model (GEM), which relies on ab initio-derived fragment electron densities to compute the components of the total interaction energy. As GEM offers the possibility of a...
328 citations
TL;DR: The present work describes the development of microscopic polarizable force fields starting with the introduction of these powerful tools and following some of the subsequent developments in the field.
Abstract: A consistent treatment of electrostatic energies is arguably the most important requirement for the realistic modeling of biological systems. An important part of electrostatic modeling is the ability to account for the polarizability of the simulated system. This can be done both macroscopically and microscopically, but the use of macroscopic models may lead to conceptual traps, which do not exist in the microscopic treatments. The present work describes the development of microscopic polarizable force fields starting with the introduction of these powerful tools and following some of the subsequent developments in the field. Special effort has been made to review a wide range of applications and emphasize cases when the use of polarizable force fields is important. Finally, a brief perspective is given on the future of this rapidly growing field.
314 citations
TL;DR: The split valence bases of the 6-31G variety provide accuracies similar to those of the more computationally expensive Dunning type basis sets, and it is observed that the hybrid-meta-GGA functionals are typically among the most accurate functionals for all of the properties examined in this work.
Abstract: The reliable prediction of molecular properties is a vital task of computational chemistry. In recent years, density functional theory (DFT) has become a popular method for calculating molecular properties for a vast array of systems varying in size from small organic molecules to large biological compounds such as proteins. In this work we assess the ability of many DFT methods to accurately determine atomic and molecular properties for small molecules containing elements commonly found in proteins, DNA, and RNA. These properties include bond lengths, bond angles, ground state vibrational frequencies, electron affinities, ionization potentials, heats of formation, hydrogen bond interaction energies, conformational energies, and reaction barrier heights. Calculations are carried out with the 3-21G*, 6-31G*, 3-21+G*, 6-31+G*, 6-31++G*, cc-pVxZ, and aug-cc-pVxZ (x=D,T) basis sets, while bond distance and bond angle calculations are also done using the cc-pVQZ and aug-cc-pVQZ basis sets. Members of the popular functional classes, namely, LSDA, GGA, meta-GGA, hybrid-GGA, and hybrid-meta-GGA, are considered in this work. For the purpose of comparison, Hartree-Fock (HF) and second order many-body perturbation (MP2) methods are also assessed in terms of their ability to determine these physical properties. Ultimately, it is observed that the split valence bases of the 6-31G variety provide accuracies similar to those of the more computationally expensive Dunning type basis sets. Another conclusion from this survey is that the hybrid-meta-GGA functionals are typically among the most accurate functionals for all of the properties examined in this work.
294 citations
TL;DR: Several practical approaches are formulated to overcome difficulties providing a reliable description of electronic excitations in nanosystems.
Abstract: Time-dependent density functional theory (TDDFT) is a powerful tool allowing for accurate description of excited states in many nanoscale molecular systems; however, its application to large molecules may be plagued with difficulties that are not immediately obvious from previous experiences of applying TDDFT to small molecules. In TDDFT, the appearance of spurious charge-transfer states below the first optical excited state is shown to have significant effects on the predicted absorption and emission spectra of several donor-acceptor substituted molecules. The same problem affects the predictions of electronic spectra of molecular aggregates formed from weakly interacting chromophores. For selected benchmark cases, we show that today's popular density functionals, such as purely local (Local Density Approximation, LDA) and semilocal (Generalized Gradient Approximation, GGA) models, are qualitatively wrong. Nonlocal hybrid approximations including both semiempirical (B3LYP) and ab initio (PBE1PBE) containing a small fraction (20-25%) of Fock-like orbital exchange are also susceptible to such problems. Functionals that contain a larger fraction (50%) of orbital exchange like the early hybrid (BHandHLYP) are shown to exhibit far fewer spurious charge-transfer (CT) states at the expense of accuracy. Based on the trends observed in this study and our previous experience we formulate several practical approaches to overcome these difficulties providing a reliable description of electronic excitations in nanosystems.
292 citations
TL;DR: It is found that including environmental point charges can lower the errors in the electrostatically embedded pairwise additive energies for a series of water clusters by as much as a factor of 10 when compared to the traditional Pairwise additive approximation.
Abstract: The use of background molecular charge to incorporate environmental effects on a molecule or active site is widely employed in quantum chemistry. In the present article we employ this practice in conjunction with many-body expansions. In particular, we present electrostatically embedded two-body and three-body expansions for calculating the energies of molecular clusters. The system is divided into fragments, and dimers or trimers of fragments are calculated in a field of point charges representing the electrostatic potential of the other fragments. We find that including environmental point charges can lower the errors in the electrostatically embedded pairwise additive (EE-PA) energies for a series of water clusters by as much as a factor of 10 when compared to the traditional pairwise additive approximation and that for the electrostatically embedded three-body (EE-3B) method the average mean unsigned error over nine different levels of theory for a set of six tetramers and one pentamer is only 0.05 kc...
TL;DR: Empirical force field parameters consistent with the CHARMM additive and classical Drude based polarizable force fields are presented for linear and cyclic ethers and are shown to be in satisfactory agreement with both pure solvent and aqueous solvation properties.
Abstract: Empirical force field parameters consistent with the CHARMM additive and classical Drude based polarizable force fields are presented for linear and cyclic ethers. Initiation of the optimization pr...
TL;DR: Comparison of experimental and back-calculated NMR spin-relaxation data provides a more objective way of assessing the quality of the trajectory than order parameters alone, and shows significantly improved agreement with experimental NMR data for the AMber99SB force field as compared to AMBER99.
Abstract: Biological function of biomolecules is accompanied by a wide range of motional behavior. Accurate modeling of dynamics by molecular dynamics (MD) computer simulations is therefore a useful approach toward the understanding of biomolecular function. NMR spin relaxation measurements provide rigorous benchmarks for assessing important aspects of MD simulations, such as the amount and time scales of conformational space sampling, which are intimately related to the underlying molecular mechanics force field. Until recently, most simulations produced trajectories that exhibited too much dynamics particularly in flexible loop regions. Recent modifications made to the backbone φ and ψ torsion angle potentials of the AMBER and CHARMM force fields indicate that these changes produce more realistic molecular dynamics behavior. To assess the consequences of these changes, we performed a series of 5−20 ns molecular dynamics trajectories of human ubiquitin using the AMBER99 and AMBER99SB force fields for different con...
TL;DR: Force field parameters for the 107 modified nucleotides currently known to be present in RNA are developed and will improve the functionality of AMBER so that simulations can now be readily performed on diverse RNAs having post-transcriptional modifications.
Abstract: Classical molecular dynamics (MD) simulations are useful for characterizing the structure and dynamics of biological macromolecules, ultimately, resulting in elucidation of biological function. The AMBER force field is widely used and has well-defined bond length, bond angle, partial charge, and van der Waals parameters for all the common amino acids and nucleotides, but it lacks parameters for many of the modifications found in nucleic acids and proteins. Presently there are 107 known naturally occurring modifications that play important roles in RNA stability, folding, and other functions. Modified nucleotides are found in almost all transfer RNAs, ribosomal RNAs of both the small and large subunits, and in many other functional RNAs. We developed force field parameters for the 107 modified nucleotides currently known to be present in RNA. The methodology used for deriving the modified nucleotide parameters is consistent with the methods used to develop the Cornell et al. force field. These parameters will improve the functionality of AMBER so that simulations can now be readily performed on diverse RNAs having post-transcriptional modifications.
TL;DR: It is envisioned that the next generation of force fields for biomolecular polymer simulations will be developed based on electronic structure theory, which can adequately define and treat many-body polarization and charge delocalization effects.
Abstract: An electronic structure-based polarization method, called the X-POL potential, has been described for the purpose of constructing an empirical force field for modeling polypeptides. The X-POL potential is a quantum mechanical model, in which the internal, bonded interactions are fully represented by an electronic structure theory augmented with some empirical torsional terms. Nonbonded interactions are modeled by an iterative, combined quantum mechanical and molecular mechanical method, in which the molecular mechanical partial charges are derived from the molecular wave functions of the individual fragments. In this paper, the feasibility of such and electronic structure based force field is illustrated by small model compounds. A method has been developed for separating a polypeptide chain into peptide units, and its parametrization procedure in the X-POL potential is documented and tested on glycine dipeptide. We envision that the next generation of force fields for biomolecular polymer simulations wil...
TL;DR: Results show that the DFTB method with the present parameters reproduces structural properties very well, but the bond energies and the relative energies of different spin states only qualitatively compared to the B3LYP/SDD+6-31G(d) density functional (DFT) results.
Abstract: Recently developed parameters for five first-row transition-metal elements (M = Sc, Ti, Fe, Co, and Ni) in combination with H, C, N, and O as well as the same metal (M−M) for the spin-polarized self-consistent-charge density-functional tight-binding (DFTB) method have been calibrated. To test their performance a couple sets of compounds have been selected to represent a variety of interactions and bonding schemes that occur frequently in transition-metal containing systems. The results show that the DFTB method with the present parameters in most cases reproduces structural properties very well, but the bond energies and the relative energies of different spin states only qualitatively compared to the B3LYP/SDD+6-31G(d) density functional (DFT) results. An application to the ONIOM(DFT:DFTB) indicates that DFTB works well as the low level method for the ONIOM calculation.
TL;DR: A new open system Monte Carlo procedure designed to overcome difficulties with insertion and deletion of molecules is introduced, and is shown to yield correct results for the volumetric properties of the Lennard-Jones fluid and water as well as the phase behavior of the CO2-ethanol binary system.
Abstract: A new open system Monte Carlo procedure designed to overcome difficulties with insertion and deletion of molecules is introduced. The method utilizes gradual insertions and deletions of molecules through the use of a continuous coupling parameter and an adaptive bias potential. The method draws upon concepts from previous open system molecular dynamics and expanded ensemble Monte Carlo techniques and is applied to both the grand canonical and osmotic ensembles. It is shown to yield correct results for the volumetric properties of the Lennard-Jones fluid and water as well as the phase behavior of the CO2-ethanol binary system.
TL;DR: This work presents a molecular dynamics simulation algorithm which is multiscale in both time and space, and supplements the potential and kinetic energy expressions with auxiliary terms in order to recover the total energy as a conserved quantity, even when the total number of degrees of freedom changes during the simulation.
Abstract: Multiscale computer simulation algorithms are required to describe complex molecular systems with events occurring over a range of time and length scales. True multiscale simulations must solve the interface, or hand-shaking, problem of coupling together different levels of description in different spatial regions of the system. If the spatial regions of different resolution move over time, or if material is allowed to flow over the inter-region boundaries, a mechanism must be introduced into the multiscale algorithm to allow material to dynamically change its representation. While such a mechanism is highly desirable in many instances, it is fraught with technical difficulties. Here, we present a molecular dynamics simulation algorithm which is multiscale in both time and space. We supplement the potential and kinetic energy expressions with auxiliary terms in order to recover the total energy as a conserved quantity, even when the total number of degrees of freedom changes during the simulation. This is crucial for a proper assessment of the quality of adaptive hybrid algorithms, and in particular, it allows us to tune the hierarchy of RESPA levels to optimize the integration scheme.
TL;DR: The results strongly suggest that for accurate modeling of ions in biomolecular systems, great care should be taken in choosing balanced ionic parameters even when using the most popular force-fields.
Abstract: Realistic all-atom simulation of biological systems requires accurate modeling of both the biomolecules and their ionic environment. Recently, ion nucleation phenomena leading to the rapid growth of KCl or NaCl clusters in the vicinity of biomolecular systems have been reported. To better understand this phenomenon, molecular dynamics simulations of KCl aqueous solutions at three (1.0, 0.25, and 0.10 M) concentrations were performed. Two popular water models (TIP3P and SPC/E) and two Lennard-Jones parameter sets (AMBER and Dang) were combined to produce a total of 80 ns of molecular dynamics trajectories. Results suggest that the use of the Dang cation Lennard-Jones parameters instead of those adopted by the AMBER force-field produces a more accurate description of the ionic solution. In the later case, formation of salt aggregates is probably indicative of an artifact resulting from misbalanced force-field parameters. Because similar results were obtained with two different water parameter sets, the simulations exclude a water model dependency in the formation of anomalous ionic clusters. Overall, the results strongly suggest that for accurate modeling of ions in biomolecular systems, great care should be taken in choosing balanced ionic parameters even when using the most popular force-fields. These results invite a reexamination of older data obtained using available force-fields and a thorough check of the quality of current parameters sets by performing simulations at finite (>0.25 M) instead of minimal salt conditions.
TL;DR: A general "confine-and-release" framework for free energy calculations that accounts for free energies of conformational change is introduced and illustrated by demonstrating that an umbrella sampling protocol can obtain converged binding free energies that are independent of the starting protein structure and include these conformational changes free energies.
Abstract: Free energy calculations are increasingly being used to estimate absolute and relative binding free energies of ligands to proteins. However, computed free energies often appear to depend on the initial protein conformation, indicating incomplete sampling. This is especially true when proteins can change conformation on ligand binding, as free energies associated with these conformational changes are either ignored or assumed to be included by virtue of the sampling performed in the calculation. Here, we show that, in a model protein system (a designed binding site in T4 lysozyme), conformational changes can make a difference of several kcal/mol in computed binding free energies and that they are neglected in computed binding free energies if the system remains kinetically trapped in a particular metastable state on simulation timescales. We introduce a general “confine-and-release” framework for free energy calculations that accounts for these free energies of conformational change. We illustrate its use...
TL;DR: This work finds that not only does this new method lead to better energetics than the original EE-MB method but also that one is able to obtain excellent agreement with full MP2 calculations by considering only a two-body expansion of the correlation energy, leading to a considerable savings in computational time as compared to the three- body expansion.
Abstract: The electrostatically embedded many-body expansion (EE-MB), previously applied to the total electronic energy, is here applied only to the electronic correlation energy (CE), combined with a Hartree-Fock calculation on the entire system. The separate treatment of the Hartree-Fock and correlation energies provides an efficient way to approximate correlation energy for extended systems. We illustrate this here by calculating accurate Moller-Plesset second-order perturbation theory (MP2) energies for a series of clusters ranging in size from 5 to 20 water molecules. In this new method, called EE-MB-CE, where MB is pairwise additive (PA) or three-body (3B), the full Hartree-Fock energy of a system of N monomers is calculated (i.e., the many-body expansion is carried out to the Nth order), while the EE-MB method is used to calculate the correlation energy of the system. We find that not only does this new method lead to better energetics than the original EE-MB method but also that one is able to obtain excellent agreement with full MP2 calculations by considering only a two-body expansion of the correlation energy, leading to a considerable savings in computational time as compared to the three-body expansion. Additionally, we propose the use of a cutoff to further reduce the number of two-body terms that must be calculated, and we show that if a cutoff of 6 A is used, then one can eliminate up to 44% of the pairs and still calculate energies to within 0.1% of the net interaction energy of the full cluster.
TL;DR: The ability of the polarizable model to yield changes in dipole moment as well as the reproduction of a range of condensedphase properties indicates its utility in the study of the properties of alcohols in a variety of condensed phase environments aswell as representing an important step in the development of a comprehensive force field for biological molecules.
Abstract: A polarizable empirical force field based on the classical Drude oscillator has been developed for the aliphatic alcohol series. The model is optimized with an emphasis on condensed-phase properties and is validated against a variety of experimental data. Transferability of the developed parameters is emphasized by the use of a single electrostatic model for the hydroxyl group throughout the alcohol series. Aliphatic moiety parameters were transferred from the polarizable alkane parameter set, with only the Lennard-Jones parameters on the carbon in methanol optimized. The developed model yields good agreement with pure solvent properties with the exception of the heats of vaporization of 1-propanol and 1-butanol, which are underestimated by approximately 6%; special LJ parameters for the oxygen in these two molecules that correct for this limitation are presented. Accurate treatment of the free energies of aqueous solvation required the use of atom-type specific Oalcohol−Owater LJ interaction terms, with ...
TL;DR: A semiempirical AM1/d Hamiltonian is developed to model phosphoryl transfer reactions catalyzed by enzymes and ribozymes for use in linear-scaling calculations and combined quantum mechanical/molecular mechanical simulations and to facilitate the design of improved next-generation multiscale quantum models.
Abstract: A semiempirical AM1/d Hamiltonian is developed to model phosphoryl transfer reactions catalyzed by enzymes and ribozymes for use in linear-scaling calculations and combined quantum mechanical/molecular mechanical simulations. The model, designated AM1/d-PhoT, is parametrized for H, O, and P atoms to reproduce high-level density-functional results from a recently constructed database of quantum calculations for RNA catalysis (http://theory.chem.umn.edu/Database/QCRNA), including geometries and relative energies of minima, transition states and reactive intermediates, dipole moments, proton affinities, and other relevant properties. The model is tested in the gas phase and in solution using a QM/MM potential. The results indicate that the method provides significantly higher accuracy than MNDO/ d, AM1, and PM3 methods and, for the transphosphorylation reactions, is in close agreement with the density-functional calculations at the B3LYP/6-311++G(3df,2p) level with a reduction in computational cost of 3-4 orders of magnitude. The model is expected to have considerable impact on the application of semiempirical QM/MM methods to transphosphorylation reactions in solution, enzymes, and ribozymes and to ultimately facilitate the design of improved next- generation multiscale quantum models.
TL;DR: New scaling parameters are presented for use in the spin-component scaled (SCS) variant of density fitted local second-order Møller-Plesset perturbation theory (DF-LMP2) that have been optimized for use for evaluating the interaction energy between nucleic acid base pairs.
Abstract: New scaling parameters are presented for use in the spin-component scaled (SCS) variant of density fitted local second-order Moller−Plesset perturbation theory (DF-LMP2) that have been optimized for use in evaluating the interaction energy between nucleic acid base pairs. The optimal set of parameters completely neglects the contribution from antiparallel-spin electron pairs to the MP2 energy while scaling the parallel contribution by 1.76. These spin-component scaled for nucleobases (SCSN) parameters are obtained by minimizing, with respect to SCS parameters, the rms interaction energy error relative to the best available literature values, over a set of ten stacked nucleic acid base pairs. The applicability of this scaling to a wide variety of noncovalent interactions is verified through evaluation of a larger set of model complexes, including those dominated by dispersion and electrostatics.
TL;DR: The combination of the methods developed here presents a comprehensive and accurate treatment for the simulation of reaction processes in solution and in enzymes with ab initio QM/MM methods.
Abstract: Structural and energetic changes are two important characteristic properties of a chemical reaction process. In the condensed phase, studying these two properties is very challenging because of the great computational cost associated with the quantum mechanical calculations and phase space sampling. Although the combined quantum mechanics/molecular mechanics (QM/MM) approach significantly reduces the amount of the quantum mechanical calculations and facilitates the simulation of solution phase and enzyme catalyzed reactions, the required quantum mechanical calculations remain quite expensive and extensive sampling can be achieved routinely only with semiempirical quantum mechanical methods. QM/MM simulations with ab initio QM methods, therefore, are often restricted to narrow regions of the potential energy surface such as the reactant, product and transition state, or the minimum energy path. Such ab initio QM/MM calculations have previously been performed with the QM/MM-Free Energy (QM/MM-FE) method of Zhang et al.1 to generate the free energy profile along the reaction coordinate using free energy perturbation calculations at fixed structures of the QM subsystems. Results obtained with the QM/MM-FE method depend on the determination of the minimum energy reaction path, which is based on local conformations of the protein/solvent environment and can be difficult to obtain in practice. To overcome the difficulties associated with the QM/MM-FE method and to further enhance the sampling of the MM environment conformations, we develop here a new method to determine the QM/MM minimum free energy path (QM/MM-MFEP) for chemical reaction processes in solution and in enzymes. Within the QM/MM framework, we express the free energy of the system as a function of the QM conformation, thus leading to a simplified potential of mean force (PMF) description for the thermodynamics of the system. The free energy difference between two QM conformations is evaluated by the QM/MM free energy perturbation method. The free energy gradients with respect to the QM degrees of freedom are calculated from molecular dynamics simulations at given QM conformations. With the free energy and free energy gradients in hand, we further implement chain-of-conformation optimization algorithms in the search for the reaction path on the free energy surface without specifying a reaction coordinate. This method thus efficiently provides a unique minimum free energy path for solution and enzyme reactions, with structural and energetic properties being determined simultaneously. To further incorporate the dynamic contributions of the QM subsystem into the simulations, we develop the reaction path potential of Lu, et al.2 for the minimum free energy path. The combination of the methods developed here presents a comprehensive and accurate treatment for the simulation of reaction processes in solution and in enzymes with ab initio QM/MM methods. The method has been demonstrated on the first step of the reaction of the enzyme triosephosphate isomerase with good agreement with previous studies.
TL;DR: There are many minima on the potential energy surfaces of the methanol clusters, the number increasing rapidly with n, and a cyclic cluster of five to six methanols appears to be sufficient to mimic liquid behavior as far as vibrational frequencies are concerned.
Abstract: The potential energy surfaces of methanol clusters, (CH3OH)n, n = 2-12, have been studied using density functional theory at the B3LYP/6-31G(d) and higher levels of theory. Cyclic clusters in which n methanol molecules are joined in a ring structure formed by n hydrogen bonds are shown to be more stable than structures of the same number of methanol molecules where one or more methanol molecules are outside the ring and are hydrogen-bonded to oxygens of methanols in rings of n - 1, n - 2, and so forth. So-called chain structures are generally even less stable. Furthermore, the hydrogen-bonding energy per methanol molecule of the n-ring clusters is shown to converge to an asymptotic value of about 27 kJ/mol at B3LYP/6-311+G(d,p)//B3LYP/6-31G(d) after five to six methanols are included in the cluster. As expected, there are many minima on the potential energy surfaces of the methanol clusters, the number increasing rapidly with n. A cyclic cluster of five to six methanol molecules appears to be sufficient to mimic liquid behavior as far as vibrational frequencies are concerned.
TL;DR: The results indicate that the present approach is useful for studying the absorption spectra and the mechanism of the color tuning in the retinal proteins.
Abstract: The excited states of the three retinal proteins, bovine rhodopsin (Rh), bacteriorhodopsin (bR), and sensory rhodopsin II (sRII) were studied using the symmetry-adapted cluster-configuration interaction (SAC-CI) and combined quantum mechanical and molecular mechanical (QM/MM) methods. The computed absorption energies are in good agreement with the experimental ones for all three proteins. The spectral tuning mechanism was analyzed in terms of three contributions: molecular structures of the chromophore in the binding pockets, electrostatic (ES) interaction of the chromophore with the surrounding protein environment, and quantum-mechanical effect between the chromophore and the counterion group. This analysis provided an insight into the mechanism of the large blue-shifts in the absorption peak position of Rh and sRII from that of bR. Protein ES effect is primarily important both in Rh and in sRII, and the structure effect is secondary important in Rh. The quantum-mechanical interaction between the chromo...
TL;DR: A new general-purpose reactivity indicator is derived that can describe the reactivity of molecules that lie between the electrostatic (or charge) control and electron-transfer (or frontier-orbital) control paradigms.
Abstract: A new general-purpose reactivity indicator is derived. Unlike existing indicators, this indicator can describe the reactivity of molecules that lie between the electrostatic (or charge) control and electron-transfer (or frontier-orbital) control paradigms. Depending on the parameters in the indicator, it describes electrostatic control (where the electrostatic potential is the appropriate indicator), electron-transfer control (where the Fukui function's potential is the appropriate indicator), and intermediate cases (where linear combinations of the electrostatic potential and the Fukui function's potential are appropriate indicators). Our analysis gives some insight into the origins of the local hard/soft-acid/base principle. The "minimum Fukui function" rule for hard reagents also emerges naturally from our analysis: if (1) a reaction is strongly electrostatically controlled and (2) there are two sites that are equally favorable from an electrostatic standpoint, then the most reactive of the electrostatically equivalent sites is the site with the smallest Fukui function. An analogous electrostatic potential rule for soft reagents is also introduced: if (1) a reaction is strongly electron-transfer-controlled and (2) there are two sites where the Fukui function's potential are equivalent, then the most reactive of the Fukui-equivalent sites will be the one with greatest electrostatic potential (for electrophilic attack on a nucleophile) or smallest electrostatic potential (for nucleophilic attack on an electrophile).
TL;DR: A new method to fit the molecular electrostatic potentials obtained in quantum mechanical calculations to a set of classical electrostatic multipoles, usually point charges located at atomic positions, shows greatly improved numerical stability with respect to the molecular positions and geometries.
Abstract: We develop here a new method to fit the molecular electrostatic potentials obtained in quantum mechanical calculations to a set of classical electrostatic multipoles, usually point charges located at atomic positions. We define an object function of fitting as an integration of the difference of electrostatic potentials in the entire 3-dimensional physical space. The object function is thus rotationally invariant with respect to the molecular orientation and varies smoothly with respect to molecular geometric fluctuations. Compared with commonly employed methods such as the Merz−Singh−Kollman and CHELPG schemes, this new method, while possessing comparable accuracy, shows greatly improved numerical stability with respect to the molecular positions and geometries. The method can be used in the fitting of electrostatic potentials for the molecular mechanics force fields and also can be applied to the calculation of electrostatic polarizabilites of molecular or atomic systems.
TL;DR: A parallel coupled cluster algorithm that combines distributed and shared memory techniques for the CCSD(T) method (singles + doubles with perturbative triples) is described, targeted at modern cluster based architectures that are comprised of multiprocessor nodes connected by a dedicated communication network.
Abstract: A parallel coupled cluster algorithm that combines distributed and shared memory techniques for the CCSD(T) method (singles + doubles with perturbative triples) is described The implementation of the massively parallel CCSD(T) algorithm uses a hybrid molecular and “direct” atomic integral driven approach Shared memory is used to minimize redundant replicated storage per compute process The algorithm is targeted at modern cluster based architectures that are comprised of multiprocessor nodes connected by a dedicated communication network Parallelism is achieved on two levels: parallelism within a compute node via shared memory parallel techniques and parallelism between nodes using distributed memory techniques The new parallel implementation is designed to allow for the routine evaluation of mid- (500−750 basis function) to large-scale (750−1000 basis function) CCSD(T) energies Sample calculations are performed on five low-lying isomers of water hexamer using the aug-cc-pVTZ basis set