Search

The Hessian matrix of the potential energy of molecular systems is employed not only in geometry optimizations or high-order molecular dynamics integrators but also in many other molecular procedures, such as instantaneous normal mode analysis, force field construction, instanton calculations, and semiclassical initial value representation molecular dynamics, to name a few. Here, we present an algorithm for the calculation of the approximated Hessian in molecular dynamics. The algorithm belongs to the family of unsupervised machine learning methods, and it is based on the neural gas idea, where neurons are molecular configurations whose Hessians are adopted for groups of molecular dynamics configurations with similar geometries. The method is tested on several molecular systems of different dimensionalities both in terms of accuracy and computational time versus calculating the Hessian matrix at each time-step, that is, without any approximation, and other Hessian approximation schemes. Finally, the method is applied to the on-the-fly, full-dimensional simulation of a small synthetic peptide (the 46 atom N-acetyl-l-phenylalaninyl-l-methionine amide) at the level of DFT-B3LYP-D/6-31G* theory, from which the semiclassical vibrational power spectrum is calculated.

The pair coupled cluster doubles (pCCD) method (where the excitation manifold is restricted to electron pairs) has a series of interesting features. Among others, it provides ground-state energies very close to what is obtained with doubly-occupied configuration interaction (DOCI), but with polynomial cost (compared with the exponential cost of the latter). Here, we address whether this similarity holds for excited states, by exploring the symmetric dissociation of the linear \ce{H4} molecule. When ground-state Hartree-Fock (HF) orbitals are employed, pCCD and DOCI excited-state energies do not match, a feature that is assigned to the poor HF reference. In contrast, by optimizing the orbitals at the pCCD level (oo-pCCD) specifically for each excited state, the discrepancies between pCCD and DOCI decrease by one or two orders of magnitude. Therefore, the pCCD and DOCI methodologies still provide comparable energies for excited states, but only if suitable, state-specific orbitals are adopted. We also assessed whether a pCCD approach could be used to directly target doubly-excited states, without having to resort to the equation-of-motion (EOM) formalism. In our $\Delta$oo-pCCD model, excitation energies were extracted from the energy difference between separate oo-pCCD calculations for the ground state and the targeted excited state. For a set comprising the doubly-excited states of \ce{CH+}, \ce{BH}, nitroxyl, nitrosomethane, and formaldehyde, we found that $\Delta$oo-pCCD provides quite accurate excitation energies, with root mean square deviations (with respect to full configuration interaction results) lower than CC3 and comparable to EOM-CCSDT, two methods with much higher computational cost., Comment: 12 pages, 4 figures

A linear-scaling local second-order Møller-Plesset (MP2) method is presented for high-spin open-shell molecules based on restricted open-shell (RO) reference functions. The open-shell local MP2 (LMP2) approach inherits the iteration- and redundancy-free formulation and the completely integral-direct, OpenMP-parallel, and memory and disk use economic algorithms of our closed-shell LMP2 implementation. By utilizing restricted local molecular orbitals for the demanding integral transformation step and by introducing a novel long-range spin-polarization approximation, the computational cost of RO-LMP2 approaches that of closed-shell LMP2. Extensive benchmarks were performed for reactions of radicals, ionization potentials, as well as spin-state splittings of carbenes and transition-metal complexes. Compared to the conventional MP2 reference for systems of up to 175 atoms, local errors of at most 0.1 kcal/mol were found, which are well below the intrinsic accuracy of MP2. RO-LMP2 computations are presented for challenging protein models of up to 601 atoms and 11 000 basis functions, which involve either spin states of a complexed iron ion or a highly delocalized singly occupied orbital. The corresponding runtimes of 9-15 h obtained with a single, many-core CPU demonstrate that MP2, as well as spin-scaled MP2 and double-hybrid density functional methods, become widely accessible for open-shell systems of unprecedented size and complexity.

An implementation of real-time time-dependent Hartree-Fock (RT-TDHF) and current-density-functional theory (RT-TDCDFT) for molecules in strong uniform magnetic fields is presented. In contrast to earlier implementations, the present work enables the use of the RT-TDCDFT formalism, which explicitly includes field dependent terms in the exchange-correlation functional. A range of current dependent exchange-correlation functionals based on the TPSS functional are considered, including a range-separated variant, which is particularly suitable for application to excited state calculations. The performance of a wide range of propagator algorithms for real-time methods is investigated in this context. A recently proposed molecular orbital pair decomposition analysis allows for assignment of electronic transitions, providing detailed information about which molecular orbitals are involved in each excitation. 1 The application of these methods is demonstrated for the electronic absorption spectra of N 2 and H 2 O both in the absence and in the presence of a magnetic field. The dependence of electronic spectra on the magnetic field strength and its orientation relative to the molecule is studied. The complex evolution of the absorption spectra with magnetic field is rationalised using the molecular orbital pair decomposition analysis, which provides crucial insight in strong fields where the spectra are radically different from their zero-field counterparts.

Protein-embedded chromophores are responsible for light harvesting, excitation energy transfer, and charge separation in photosynthesis. A critical part of the photosynthetic apparatus are reaction centers (RCs), which comprise groups of (bacterio)chlorophyll and (bacterio)pheophytin molecules that transform the excitation energy derived from light absorption into charge separation. The lowest excitation energies of individual pigments (site energies) are key for understanding photosynthetic systems, and form a prime target for quantum chemistry. A major theoretical challenge is to accurately describe the electrochromic (Stark) shifts in site energies produced by the inhomogeneous electric field of the protein matrix. Here, we present large-scale quantum mechanics/molecular mechanics calculations of electrochromic shifts for the RC chromophores of photosystem II (PSII) using various quantum chemical methods evaluated against the domain-based local pair natural orbital (DLPNO) implementation of the similarity-transformed equation of motion coupled cluster theory with single and double excitations (STEOM-CCSD). We show that certain range-separated density functionals (ωΒ97, ωΒ97X-V, ωΒ2PLYP, and LC-BLYP) correctly reproduce RC site energy shifts with time-dependent density functional theory (TD-DFT). The popular CAM-B3LYP functional underestimates the shifts and is not recommended. Global hybrid functionals are too insensitive to the environment and should be avoided, while nonhybrid functionals are strictly nonapplicable. Among the applicable approximate coupled cluster methods, the canonical versions of CC2 and ADC(2) were found to deviate significantly from the reference results both for the description of the lowest excited state and for the electrochromic shifts. By contrast, their spin-component-scaled (SCS) and particularly the scale-opposite-spin (SOS) variants compare well with the reference DLPNO-STEOM-CCSD and the best range-separated DFT methods. The emergence of RC excitation asymmetry is discussed in terms of intrinsic and protein electrostatic potentials. In addition, we evaluate a minimal structural scaffold of PSII, the D1-D2-CytB559 RC complex often employed in experimental studies, and show that it would have the same site energy distribution of RC chromophores as the full PSII supercomplex, but only under the unlikely conditions that the core protein organization and cofactor arrangement remain identical to those of the intact enzyme.

Funding Information: This work was supported by the Academy of Finland (project numbers 308647, 314298, and 316601) and through their Flagship program: Finnish Center for Artificial Intelligence FCAI. We thank CSC, the Finnish IT Center for Science, and Aalto Science IT for computational resources. This work is supported by COST (European Cooperation in Science and Technology) Action 18234. L.F. thanks Guoxu Zhang, Marc Dvorak, Jingrui Li, and Annika Stuke for the help with FHI-Aims. He also acknowledges financial support from the Chinese Scholarship Council (grant no. [2017]3109). 1 Publisher Copyright: © 2020 American Chemical Society. All rights reserved. Copyright: Copyright 2021 Elsevier B.V., All rights reserved. Finding low-energy molecular conformers is challenging due to the high dimensionality of the search space and the computational cost of accurate quantum chemical methods for determining conformer structures and energies. Here, we combine active-learning Bayesian optimization (BO) algorithms with quantum chemistry methods to address this challenge. Using cysteine as an example, we show that our procedure is both efficient and accurate. After only 1000 single-point calculations and approximately 80 structure relaxations, which is less than 10% computational cost of the current fastest method, we have found the low-energy conformers in good agreement with experimental measurements and reference calculations. To test the transferability of our method, we also repeated the conformer search of serine, tryptophan, and aspartic acid. The results agree well with previous conformer search studies.

The popular MARTINI coarse-grained model is used as a test case to analyze the adherence of top-down coarse-grained molecular dynamics models (i.e., models primarily parametrized to match experimental results) to the known features of statistical mechanics for the underlying all-atom representations. Specifically, the temperature dependence of various pair distribution functions, and hence their underlying potentials of mean force via the reversible work theorem, are compared between MARTINI 2.0, Dry MARTINI, and all-atom simulations mapped onto equivalent coarse-grained sites for certain lipid bilayers. It is found that the MARTINI models do not completely capture the lipid structure seen in atomistic simulations as projected onto the coarse-grained mappings and that issues of accuracy and temperature transferability arise due to an incorrect enthalpy-entropy decomposition of these potentials of mean force. The potential of mean force for the association of two amphipathic helices in a lipid bilayer is also calculated, and especially at shorter ranges, the MARTINI and all-atom projection results differ substantially. The former is much less repulsive and hence will lead to a higher probability of MARTINI helix association in the MARTINI bilayer than occurs in the actual all-atom case. Additionally, the bilayer height fluctuation spectra are calculated for the MARTINI model, and compared to the all-atom results, it is found that the magnitude of thermally averaged amplitudes at intermediate length scales are quite different, pointing to a number of possible consequences for realistic modeling of membrane processes. Taken as a whole, the results presented here show disagreement in the enthalpic and entropic driving forces driving lateral structure in lipid bilayers as well as quantitative differences in association of embedded amphipathic helices, which can help direct future efforts to parametrize CG models with better agreement to the all-atom systems they aspire to represent.

A simple and robust range-separated (RS) double-hybrid (DH) time-dependent density functional approach is presented for the accurate calculation of excitation energies of molecules within the Tamm-Dancoff approximation. The scheme can be considered as an excited-state extension of the ansatz proposed by Toulouse and co-workers [J. Chem. Phys. 2018, 148, 164105], which is based on the two-parameter decomposition of the Coulomb potential, for which both the exchange and correlation contributions are range-separated. A flexible and efficient implementation of the new scheme is also presented, which facilitates its extension to any combination of exchange and correlation functionals. The performance of the new approximation is tested for singlet excitations on several benchmark compilations and thoroughly compared to that of representative DH, RS hybrid, and RS DH functionals. The one-electron basis set dependence and computation times are also assessed. Our results show that the new approach improves on standard DHs in most cases, and it can provide a more robust and accurate alternative. In addition, on average, it noticeably surpasses the existing RS hybrid and RS DH functionals.

Thermodynamic integration (TI) is a commonly used method to determine free-energy differences. One of its disadvantages is that many intermediate λ-states need to be sampled in order to be able to integrate accurately over ⟨∂H/∂λ⟩. Here, we use the recently introduced extended TI to study alternative parameterizations of H(λ) and its influence on the smoothness of the ⟨∂H/∂λ⟩ curves as well as the efficiency of the simulations. We find that the extended TI approach can be used to select curves of low curvature. An optimal parameterization is suggested for the calculation of hydration free energies. For calculations of relative binding free energies, we show that optimized parameterizations of the Hamiltonian in the unbound state also effectively lower the curvature in the bound state of the ligand.

Because of their anisotropic electron distribution and electron deficiency, halonium ions are unusually strong halogen-bond donors that form strong and directional three-center, four-electron halogen bonds. These halogen bonds have received considerable attention owing to their applicability in supramolecular and synthetic chemistry and have been intensely studied using spectroscopic and crystallographic techniques over the past decade. Their computational treatment faces different challenges to those of conventional weak and neutral halogen bonds. Literature studies have used a variety of wave functions and DFT functionals for prediction of their geometries and NMR chemical shifts, however, without any systematic evaluation of the accuracy of these methods being available. In order to provide guidance for future studies, we present the assessment of the accuracy of 12 common DFT functionals along with the Hartree–Fock (HF) and the second-order Møller–Plesset perturbation theory (MP2) methods, selected from an initial set of 36 prescreened functionals, for the prediction of 1H, 13C, and 15N NMR chemical shifts of [N–X–N]+ halogen-bond complexes, where X = F, Cl, Br, and I. Using a benchmark set of 14 complexes, providing 170 high-quality experimental chemical shifts, we show that the choice of the DFT functional is more important than that of the basis set. The M06 functional in combination with the aug-cc-pVTZ basis set is demonstrated to provide the overall most accurate NMR chemical shifts, whereas LC-ωPBE, ωB97X-D, LC-TPSS, CAM-B3LYP, and B3LYP to show acceptable performance. Our results are expected to provide a guideline to facilitate future developments and applications of the [N–X–N]+ halogen bond.

In this work, we explore the accuracy of post-Hartree–Fock (HF) methods and double-hybrid density functional theory (DFT) for the computation of solid-state NMR chemical shifts. We apply an embedded cluster approach and investigate the convergence with cluster size and embedding for a series of inorganic solids with long-range electrostatic interactions. In a systematic study, we discuss the cluster design, the embedding procedure, and basis set convergence using gauge-including atomic orbital (GIAO) NMR calculations at the DFT and MP2 levels of theory. We demonstrate that the accuracy obtained for the prediction of NMR chemical shifts, which can be achieved for molecular systems, can be carried over to solid systems. An appropriate embedded cluster approach allows one to apply methods beyond standard DFT even for systems for which long-range electrostatic effects are important. We find that an embedded cluster should include at least one sphere of explicit neighbors around the nuclei of interest, given that a sufficiently large point charge and boundary effective potential embedding is applied. Using the pcSseg-3 basis set and GIAOs for the computation of nuclear shielding constants, accuracies of 1.6 ppm for 7Li, 1.5 ppm for 23Na, and 5.1 ppm for 39K as well as 9.3 ppm for 19F, 6.5 ppm for 35Cl, 7.4 ppm for 79Br, and 7.5 ppm for 25Mg as well as 3.8 ppm for 67Zn can be achieved with MP2. Comparing various DFT functionals with HF and MP2, we report the superior quality of results for methods that include post-HF correlation like MP2 and double-hybrid DFT.

We systematically investigate an impact of the size and content of a quantum (QM) region, treated at the density functional theory level, in embedding calculations on one- (OPA) and two-photon absorption (TPA) spectra of the following fluorescent proteins (FPs) models: Aequorea victoria green FP (avGFP) with neutral (avGFP-n) and anionic (avGFP-a) chromophore as well as Citrine FP. We find that amino acid (a.a.) residues as well as water molecules hydrogen-bonded (h-bonded) to the chromophore usually boost both OPA and TPA processes intensity. The presence of hydrophobic a.a. residues in the quantum region also non-negligibly affects both absorption spectra but decreases absorption intensity. We conclude that to reach a quantitative description of OPA and TPA spectra in multiscale modeling of FPs, the quantum region should consist of a chromophore and most of a.a. residues and water molecules in a radius of 0.30-0.35 nm (ca. 200-230 atoms) when the remaining part of the system is approximated by the electrostatic point-charges. The optimal size of the QM region can be reduced to 80-100 atoms by utilizing a more advanced polarizable embedding model. We also find components of the QM region that are specific to a FP under study. We propose that the F165 a.a. residue is important in tuning the TPA spectrum of avGFP-n but not other investigated FPs. In the case of Citrine, Y203 and M69 a.a. residues must definitely be part of the QM subsystem. Furthermore, we find that long-range electrostatic interactions between the QM region and the rest of the protein cannot be neglected even for the most extensive QM regions (ca. 350 atoms).

Orbital-free approaches might offer a way to boost the applicability of density functional theory by orders of magnitude in system size. An important ingredient for this endeavor is the kinetic energy density functional. Snyder et al. [Phys. Rev. Lett.2012, 108, 25300223004593] presented a machine learning approximation for this functional achieving chemical accuracy on a one-dimensional model system. However, a poor performance with respect to the functional derivative, a crucial element in iterative energy minimization procedures, enforced the application of a computationally expensive projection method. In this work we circumvent this issue by including the functional derivative into the training of various machine learning models. Besides kernel ridge regression, the original method of choice, we also test the performance of convolutional neural network techniques borrowed from the field of image recognition.

Following our recent work to reduce a dimension of a set of reference structures along the intrinsic reaction coordinate (IRC) by a classical multidimensional scaling (CMDS) approach (J. Chem. Theory Comput. 2018, 14, 4263-4270), we propose the method to project on-the-fly trajectories into a reduced-dimension subspace determined by the IRC network, using the out-of-sample extension of CMDS. The method was applied to the S(N)2 reaction, OH- + CH3F, in which trajectories show a bifurcating nature around the highly curved region of the IRC path, and to the structural transformation of Au-5 cluster in which the global reaction path network consists of five equilibrium structures and 14 IRCs. It was demonstrated that the present analysis can visualize the dynamics effect by showing a dynamic reaction route on the basis of the static reaction paths.

A fundamental requirement to predict the native conformation, address questions of sequence design and optimization, and gain insights into the folding mechanisms of proteins lies in the definition of an unbiased reaction coordinate that reports on the folding state without the need to compare it to reference values, which might be unavailable for new (designed) sequences. Here, we introduce such a reaction coordinate, which does not depend on previous structural knowledge of the native state but relies solely on the energy partition within the protein: the spectral gap of the pair nonbonded energy matrix (ENergy Gap, ENG). This quantity can be simply calculated along unbiased MD trajectories. We show that upon folding the gap increases significantly, while its fluctuations are reduced to a minimum. This is consistently observed for a diverse set of systems and trajectories. Our approach allows one to promptly identify residues that belong to the folding core as well as residues involved in non-native contacts that need to be disrupted to guide polypeptides to the folded state. The energy gap and fluctuations criteria are then used to develop an automatic detection system which allows us to extract and analyze folding transitions from a generic MD trajectory. We speculate that our method can be used to detect conformational ensembles in dynamic and intrinsically disordered proteins, revealing potential preorganization for binding.

We present a new coarse-grained Cα-based protein model with a nonradial multibody pseudo-improper-dihedral potential that is transferable, time-independent, and suitable for molecular dynamics. It captures the nature of backbone and side-chain interactions between amino acid residues by adapting a simple improper dihedral term for a one-bead-per-residue model. It is parameterized for intrinsically disordered proteins and applicable to simulations of such proteins and their assemblies on millisecond time scales.

The theoretical description of the electronic structure of magnetic insulators and, in particular, of transition-metal oxides (TMOs), MnO, FeO, CoO, NiO, and CuO, poses several problems due to their highly correlated nature. Particularly challenging is the determination of the band gap. The most widely used approach is based on density functional theory (DFT) Kohn-Sham energy levels using self-interaction-corrected functionals (such as hybrid functionals). Here, we present a different approach based on the assumption that the band gap in some TMOs can have a partial Mott-Hubbard character and can be defined as the energy associated with the process Mm+(3dn) + Mm+(3dn) → M(m+1)+(3dn-1) + M(m-1)+(3dn+1). The band gap is thus associated with the removal (ionization potential, I) and addition (electron affinity, A) of one electron to an ion of the lattice. In fact, due to the hybridization of metal with ligand orbitals, these energy contributions are not purely atomic in nature. I and A can be computed accurately using the charge transition level (CTL) scheme. This procedure is based on the calculation of energy levels of charged states and goes beyond the approximations inherent to the Kohn-Sham (KS) approach. The novel and relevant aspect of this work is the extension of CTLs from the domain of point defects to a bulk property such as the band gap. The results show that the calculation based on CTLs provides band gaps in better agreement with experiments than the KS approach, with direct insight into the nature of the gap in these complex systems.

Free-energy calculations based on atomistic Hamiltonians and sampling are key to a first-principles understanding of biomolecular processes, material properties, and macromolecular chemistry. Here, we generalize the free-energy perturbation method and derive nonlinear Hamiltonian transformation sequences yielding free-energy estimates with minimal mean squared error with respect to the exact values. Our variational approach applies to finite sampling and holds for any finite number of intermediate states. We show that our sequences are also optimal for the Bennett acceptance ratio (BAR) method, thereby generalizing BAR to small sampling sizes and non-Gaussian error distributions.

Most classical nonpolarizable ion potential models underestimate the solubility values of NaCl and KCl in water significantly. We determine Lennard-Jones parameters of Na+, K+, and Cl- that reproduce the solubility as well as the hydration free energy in dilute aqueous solutions for three water potential models, SPC/E, TIP3P, and TIP4P/2005. The ion-oxygen distance in the solution and the cation-anion distance in salt are also considered in the parametrization. In addition to the target properties, the hydration enthalpy, hydration entropy, self-diffusion coefficient, coordination number, lattice energy, enthalpy of solution, density, viscosity, and number of contact ion pairs are calculated for comparison with 17 frequently used or recently developed ion potential models. The overall performance of each ion model is represented by a global score using a scheme that was originally developed for comparison of water potential models. The global score is better for our models than for the other 17 models not only because of the quite good prediction for the solubility but also because of the relatively small deviation from the experimental value for many of the other properties.

Electron and charge transfers are part of many vital processes in nature and technology. Ab initio descriptions of these processes provide useful insights that can be utilized for applications. A combination of the embedded cluster material model and nonorthogonal configuration interaction (NOCI), in which the cluster wave functions are expanded in many-electron basis functions (MEBFs) consisting of spin-adapted, antisymmetrized products of multiconfigurational wave functions of fragments (which are usually molecules) in the cluster, appears to provide a compromise between accuracy and calculation time. Additional advantages of this NOCI-Fragments approach are the chemically convenient interpretation of the wave function in terms of molecular states, and the direct accessibility of electronic coupling between diabatic states to describe energy and electron transfer processes. Bottlenecks in this method are the large number of two-electron integrals that have to be handled for the calculation of an electronic coupling matrix element and the enormous number of matrix elements over determinant pairs that have to be evaluated for the calculation of one matrix element between the MEBFs. We show here how we created a reduced common molecular orbital basis that is utilized to significantly reduce the number of two-electron integrals that need to be handled. The results obtained with this basis do not show any loss of accuracy in relevant quantities like electronic couplings and vertical excitation energies. We also show a significant reduction in computation time without loss in accuracy when matrix elements over determinant pairs with small weights are neglected in the NOCI. These improvements in the methodology render NOCI-Fragments to be also applicable to treat clusters of larger molecular systems with larger atomic basis sets and larger active spaces, as the computation time becomes dependent on the number of occupied orbitals and less dependent on the size of the active space.

The system-size dependence of computed mutual diffusion coefficients of multicomponent mixtures is investigated, and a generalized correction term is derived. The generalized finite-size correction term was validated for the ternary molecular mixture chloroform/acetone/methanol as well as 28 ternary LJ systems. It is shown that only the diagonal elements of the Fick matrix show system-size dependency. The finite-size effects of these elements can be corrected by adding the term derived by Yeh and Hummer (J. Phys. Chem. B2004, 108, 15873–15879). By performing an eigenvalue analysis of the finite-size effects of the matrix of Fick diffusivities we show that the eigenvector matrix of Fick diffusivities does not depend on the size of the simulation box. Only eigenvalues, which describe the speed of diffusion, depend on the size of the system. An analytic relation for finite-size effects of the matrix of Maxwell–Stefan diffusivities was developed. All Maxwell–Stefan diffusivities depend on the system size, and the required correction depends on the matrix of thermodynamic factors., Journal of Chemical Theory and Computation, 16 (6), ISSN:1549-9618, ISSN:1549-9626

Calculating free energies of binding (ΔGbind) between ligands and their target protein is of major interest to drug discovery and safety, yet it is still associated with several challenges and difficulties. Linear interaction energy (LIE) is an efficient in silico method for ΔGbind computation. LIE models can be trained and used to directly calculate binding affinities from interaction energies involving ligands in the bound and unbound states only, and LIE can be combined with statistical weighting to calculate ΔGbind for flexible proteins that may bind their ligands in multiple orientations. Here, we investigate if LIE predictions can be effectively improved by explicitly including the entropy of (de)solvation into our free-energy calculations. For that purpose, we combine LIE calculations for the protein-ligand-bound state with explicit free-energy perturbation to rigorously compute the unbound ligand's solvation free energy. We show that for 28 Cytochrome P450 2A6 (CYP2A6) ligands, coupling LIE with alchemical solvation free-energy calculation helps to improve obtained correlation between computed and reference (experimental) binding data.

The prototype Diels-Alder (DA) reaction between butadiene and ethene (system 1) and the DA reaction involving 1-methoxy-butadiene and cyano-ethylene (system 2) are investigated with an explicit-time-dependent Density Functional Theory approach. Bond orders and atomic net charges obtained in the dynamics at the transition state geometry and along the reaction coordinate toward reactants are used to provide a picture of the process in terms of VB/Lewis resonance structures that contribute to a resonance hybrid. The entire dynamics can be divided into different domains (reactant-like, product-like, and transition state domains) where different Lewis resonance structures contribute with different weights. The relative importance of these three domains varies along the reaction coordinate. In addition to the usual reactant-like and product-like covalent Lewis structures, ionic Lewis structures have non-negligible weights. In system 2, the electron-donor OCH3 on the diene and the electron-acceptor CN on the dienophile make more important the contributions of ionic Lewis structures that stabilize the transition state and determine the decrease of the reaction barrier with respect to system 1.

Face-sharing octahedral dinuclear Cr(III) compounds with d3-d3 electronic configurations represent nontrivial examples of electronic complexity, posing particular challenges for theoretical and computational studies. A tris-hydroxy-bridged Cr(III)-Cr(III) system has proven to be a richly rewarding target for studies of magnetism and electron paramagnetic resonance spectroscopy. It was also reported to be a peculiarly difficult system to treat with density functional theory (DFT). In this work the magnetic coupling problem for this dimer is approached with broken-symmetry (BS)-DFT and multireference calculations that utilize the density matrix renormalization group (DMRG) to handle full-valence active spaces. BS-DFT is shown to recover the correct ordering and energy spacing of Heisenberg spin states if used in conjunction with appropriate spin projection procedures, albeit with pronounced functional sensitivity. The contrasting conclusions of previous studies are traced to incorrect inclusion of electronically excited configurations. Analysis of the direct and differential overlap of corresponding orbital pairs from the BS-DFT solution indicates that metal-metal through-space interaction is the dominant contributor to antiferromagnetic coupling. At the DFT level a procedure that utilizes pseudopotential substitution is demonstrated that allows evaluation of the direct exchange vs superexchange contributions. A complementary description is obtained with DMRG-SCF calculations that enable state-averaged CASSCF calculations with both metal and bridge orbitals in the active space. A localized orbital subspace analysis supports the DFT conclusions that in contrast to doubly bridged isoelectronic analogues, antiferromagnetic coupling in the chromium dimer arises primarily from direct metal-metal interaction but is significantly enhanced by ligand-mediated superexchange.

Knowledge on thermodynamic and transport properties of aqueous solutions of carbohydrates is of great interest for process and product design in the food, pharmaceutical, and biotechnological industries. Molecular simulation is a powerful tool to calculate these properties, but current classical force fields cannot provide accurate estimates for all properties of interest. The poor performance of the force fields is mainly observed for concentrated solutions, where solute-solute interactions are overestimated. In this study, we propose a method to refine force fields, such that solute-solute interactions are more accurately described. The OPLS force field combined with the SPC/Fw water model is used as a basis. We scale the nonbonded interaction parameters of sucrose, a disaccharide. The scaling factors are chosen in such a way that experimental thermodynamic and transport properties of aqueous solutions of sucrose are accurately reproduced. Using a scaling factor of 0.8 for Lennard-Jones energy parameters (ϵ) and a scaling factor of 0.95 for partial atomic charges (q), we find excellent agreement between experiments and computed liquid densities, thermodynamic factors, shear viscosities, self-diffusion coefficients, and Fick (mutual) diffusion coefficients. The transferability of these optimum scaling factors to other carbohydrates is verified by computing thermodynamic and transport properties of aqueous solutions of d-glucose, a monosaccharide. The good agreement between computed properties and experiments suggests that the scaled interaction parameters are transferable to other carbohydrates, especially for concentrated solutions.

Warfarin, a widely used oral anticoagulant, is prescribed as a racemic mixture. Each enantiomer of neutral Warfarin can exist in 20 possible tautomeric states leading to complex pharmacokinetics and uncertainty as to the relevant species under different conditions. Here, the ability of alternative computational approaches to predict the preferred tautomeric form(s) of neutral Warfarin in different solvents is examined. It is shown that varying the method used to estimate the heat of formation in vacuum (direct or via homodesmic reactions), whether entropic corrections were included, and the method used to estimate the free enthalpy of solvation (i.e., PCM, COSMO, or SMD implicit models or explicit solvent) lead to large differences in the predicted rank and relative populations of the tautomers. In this case, only a combination of the enthalpy of formation using homodesmic reactions and explicit solvent to estimate the free enthalpy of solvation yielded results compatible with the available experimental data. The work also suggests that a small but significant subset of the possible Warfarin tautomers are likely to be physiologically relevant.

Nucleophilic addition onto a carbonyl moiety is strongly affected by solvent, and correctly simulating this solvent effect is often beyond the capability of single-scale quantum mechanical (QM) models. This work explores multiscale approaches for the description of the reversible and highly solvent-sensitive nucleophilic N|···C=O bond formation in an Me2N–(CH2)3–CH=O molecule. In the first stage of this work, we rigorously compare and test four recent quantum mechanical/molecular mechanical (QM/MM) explicit solvation models, employing a QM description of water molecules in spherical regions around both the oxygen and the nitrogen atom of the solute. The accuracy of the models is benchmarked against a reference QM simulation, focusing on properties of the solvated Me2N–(CH2)3–CH=O molecule in its ring-closed form. In the second stage, we select one of the models (continuous adaptive QM/MM) and use it to obtain a reliable free energy profile for the N|···C bond formation reaction. We find that the dual-sphere approach allows the model to accurately account for solvent reorganization along the entire reaction path. In contrast, a simple microsolvation model cannot adapt to the changing conditions and provides an incorrect description of the reaction process.

Computational prediction of noncovalent binding free energies with methods based on molecular mechanical force fields has become increasingly routine in drug discovery projects, where they promise to speed the discovery of small molecule ligands to bind targeted proteins with high affinity. Because the reliability of free energy methods still has significant room for improvement, new force fields, or modifications of existing ones, are regularly introduced with the aim of improving the accuracy of molecular simulations. However, comparatively little work has been done to systematically assess how well force fields perform, particularly in relation to the calculation of binding affinities. Hardware advances have made these calculations feasible, but comprehensive force field assessments for protein-ligand sized systems still remain costly. Here, we turn to cyclodextrin host-guest systems, which feature many hallmarks of protein-ligand binding interactions but are generally much more tractable due to their small size. We present absolute binding free energy and enthalpy calculations, using the attach-pull-release (APR) approach, on a set of 43 cyclodextrin-guest pairs for which experimental ITC data are available. The test set comprises both α- and β-cyclodextrin hosts binding a series of small organic guests, each with one of three functional groups: ammonium, alcohol, or carboxylate. Four water models are considered (TIP3P, TIP4Pew, SPC/E, and OPC), along with two partial charge assignment procedures (RESP and AM1-BCC) and two cyclodextrin host force fields. The results suggest a complex set of considerations when choosing a force field for biomolecular simulations. For example, some force field combinations clearly outperform others at the binding enthalpy calculations but not for the binding free energy. Additionally, a force field combination which we expected to be the worst performer gave the most accurate binding free energies - but the least accurate binding enthalpies. The results have implications for the development of improved force fields, and we propose this test set, and potential future elaborations of it, as a powerful validation suite to evaluate new force fields and help guide future force field development.

The nonadditive many-body interactions are significant for structural and thermodynamic properties of condensed phase systems. In this work we examined the many-body interaction energy of a large number of common organic/biochemical molecular clusters, which consist of 18 chemical species and cover nine common organic elements, using the Møller-Plesset perturbation theory to the second order (MP2) [ Møller et al. Phys. Rev. 1934 , 46 , 618 . ]. We evaluated the capability of Thole-based dipole induction models to capture the many-body interaction energy. Three models were compared: the original model and parameters used by the AMOEBA force field, a variation of this original model where the damping parameters have been reoptimized to MP2 data, and a third model where the damping function form applied to the permanent electric field is modified. Overall, we find the simple classical atomic dipole models are able to capture the 3- and 4-body interaction energy across a wide variety of organic molecules in various intermolecular configurations. With modified Thole models, it is possible to further improve the agreement with MP2 results. These models were also tested on systems containing metal/halogen ions to examine the accuracy and transferability. This work suggests that the form of damping function applied to the permanent electrostatic field strongly affects the distance dependence of polarization energy at short intermolecular separations.

The protein mediated hydrolysis of nucleoside triphosphates such as ATP or GTP is one of the most important and challenging biochemical reactions in nature. The chemical environment (water structure, catalytic metal, and amino acid residues) adjacent to the hydrolysis site contains hundreds of atoms, usually greatly limiting the amount of the free energy sampling that one can achieve from computationally demanding electronic structure calculations such as QM/MM simulations. Therefore, the combination of QM/MM molecular dynamics with the recently developed transition-tempered metadynamics (TTMetaD), an enhanced sampling method that can provide a high-quality free energy estimate at an early stage in a simulation, is an ideal approach to address the biomolecular nucleoside triphosphate hydrolysis problem. In this work the ATP hydrolysis process in monomeric and filamentous actin is studied as an example application of the combined methodology. The performance of TTMetaD in these demanding QM/MM simulations is compared with that of the more conventional well-tempered metadynamics (WTMetaD). Our results show that TTMetaD exhibits much better exploration of the hydrolysis reaction free energy surface in two key collective variables (CVs) during the early stages of the QM/MM simulation than does WTMetaD. The TTMetaD simulations also reveal that a key third degree of freedom, the O-H bond-breaking and proton transfer from the lytic water, must be biased for TTMetaD to converge fully. To perturb the NTP hydrolysis dynamics to the least extent and to properly focus the MetaD free energy sampling, we also adopt here the recently developed metabasin metadynamics (MBMetaD) to construct a self-limiting bias potential that only applies to the lytic water after its nucleophilic attack of the phosphate of ATP. With these new, state-of-the-art enhanced sampling metadynamics techniques, we present an effective and accurate computational strategy for combining QM/MM molecular dynamics simulation with free energy sampling methodology, including a means to analyze the convergence of the calculations through robust numerical criteria.

Aromatic molecules with π electrons are commonly involved in chemical and biological recognitions. For example, nucleobases play central roles in DNA/RNA structure and their interactions with proteins. The delocalization of the π electrons is responsible for the high polarizability of aromatic molecules. In this work, the AMOEBA force field has been developed and applied to 5 regular nucleobases and 12 aromatic molecules. The permanent electrostatic energy is expressed as atomic multipole interactions between atom pairs, and many-body polarization is accounted for by mutually induced atomic dipoles. We have systematically investigated aromatic ring stacking and aromatic-water interactions for nucleobases and aromatic molecules, as well as base-base hydrogen-bonding pair interactions, all at various distances and orientations. van der Waals parameters were determined by comparison to the quantum mechanical interaction energy of these dimers and fine-tuned using condensed phase simulation. By comparing to quantum mechanical calculations, we show that the resulting classical potential is able to accurately describe molecular polarizability, molecular vibrational frequency, and dimer interaction energy of these aromatic systems. Condensed phase properties, including hydration free energy, liquid density, and heat of vaporization, are also in good overall agreement with experimental values. The structures of benzene liquid phase and benzene-water solution were also investigated by simulation and compared with experimental and PDB structure derived statistical results.

We describe a general scheme to obtain force-field parameters for classical molecular dynamics simulations of conjugated polymers. We identify a computationally inexpensive methodology for calculation of accurate intermonomer dihedral potentials and partial charges. Our findings indicate that the use of a two-step methodology of geometry optimization and single-point energy calculations using DFT methods produces potentials which compare favorably to high level theory calculation. We also report the effects of varying the conjugated backbone length and alkyl side-chain lengths on the dihedral profiles and partial charge distributions and determine the existence of converged lengths above which convergence is achieved in the force-field parameter sets. We thus determine which calculations are required for accurate parametrization and the scope of a given parameter set for variations to a given molecule. We perform simulations of long oligomers of dioctylfluorene and hexylthiophene in explicit solvent and find peristence lengths and end-length distributions consistent with experimental values.

We present the AMBER ff15ipq force field for proteins, the second-generation force field developed using the Implicitly Polarized Q (IPolQ) scheme for deriving implicitly polarized atomic charges in the presence of explicit solvent. The ff15ipq force field is a complete rederivation including more than 300 unique atomic charges, 900 unique torsion terms, 60 new angle parameters, and new atomic radii for polar hydrogens. The atomic charges were derived in the context of the SPC/Eb water model, which yields more-accurate rotational diffusion of proteins and enables direct calculation of nuclear magnetic resonance (NMR) relaxation parameters from molecular dynamics simulations. The atomic radii improve the accuracy of modeling salt bridge interactions relative to contemporary fixed-charge force fields, rectifying a limitation of ff14ipq that resulted from its use of pair-specific Lennard-Jones radii. In addition, ff15ipq reproduces penta-alanine J-coupling constants exceptionally well, gives reasonable agreement with NMR relaxation rates, and maintains the expected conformational propensities of structured proteins/peptides, as well as disordered peptides-all on the microsecond (μs) time scale, which is a critical regime for drug design applications. These encouraging results demonstrate the power and robustness of our automated methods for deriving new force fields. All parameters described here and the mdgx program used to fit them are included in the AmberTools16 distribution.

We compare the performance of five magnesium (Mg(2+)) ion models in simulations of an RNA stem loop which has an experimentally determined divalent ion dependent conformational shift. We show that despite their differences in parametrization and resulting van der Waals terms, including differences in the functional form of the nonbonded potential, when the RNA adopts its folded conformation, all models behave similarly across ten independent microsecond length simulations with each ion model. However, when the entire structure ensemble is accounted for, chelation of Mg(2+) to RNA is seen in three of the five models, most egregiously and likely artifactual in simulations using a 12-6-4 model for the Lennard-Jones potential. Despite the simple nature of the fixed point-charge and van der Waals sphere models employed, and with the exception of the likely oversampled directed chelation of the 12-6-4 potential models, RNA-Mg(2+) interactions via first shell water molecules are surprisingly well described by modern parameters, allowing us to observe the spontaneous conformational shift from Mg(2+) free RNA to Mg(2+) associated RNA structure in unrestrained molecular dynamics simulations.

Because standard molecular dynamics (MD) simulations are unable to access time scales of interest in complex biomolecular systems, it is common to “stitch together” information from multiple shorter trajectories using approximate Markov state model (MSM) analysis. However, MSMs may require significant tuning and can yield biased results. Here, by analyzing some of the longest protein MD data sets available (>100 μs per protein), we show that estimators constructed based on exact non-Markovian (NM) principles can yield significantly improved mean first-passage times (MFPTs) for protein folding and unfolding. In some cases, MSM bias of more than an order of magnitude can be corrected when identical trajectory data are reanalyzed by non-Markovian approaches. The NM analysis includes “history” information, higher order time correlations compared to MSMs, that is available in every MD trajectory. The NM strategy is insensitive to fine details of the states used and works well when a fine time-discretization (i.e., small “lag time”) is used.

Proper treatment of nonbonded interactions is essential for the accuracy of molecular dynamics (MD) simulations, especially in studies of lipid bilayers. The use of the CHARMM36 force field (C36 FF) in different MD simulation programs can result in disagreements with published simulations performed with CHARMM due to differences in the protocols used to treat the long-range and 1-4 nonbonded interactions. In this study, we systematically test the use of the C36 lipid FF in NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM. A wide range of Lennard-Jones (LJ) cutoff schemes and integrator algorithms were tested to find the optimal simulation protocol to best match bilayer properties of six lipids with varying acyl chain saturation and head groups. MD simulations of a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer were used to obtain the optimal protocol for each program. MD simulations with all programs were found to reasonably match the DPPC bilayer properties (surface area per lipid, chain order parameters, and area compressibility modulus) obtained using the standard protocol used in CHARMM as well as from experiments. The optimal simulation protocol was then applied to the other five lipid simulations and resulted in excellent agreement between results from most simulation programs as well as with experimental data. AMBER compared least favorably with the expected membrane properties, which appears to be due to its use of the hard-truncation in the LJ potential versus a force-based switching function used to smooth the LJ potential as it approaches the cutoff distance. The optimal simulation protocol for each program has been implemented in CHARMM-GUI. This protocol is expected to be applicable to the remainder of the additive C36 FF including the proteins, nucleic acids, carbohydrates, and small molecules.

Despite decades of investigations, the principal mechanisms responsible for the high affinity and specificity of proteins for key physiological cations K(+), Na(+), and Ca(2+) remain a hotly debated topic. At the core of the debate is an apparent need (or lack thereof) for an accurate description of the electrostatic response of the charge distribution in a protein to the binding of an ion. These effects range from partial electronic polarization of the directly ligating atoms to long-range effects related to partial charge transfer and electronic delocalization effects. While accurate modeling of cation recognition by metalloproteins warrants the use of quantum-mechanics (QM) calculations, the most popular approximations used in major biomolecular simulation packages rely on the implicit modeling of electronic polarization effects. That is, high-level QM computations for ion binding to proteins are desirable, but they are often unfeasible, because of the large size of the reactive-site models and the need to sample conformational space exhaustively at finite temperature. Several solutions to this challenge have been proposed in the field, ranging from the recently developed Drude polarizable force-field for simulations of metalloproteins to approximate tight-binding density functional theory (DFTB). To delineate the usefulness of different approximations, we examined the accuracy of three recent and commonly used theoretical models and numerical algorithms, namely, CHARMM C36, the latest developed Drude polarizable force fields, and DFTB3 with the latest 3OB parameters. We performed MD simulations for 30 cation-selective proteins with high-resolution X-ray structures to create ensembles of structures for analysis with different levels of theory, e.g., additive and polarizable force fields, DFTB3, and DFT. The results from DFT computations were used to benchmark CHARMM C36, Drude, and DFTB3 performance. The explicit modeling of quantum effects unveils the key electrostatic properties of the protein sites and the importance of specific ion-protein interactions. One of the most interesting findings is that secondary coordination shells of proteins are noticeably perturbed in a cation-dependent manner, showing significant delocalization and long-range effects of charge transfer and polarization upon binding Ca(2+).

This article reports a novel theoretical approach related to the analysis of molecular excited states. The strategy introduced here involves gathering two pieces of physical information, coming from Hilbert and direct space operations, into a general, unique quantum mechanical descriptor of electronic transitions locality. Moreover, the projection of Hilbert and direct space-derived indices in an Argand plane delivers a straightforward way to visually probe the ability of a dye to undergo a long- or short-range charge-transfer. This information can be applied, for instance, to the analysis of the electronic response of families of dyes to light absorption by unveiling the trend of a given push-pull chromophore to increase the electronic cloud polarization magnitude of its main transition with respect to the size extension of its conjugated spacer. We finally demonstrate that all the quantities reported in this article can be reliably approximated by a linear algebraic derivation, based on the contraction of detachment/attachment density matrices from canonical to atomic space. This alternative derivation has the remarkable advantage of a very low computational cost with respect to the previously used numerical integrations, making fast and accurate characterization of large molecular systems excited states easily affordable.

We used microsecond time scale molecular dynamics simulations to compute, at high precision, binding enthalpies for cucurbit[7]uril (CB7) with eight guests in aqueous solution. The results correlate well with experimental data from previously published isothermal titration calorimetry studies, and decomposition of the computed binding enthalpies by interaction type provides plausible mechanistic insights. Thus, dispersion interactions appear to play a key role in stabilizing these complexes, due at least in part to the fact that their packing density is greater than that of water. On the other hand, strongly favorable Coulombic interactions between the host and guests are compensated by unfavorable solvent contributions, leaving relatively modest electrostatic contributions to the binding enthalpies. The better steric fit of the aliphatic guests into the circular host appears to explain why their binding enthalpies tend to be more favorable than those of the more planar aromatic guests. The present calculations also bear on the validity of the simulation force field. Somewhat unexpectedly, the TIP3P water yields better agreement with experiment than the TIP4P-Ew water model, although the latter is known to replicate the properties of pure water more accurately. More broadly, the present results demonstrate the potential for computational calorimetry to provide atomistic explanations for thermodynamic observations.

The epidermal growth factor receptor (EGFR) is aberrantly activated in various cancer cells and an important target for cancer treatment. Deep understanding of EGFR conformational changes between the active and inactive states is of pharmaceutical interest. Here we present a strategy combining multiply targeted molecular dynamics simulations, unbiased molecular dynamics simulations, and Bayesian clustering to investigate transition pathways during the activation/inactivation process of EGFR kinase domain. Two distinct pathways between the active and inactive forms are designed, explored, and compared. Based on Bayesian clustering and rough two-dimensional free energy surfaces, the energy-favorable pathway is recognized, though DFG-flip happens in both pathways. In addition, another pathway with different intermediate states appears in our simulations. Comparison of distinct pathways also indicates that disruption of the Lys745-Glu762 interaction is critically important in DFG-flip while movement of the A-loop significantly facilitates the conformational change. Our simulations yield new insights into EGFR conformational transitions. Moreover, our results verify that this approach is valid and efficient in sampling of protein conformational changes and comparison of distinct pathways.

Owing to recent developments in computational algorithms and architectures, it is now computationally tractable to explore biologically relevant, equilibrium dynamics of realistically sized functional proteins using all-atom molecular dynamics simulations. Molecular dynamics simulations coupled with Markov state models is a nascent but rapidly growing technology that is enabling robust exploration of equilibrium dynamics. The objective of this work is to explore the challenges of coupling molecular dynamics simulations and Markov state models in the study of functional proteins. Using recent studies as a framework, we explore progress in sampling, model building, model selection, and coarse-grained analysis of models. Our goal is to highlight some of the current challenges in applying Markov state models to realistically sized proteins and spur discussion on advances in the field.

Continuing our effort to introduce d-amino-acid residues in the united residue (UNRES) force field developed in our laboratory, in this work the Cα ··· Cα ··· Cα backbone-virtual-bond-valence-angle (θ) potentials for systems containing d-amino-acid residues have been developed. The potentials were determined by integrating the combined energy surfaces of all possible triplets of terminally blocked glycine, alanine, and proline obtained with ab initio molecular quantum mechanics at the MP2/6-31G(d,p) level to calculate the corresponding potentials of mean force (PMFs). Subsequently, analytical expressions were fitted to the PMFs to give the virtual-bond-valence potentials to be used in UNRES. Alanine represented all types of amino-acid residues except glycine and proline. The blocking groups were either the N-acetyl and N′,N′-dimethyl or N-acetyl and pyrrolidyl group, depending on whether the residue next in sequence was an alanine-type or a proline residue. A total of 126 potentials (63 symmetry-unrelated potentials for each set of terminally blocking groups) were determined. Together with the torsional, double-torsional, and side-chain-rotamer potentials for polypeptide chains containing d-amino-acid residues determined in our earlier work (Sieradzan et al. J. Chem. Theory Comput., 2012, 8, 4746), the new virtual-bond-angle (θ) potentials now constitute the complete set of physics-based potentials with which to run coarse-grained simulations of systems containing d-amino-acid residues. The ability of the extended UNRES force field to reproduce thermodynamics of polypeptide systems with d-amino-acid residues was tested by comparing the experimentally measured and the calculated free energies of helix formation of model KLALKLALxxLKLALKLA peptides, where x denotes any d- or l- amino-acid residue. The obtained results demonstrate that the UNRES force field with the new potentials reproduce the changes of free energies of helix formation upon d-substitution but overestimate the free energies of helix formation. To test the ability of UNRES with the new potentials to reproduce the structures of polypeptides with d-amino-acid residues, an ab initio replica-exchange folding simulation of thurincin H from Bacillus thuringiensis, which has d-amino-acid residues in the sequence, was carried out. UNRES was able to locate the native α-helical hairpin structure as the dominant structure even though no native sulfide–carbon bonds were present in the simulation.

Broken-symmetry density functional theory (BS-DFT) calculations are assessed for redox energetics [Cu(SCH3)2]1–/0, [Cu(NCS)2]1–/0, [FeCl4]1–/0, and [Fe(SCH3)4]1–/0 against vertical detachment energies (VDE) from valence photoelectron spectroscopy (PES), as a prelude to studies of metalloprotein analogs. The M06 and B3LYP hybrid functionals give VDE that agree with the PES VDE for the Fe complexes, but both underestimate it by ∼400 meV for the Cu complexes; other hybrid functionals give VDEs that are an increasing function of the amount of Hartree–Fock (HF) exchange and so cannot show good agreement for both Cu and Fe complexes. Range-separated (RS) functionals appear to give a better distribution of HF exchange since the negative HOMO energy is approximately equal to the VDEs but also give VDEs dependent on the amount of HF exchange, sometimes leading to ground states with incorrect electron configurations; the LRC-ωPBEh functional reduced to 10% HF exchange at short-range give somewhat better values for both, although still ∼150 meV too low for the Cu complexes and ∼50 meV too high for the Fe complexes. Overall, the results indicate that while HF exchange compensates for self-interaction error in DFT calculations of both Cu and Fe complexes, too much may lead to more sensitivity to nondynamical correlation in the spin-polarized Fe complexes.

Molecular mechanics with all-atom models was used to understand the conformational preference of tandem guanine-adenine (GA) noncanonical pairs in RNA. These tandem GA pairs play important roles in determining stability, flexibility, and structural dynamics of RNA tertiary structures. Previous solution structures showed that these tandem GA pairs adopt either imino (cis Watson–Crick/Watson–Crick A-G) or sheared (trans Hoogsteen/sugar edge A-G) conformations depending on the sequence and orientation of the adjacent closing base pairs. The solution structures (GCGGACGC)2 [Biochemistry, 1996, 35, 9677–9689] and (GCGGAUGC)2 [Biochemistry, 2007, 46, 1511–1522] demonstrate imino and sheared conformations for the two central GA pairs, respectively. These systems were studied using molecular dynamics and free energy change calculations for conformational changes, using umbrella sampling. For the structures to maintain their native conformations during molecular dynamics simulations, a modification to the standard Amber ff10 force field was required, which allowed the amino group of guanine to leave the plane of the base [J. Chem. Theory Comput., 2009, 5, 2088–2100] and form out-of-plane hydrogen bonds with a cross-strand cytosine or uracil. The requirement for this modification suggests the importance of out-of-plane hydrogen bonds in stabilizing the native structures. Free energy change calculations for each sequence demonstrated the correct conformational preference when the force field modification was used, but the extent of the preference is underestimated.

Density functional theory (DFT) provides a formally exact framework for quantum embedding. The appearance of nonadditive kinetic energy contributions in this context poses significant challenges, but using optimized effective potential (OEP) methods, various groups have devised DFT-in-DFT methods that are equivalent to Kohn–Sham (KS) theory on the whole system. This being the case, we note that a very considerable simplification arises from doing KS theory instead. We then describe embedding schemes that enforce Pauli exclusion via a projection technique, completely avoiding numerically demanding OEP calculations. Illustrative applications are presented using DFT-in-DFT, wave-function-in-DFT, and wave-function-in-Hartree–Fock embedding, and using an embedded many-body expansion.

The design of ligands with high affinity and specificity remains a fundamental challenge in understanding molecular recognition and developing therapeutic interventions. Charge optimization theory addresses this problem by determining ligand charge distributions that produce the most favorable electrostatic contribution to the binding free energy. The theory has been applied to the design of binding specificity as well. However, the formulations described only treat a rigid ligand—one that does not change conformation upon binding. Here, we extend the theory to treat induced-fit ligands for which the unbound ligand conformation may differ from the bound conformation. We develop a thermodynamic pathway analysis for binding contributions relevant to the theory, and we illustrate application of the theory using HIV-1 protease with our previously designed and validated subnanomolar inhibitor. Direct application of rigid charge optimization approaches to nonrigid cases leads to very favorable intramolecular electrostatic interactions that are physically unreasonable, and analysis shows the ligand charge distribution massively stabilizes the preconformed (bound) conformation over the unbound. After analyzing this case, we provide a treatment for the induced-fit ligand charge optimization problem that produces physically realistic results. The key factor is introducing the constraint that the free energy of the unbound ligand conformation be lower or equal to that of the preconformed ligand structure, which corresponds to the notion that the unbound structure is the ground unbound state. Results not only demonstrate the applicability of this methodology to discovering optimized charge distributions in an induced-fit model, but also provide some insights into the energetic consequences of ligand conformational change on binding. Specifically, the results show that, from an electrostatic perspective, induced-fit binding is not an adaptation designed to enhance binding affinity; at best, it can only achieve the same affinity as optimized rigid binding.

We present an extension of the coarse-grained MARTINI model for proteins and apply this extension to amyloid- and elastin-like peptides. Atomistic simulations of tetrapeptides, octapeptides, and longer peptides in solution are used as a reference to parametrize a set of pseudodihedral potentials that describe the internal flexibility of MARTINI peptides. We assess the performance of the resulting model in reproducing various structural properties computed from atomistic trajectories of peptides in water. The addition of new dihedral angle potentials improves agreement with the contact maps computed from atomistic simulations significantly. We also address the question of which parameters derived from atomistic trajectories are transferable between different lengths of peptides. The modified coarse-grained model shows reasonable transferability of parameters for the amyloid- and elastin-like peptides. In addition, the improved coarse-grained model is also applied to investigate the self-assembly of β-sheet forming peptides on the microsecond time scale. The octapeptides SNNFGAIL and (GV)(4) are used to examine peptide aggregation in different environments, in water, and at the water-octane interface. At the interface, peptide adsorption occurs rapidly, and peptides spontaneously aggregate in favor of stretched conformers resembling β-strands.

Prolyl oligopeptidase (POP), a member of the prolyl endopeptidase family, is known to play a role in several neurological disorders. Its primary function is to cleave a wide range of small oligopeptides, including neuroactive peptides. We have used force biased molecular dynamics simulation to study the binding mechanism of POP. We examined three possible binding pathways using Steered Molecular Dynamics (SMD) and Umbrella Sampling (US) on a crystal structure of porcine POP with bound Z-pro-prolinal (ZPP). Using SMD, an exit pathway between the first and seventh blade of the β-propeller domain of POP was found to be a nonviable route. US on binding pathways through the β-propeller tunnel and the TYR190-GLN208 flexible loop at the interface between both POP domains allowed us to isolate the flexible loop pathway as the most probable. Further analysis of that pathway suggests a long-range covariation of the interdomain H-bond network, which indicates the possibility of large-scale domain reorientation observed in bacterial homologues and hypothesized to also occur in human POP.

Accelerated molecular dynamics (aMD) has been shown to enhance conformational space sampling relative to classical molecular dynamics; however, the exponential reweighting of aMD trajectories, which is necessary for the calculation of free energies relating to the classical system, is oftentimes problematic, especially for systems larger than small poly peptides. Here, we propose a method of accelerating only the degrees of freedom most pertinent to sampling, thereby reducing the total acceleration added to the system and improving the convergence of calculated ensemble averages, which we term selective aMD. Its application is highlighted in two biomolecular cases. First, the model system alanine dipeptide is simulated with classical MD, all-dihedral aMD, and selective aMD, and these results are compared to the infinite sampling limit as calculated with metadynamics. We show that both forms of aMD enhance the convergence of the underlying free energy landscape by 5-fold relative to classical MD; however, selective aMD can produce improved statistics over all-dihedral aMD due to the improved reweighting. Then we focus on the pharmaceutically relevant case of computing the free energy of the decoupling of oseltamivir in the active site of neuraminidase. Results show that selective aMD greatly reduces the cost of this alchemical free energy transformation, whereas all-dihedral aMD produces unreliable free energy estimates.

The convergence properties of the absolute single-molecule configurational entropy and the correction terms used to estimate it are investigated using microsecond molecular dynamics simulation of a peptide test system and an improved methodology. The results are compared with previous applications for systems of diverse chemical nature. It is shown that (i) the effect of anharmonicity is small, (ii) the effect of pairwise correlation is typically large, and (iii) the latter affects to a larger extent the entropy estimate of thermodynamic states characterized by a higher motional correlation. The causes of such deviations from a quasi-harmonic behavior are explained. This improved approach provides entropies also for molecular systems undergoing conformational transitions and characterized by highly frustrated energy surfaces, thus not limited to systems sampling a single quasi-harmonic basin. Overall, this study emphasizes the need for extensive phase-space sampling in order to obtain a reliable estimation of entropic contributions.