Search

We use local diffusion maps to assess the quality of two types of collective variables (CVs) for a recently published hydrogen combustion benchmark dataset~\cite{guan2022benchmark} that contains ab initio molecular dynamics trajectories and normal modes along minimum energy paths. This approach was recently advocated in~\cite{tlldiffmap20} for assessing CVs and analyzing reactions modeled by classical molecular dynamics simulations. We report the effectiveness of this approach to molecular systems modeled by quantum ab initio molecular dynamics. In addition to assessing the quality of CVs, we also use global diffusion maps to perform committor analysis as proposed in~\cite{tlldiffmap20}. We show that the committor function obtained from the global diffusion map allows us to identify transition regions of interest in several hydrogen combustion reaction channels.

We present an investigation into the transferability of pseudopotentials (PPs) with a nonlinear core correction (NLCC) using the Goedecker, Teter, and Hutter (GTH) protocol across a range of pure GGA, meta-GGA and hybrid functionals, and their impact on thermochemical and non-thermochemical properties. The GTH-NLCC PP for the PBE density functional demonstrates remarkable transferability to the PBE0 and $\omega$B97X-V exchange-correlation functionals, and relative to no NLCC, improves agreement significantly for thermochemical benchmarks compared to all-electron calculations. On the other hand, the B97M-rV meta-GGA functional performs poorly with the PBE-derived GTH-NLCC PP, which is mitigated by reoptimizing the NLCC parameters for this specific functional. The findings reveal that atomization energies exhibit the greatest improvements from use of the NLCC, which thus provides an important correction needed for covalent interactions relevant to applications involving chemical reactivity. Finally we test the NLCC-GTH PPs when combined with medium-size TZV2P molecularly optimized (MOLOPT) basis sets which are typically utilized in condensed phase simulations, and show that they lead to consistently good results when compared to all-electron calculations for atomization energies, ionization potentials, barrier heights, and non-covalent interactions, but lead to somewhat larger errors for electron affinities.

Theoretical predictions of NMR chemical shifts from first-principles can greatly facilitate experimental interpretation and structure identification. However, accurate prediction of chemical shifts using the best coupled cluster methods can be prohibitively expensive for systems larger than ten to twenty non-hydrogen atoms on today's computers. By contrast machine learning methods offer inexpensive alternatives but are hampered by generalization to molecules outside the original training set. Here we propose a novel machine learning feature representation informed by intermediate calculations of atomic chemical shielding tensors within a molecular environment using an inexpensive quantum mechanics method, and training it to predict NMR chemical shieldings of a high-level composite theory that is comparable to CCSD(T) in the complete basis set limit. The inexpensive shift machine learning (iShiftML) algorithm is trained through a new progressive active learning workflow that reduces the total number of expensive calculations required when constructing the dataset, while allowing the model to continuously improve on data it has never seen. Furthermore, we show that the error estimations from our model correlate quite well with actual errors to provide confidence values on new predictions. We illustrate the predictive capacity of iShiftML across gas phase experimental chemical shifts for small organic molecules and much larger and more complex natural products in which we can accurately differentiate between subtle diastereomers based on chemical shift assignments.

The electrostatic potential (ESP) is a powerful property for understanding and predicting electrostatic charge distributions that drive interactions between molecules. In this study, we compare various charge partitioning schemes including fitted charges, density-based quantum mechanical (QM) partitioning schemes, charge equilibration methods, and our recently introduced coarse-grained electron model, C-GeM, to describe the ESP for protein systems. When benchmarked against high quality density functional theory calculations of the ESP for tripeptides and the crambin protein, we find that the C-GeM model is of comparable accuracy to ab initio charge partitioning methods, but with orders of magnitude improvement in computational efficiency since it does not require either the electron density or the electrostatic potential as input.

Reactive force fields provide an affordable model for simulating chemical reactions at a fraction of the cost of quantum mechanical approaches. However, classically accounting for chemical reactivity often comes at the expense of accuracy and transferability, while computational cost is still large relative to nonreactive force fields. In this Perspective, we summarize recent efforts for improving the performance of reactive force fields in these three areas with a focus on the ReaxFF theoretical model. To improve accuracy, we describe recent reformulations of charge equilibration schemes to overcome unphysical long-range charge transfer, new ReaxFF models that account for explicit electrons, and corrections for energy conservation issues of the ReaxFF model. To enhance transferability we also highlight new advances to include explicit treatment of electrons in the ReaxFF and hybrid nonreactive/reactive simulations that make it possible to model charge transfer, redox chemistry, and large systems such as reverse micelles within the framework of a reactive force field. To address the computational cost, we review recent work in extended Lagrangian schemes and matrix preconditioners for accelerating the charge equilibration method component of ReaxFF and improvements in its software performance in LAMMPS.

Density Functional Theory (DFT) is currently the most tractable choice of theoretical model used to understand the mechanistic pathways for electrocatalytic processes such as CO2 or CO reduction. Here, we assess the performance of two DFT functionals designed specifically to describe surface interactions, RTPSS and RPBE, as well as two popular meta-GGA functionals, SCAN and B97M-rV, that have not been a priori optimized for better interfacial properties. We assess all four functionals against available experimental data for prediction of bulk and bare surface properties on four electrocatalytically relevant metals, Au, Ag, Cu, and Pt, and for binding CO to surfaces of these metals. To partially mitigate issues such as thermal and anharmonic corrections associated with comparing computations with experiments, molecular benchmarks against high level quantum chemistry are reported for CO complexes with Au, Ag, Cu and Pt atoms, as well as the CO-water complex and the water dimer. Overall, we find that the surface modified RPBE functional performs reliably for many of the benchmarks examined here, and the meta-GGA functionals also show promising results. Specifically B97M-rV predicts the correct site preference for CO binding on Ag and Au (the only functional tested here to do so), while RTPSS performs well for surface relaxations and binding of CO on Pt and Cu.

Magnetic properties of molecules such as magnetizabilities represent second order derivatives of the energy with respect to external perturbations. To avoid the need for analytic second derivatives and thereby permit evaluation of the performance of methods where they are not available, a new implementation of quantum chemistry calculations in finite applied magnetic fields is reported. This implementation is employed for a collection of small molecules with the aug-cc-pVTZ basis set to assess orbital optimized (OO) MP2 and a recently proposed regularized variant of OOMP2, called κ-OOMP2. κ-OOMP2 performs significantly better than conventional second order Møller-Plesset (MP2) theory, by reducing MP2's exaggeration of electron correlation effects. As a chemical application, we revisit an old aromaticity criterion called magnetizability exaltation. In lieu of empirical tables or increment systems to generate references, we instead use straight chain molecules with the same formal bond structure as the target cyclic planar conjugated molecules. This procedure is found to be useful for qualitative analysis, yielding exaltations that are typically negative for aromatic species and positive for antiaromatic molecules. One interesting species, N2S2, shows a positive exaltation despite having aromatic characteristics.

We present a new stochastic extended Lagrangian molecular dynamics solution to charge equilibration that eliminates self-consistent field (SCF) calculations, thus eliminating the computational bottleneck in solving the charge distribution with standard SCF solvers. By formulating both charges and chemical potential as latent variables and introducing a holonomic constraint that satisfies charge conservation, the SC-XLMD method accurately reproduces thermodynamic, dynamic, and structural properties within the framework of ReaxFF for a bulk water system and highly reactive RDX molecules simulated at high temperature. The SC-XLMD method shows excellent computational performance and is available in the publicly available LAMMPS package.

We report the adiabatic energy decomposition analysis (EDA) of density functional theory (DFT) results, shedding light on the physical content of binding energies and carbon monoxide (CO) frequency (υCO) shifts in select first-row transition metal monocarbonyls (MCOs; M = Ti-, V-, Cr-, Co-, Ni-, Cu-, V, Cr, Mn, Ni, Cu, Zn, Cr+, Mn+, Fe+, Cu+, and Zn+). This approach allows for the direct decomposition of υCO, in contrast to previous studies of these systems. Neutral, anionic, and cationic systems are compared, and our results indicate that the relative importance of electrostatic interactions, intramolecular orbital polarization, and charge transfer can vary significantly with the charge and electron configuration of the metal participating in binding. Various anomalous systems are also discussed and incorporated into a general model of MCO binding. Electrostatic interactions and orbital polarization are found to promote blue shifts in υCO, while charge transfer effects encourage υCO red-shifting; previously reported values of υCO are found to be a result of a complex but quantifiable interplay between these physical components. Our computations indicate that CuCO- and ZnCO possess triplet ground states, and also that CrCO- exhibits a non-linear geometry, all in contrast to previous computational results. Advantages and limitations of this model as an approximation to more complicated systems, like those implicated in heterogeneous catalysis, are discussed. We also report benchmark results for MCO geometries, binding energies, and harmonic CO frequencies, and discuss the validity of single-reference wave function and DFT approaches to the study of these transition metal systems.

To study intermolecular interactions involving radicals at the correlated level, the energy decomposition analysis scheme for second-order Mo̷ller-Plesset perturbation theory based on absolutely localized molecular orbitals (ALMO-MP2-EDA) is generalized to unrestricted and restricted open-shell MP2. The benefit of restricted open-shell MP2 is that it can provide accurate binding energies for radical complexes where density functional theory can be error-prone due to delocalization errors. As a model application, the open-shell ALMO-MP2-EDA is applied to study the first solvation step of halogenated benzene radical cations, where both halogen- and hydrogen-bonded isomers are possible. We determine that the lighter halogens favor the hydrogen-bonded form, while the iodine-substituted species prefers halogen bonding due to larger polarizability and charge transfer at the halogen. As a second application, relevant to the activation of CO2 in photoelectrocatalysis, complexes of CO2-· interacting with both pyridine and imidazole are analyzed with ALMO-MP2-EDA. The results reveal the importance of charge transfer into the π* orbital of the heterocycle in controlling the stability of the carbamate binding mode, which is favored for pyridine but not for imidazole.

We review the standard model for de novo computational design of enzymes, which primarily focuses on the development of an active-site geometry, composed of protein functional groups in orientations optimized to stabilize the transition state, for a novel chemical reaction not found in nature. Its emphasis is placed on the structure and energetics of the active site embedded in an accommodating protein that serves as a physical support that shields the reaction chemistry from solvent, which is typically improved upon by laboratory-directed evolution. We also provide a review of design strategies that move beyond the standard model, by placing more emphasis on the designed enzyme as a whole catalytic construct. Starting with complete de novo enzyme design examples, we consider additional design factors such as entropy of individual residues, correlated motion between side chains (mutual information), dynamical correlations of the enzyme motions that could aid the reaction, reorganization energy, and electric fields as ways to exploit the entire protein scaffold to improve upon the catalytic rate, thereby providing directed evolution with better starting sequences for increasing biocatalytic performance.

We have developed a deep learning algorithm for chemical shift prediction for atoms in molecular crystals that utilizes an atom-centered Gaussian density model for the 3D data representation of a molecule. We define multiple channels that describe different spatial resolutions for each atom type that utilizes cropping, pooling, and concatenation to create a multi-resolution 3D-DenseNet architecture (MR-3D-DenseNet). Because the training and testing time scale linearly with the number of samples, the MR-3D-DenseNet can exploit data augmentation that takes into account the property of rotational invariance of the chemical shifts, thereby also increasing the size of the training dataset by an order of magnitude without additional cost. We obtain very good agreement for 13C, 15N, and 17O chemical shifts, with the highest accuracy found for 1H chemical shifts that is equivalent to the best predictions using ab initio quantum chemistry methods.

Alzheimer's disease (AD) is a progressive neurodegenerative disease affecting millions of people worldwide. AD is characterized by the presence of extracellular plaques composed of aggregated/oligomerized β-amyloid peptides with Aβ42 peptide representing a major isoform in the senile plaques. Given the pathological significance of Aβ42 in the progression of AD, there is considerable interest in understanding the structural ensembles for soluble monomer and oligomeric forms of Aβ42. This report describes an efficient method to express and purify high quality (15)N isotope-labeled Aβ42 for structural studies by NMR. The protocol involves utilization of an auto induction system with (15)N isotope labeled medium, for high-level expression of Aβ42 as a fusion with IFABP. After the over-expression of the (15)N isotope-labeled IFABP-Aβ42 fusion protein in the inclusion bodies, pure (15)N isotope-labeled Aβ42 peptide is obtained following a purification method that is streamlined and improved from the method originally developed for the isolation of unlabeled Aβ42 peptide (Garai et al., 2009). We obtain a final yield of ∼ 6 mg/L culture for (15)N isotope-labeled Aβ42 peptide. Mass spectrometry and (1)H-(15)N HSQC spectra of monomeric Aβ42 peptide validate the uniform incorporation of the isotopic label. The method described here is equally applicable for the uniform isotope labeling with (15)N and (13)C in Aβ42 peptide as well as its other variants including any Aβ42 peptide mutants.

We present a reinvestigation of sulfate-water clusters SO4(2-) (H2O)(n=3-7), which involves several new aspects. Using a joint molecular mechanics/first principles approach, we perform exhaustive searches for stable cluster geometries, showing that the sulfate-water landscape is much richer than anticipated previously. We check the compatibility of the new structures with experiment by comparing vertical detachment energies (VDEs) calculated at the B3LYP/6-311++G** level of theory and determine the energetic ordering of the isomers at the RI-MP2/aug-cc-pVTZ level. Our results are bench-marked carefully against reference energies of estimated CCSD(T)/aug-cc-VTZ quality and VDEs of CCSD(T)/aug-cc-pVDZ quality. Furthermore, we calculate anharmonic vibrational corrections for up to the n = 6 clusters, which are shown to be significant for isomer energy ordering. We use energy decomposition analysis (EDA) based on the absolutely localized fragment (ALMO) expansion to gain chemical insight into the binding motifs.

We present the simulated photoelectron spectrum (PES) for cyanide-water CN(H2O)-based on quasiclassical trajectory molecular dynamics (QCT-MD). Using density functional theory to generate trajectories and to calculate vertical detachment energies, we obtain simulated spectra that are in qualitative agreement with experiment. We obtain a theoretical 12 → 300 K temperature red shift of 0.1 eV as compared to an experimental redshift of 0.25 eV. The calculated linewidths of 0.3 eV are in excellent agreement with experiment. Our trajectories show that the temperature red shift as being dominated by dynamics within the basin of the N-bound minimum, however, at 300 K we predict conversion into the basin of the C-bound minimum, equilibrating at a 80:20 ratio of N- vs C-bound mixture. We discuss the potential advantages of QCT-MD over anharmonic Franck-Condon analysis such as natural incorporation of anharmonicity (as necessary for weakly bound systems), and reduced computational scaling, but also drawbacks such as neglect of final-state (e.g., Duschinsky) effects. © 2010 American Chemical Society.

We present what we believe to be a novel statistical contact potential based on solved structures of transmembrane (TM) alpha-helical bundles, and we use this contact potential to investigate the amino acid likelihood of stabilizing helix-helix interfaces. To increase statistical significance, we have reduced the full contact energy matrix to a four-flavor alphabet of amino acids, automatically determined by our methodology, in which we find that polarity is a more dominant factor of group identity than is size, with charged or polar groups most often occupying the same face, whereas polar/apolar residue pairs tend to occupy opposite faces. We found that the most polar residues strongly influence interhelical contact formation, although they occur rarely in TM helical bundles. Two-body contact energies in the reduced letter code are capable of determining native structure from a large decoy set for a majority of test TM proteins, at the same time illustrating that certain higher-order sequence correlations are necessary for more accurate structure predictions.

We derive a new numerical approach to solving the linearized Poisson Boltzmann equation (PBE) by representing the protein surface as a collection of spheres in which the surface charges can then be iteratively solved by new analytical multipole methods previously introduced by us [Lotan & Head-Gordon, 2006]. We show that our Poisson Boltzmann semi-analytical method, PB-SAM, realizes better accuracy, more flexible memory management, and at reduced cost relative to either finite difference or boundary element method PBE solvers. We provide two new benchmarks of PBE solution accuracy to test the numerical PBE solutions based on (1) arrays of up to hundreds of spherical low dielectric geometries with asymmetric charges in which mutual polarization is treated exactly, and (2) two overlapping spheres with increasing charge asymmetry by solving the PB-SAM method to very high pole order. We illustrate the strength of the PB-SAM approach by computing the potential profile of an array of 60 T1-particle forming monomers of the bromine mosaic virus.

We combine molecular dynamics simulations and new high-field NMR experiments to describe the solution structure of the Abeta(21-30) peptide fragment that may be relevant for understanding structural mechanisms related to Alzheimer's disease. By using two different empirical force-field combinations, we provide predictions of the three-bond scalar coupling constants ((3)J(H(N)H(alpha))), chemical-shift values, (13)C relaxation parameters, and rotating-frame nuclear Overhauser effect spectroscopy (ROESY) crosspeaks that can then be compared directly to the same observables measured in the corresponding NMR experiment of Abeta(21-30). We find robust prediction of the (13)C relaxation parameters and medium-range ROESY crosspeaks by using new generation TIP4P-Ew water and Amber ff99SB protein force fields, in which the NMR validates that the simulation yields both a structurally and dynamically correct ensemble over the entire Abeta(21-30) peptide. Analysis of the simulated ensemble shows that all medium-range ROE restraints are not satisfied simultaneously and demonstrates the structural diversity of the Abeta(21-30) conformations more completely than when determined from the experimental medium-range ROE restraints alone. We find that the structural ensemble of the Abeta(21-30) peptide involves a majority population (approximately 60%) of unstructured conformers, lacking any secondary structure or persistent hydrogen-bonding networks. However, the remaining minority population contains a substantial percentage of conformers with a beta-turn centered at Val24 and Gly25, as well as evidence of the Asp23 to Lys28 salt bridge important to the fibril structure. This study sets the stage for robust theoretical work on Abeta(1-40) and Abeta(1-42), for which collection of detailed NMR data on the monomer will be more challenging because of aggregation and fibril formation on experimental timescales at physiological conditions. In addition, we believe that the interplay of modern molecular simulation and high-quality NMR experiments has reached a fruitful stage for characterizing structural ensembles of disordered peptides and proteins in general.

Using a coarse-grained model of the Abeta peptide, we analyze the Arctic (E22G), Dutch (E22Q), and Flemish (A21G) familial Alzheimer's disease (FAD) mutants for any changes in the stability of amyloid assemblies with respect to the wild-type (WT) sequence. Based on a structural reference state of two protofilaments aligned to create the "agitated" protofibril as determined by solid-state NMR, we determine free energy trends for Abeta assemblies for the WT and FAD familial sequences. We find that the structural characteristics and oligomer size of the critical nucleus vary dramatically among the hereditary mutants. The Arctic mutant's disorder in the turn region introduces new stabilizing interactions that better align the two protofilaments, yielding a well-defined protofibril axis at relatively small oligomer sizes with respect to WT. By contrast, the critical nucleus for the Flemish mutant is beyond the 20 chains characterized in this study, thereby showing a strong shift in the equilibrium toward monomers with respect to larger protofibril assemblies. The Dutch mutant forms more ordered protofilaments than WT, but exhibits greater disorder in protofibril structure that includes an alternative polymorph of the WT fibril. An important conclusion of this work is that the Dutch mutant does not support the agitated protofibril assembly. We discuss the implications of the structural ensembles and free energy profiles for the FAD mutants in regards to interpretation of the kinetics of fibril assembly using chromatography and dye-binding experiments.

We develop a sequence based alpha-carbon model to incorporate a mean field estimate of the orientation dependence of the polypeptide chain that gives rise to specific hydrogen bond pairing to stabilize alpha-helices and beta-sheets. We illustrate the success of the new protein model in capturing thermodynamic measures and folding mechanism of proteins L and G. Compared to our previous coarse-grained model, the new model shows greater folding cooperativity and improvements in designability of protein sequences, as well as predicting correct trends for kinetic rates and mechanism for proteins L and G. We believe the model is broadly applicable to other protein folding and protein-protein co-assembly processes, and does not require experimental input beyond the topology description of the native state. Even without tertiary topology information, it can also serve as a mid-resolution protein model for more exhaustive conformational search strategies that can bridge back down to atomic descriptions of the polypeptide chain.

The molecular mechanism of the solvent motion that is required to instigate the protein structural relaxation above a critical hydration level or transition temperature has yet to be determined. In this work we use quasi-elastic neutron scattering (QENS) and molecular dynamics simulation to investigate hydration water dynamics near a greatly simplified protein surface. We consider the hydration water dynamics near the completely deuterated N-acetyl-leucine-methylamide (NALMA) solute, a hydrophobic amino acid side chain attached to a polar blocked polypeptide backbone, as a function of concentration between 0.5M-2.0M, under ambient conditions. In this Communication, we focus our results of hydration dynamics near a model protein surface on the issue of how enzymatic activity is restored once a critical hydration level is reached, and provide a hypothesis for the molecular mechanism of the solvent motion that is required to trigger protein structural relaxation when above the hydration transition., Comment: 2 pages, 2 figures, Communication