Journal of Chemical Theory and Computation最新文献

Accurate computation of protein-ligand binding free energy remains an elusive goal due to inherent difficulties involved in the accurate calculation of gas-phase protein-ligand interaction energy, the entropy, and the solvation energy. In this study, we explore the use of fragment quantum chemical calculations for improved accuracy in protein-ligand binding free energy calculations. The present work demonstrated that the gas-phase protein-ligand interaction energies can be accurately calculated by the molecular fractionation with conjugate caps method as verified by comparison with the full quantum calculations for several protein-ligand systems. The m06-2x/6-31+G* level of density functional theory calculation with basis set superposition error correction is found to give excellent protein-ligand interaction energies. The quantum calculated protein-ligand interaction energies are then combined with implicit solvation methods to obtain absolute binding free energies and the results are shown to be sensitive to the specific solvation models used. In particular, the accuracy of the quantum calculated binding free energies is significantly improved over that of the force field calculations using the same solvation models in terms of mean absolute errors. However, the correlation coefficients with respect to the experimental data do not show improvement over the corresponding result computed from the force field. Such result and the related analysis underscore the critical importance of solvation energies to the binding free energies and the need for developing new methods to calculate solvation energies more accurately in the future.

Low-dimensional water, despite the relative simplicity of its constituents, exhibits a vast range of phenomena that are of central importance in natural sciences. A large number of bulk as well as nanoscale polymorphs offer engineering possibilities for technological applications such as desalinization, drug delivery, or biological interfacing. However, little is known about the stability of such structures. Therefore, in this study, we employ an array of state-of-the-art computational techniques to study the vibrational properties of ice Ih and XI in their bulk and thin film forms in order to elucidate their structural stability and dynamic behavior. An efficient workflow, consisting of quantum mechanical simulations (based on density functional theory) and machine learning interatomic potentials (MTPs) coupled to temperature-dependent effective potentials (TDEP) and classical molecular dynamics, was verified necessary to capture the temperature-dependent stabilization of the phonons in bulk ice Ih and XI. Anharmonicity and nuclear quantum effects, incorporated in an efficient way through a quantum thermal bath technique, were found crucial to dynamically stabilize low-frequency lattice modes and high-frequency vibrational stretching involving hydrogen. We have identified three novel thin film structures that retain their stability up to at least 250 K and have shed light on their phonon characteristics. In addition, our examination of the Raman spectrum of ice underscores the shortcomings of predicting vibrational properties when relying entirely on the harmonic approximation or purely anharmonic effects. The corrected redistribution of vibrational intensities is found to be achieved only upon inclusion of quantum nuclear vibrations. This was found to be even more crucial for low-dimensional thin film (2D) structures. Overall, our findings demonstrate the significance of joining advanced computational methodologies in unraveling the intricate vibrational dynamics of crystalline ice materials, offering valuable insights into their thermodynamic and structural properties. Furthermore, we suggest a procedure based on MTPs coupled to a quantum thermal bath for the computationally efficient probing of nuclear effects in ice structures, although equally applicable to any other system.

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or complexes thereof, are strongly influenced by protonation changes of their typically many titratable groups, which explains their sensitivity to pH changes. Conversely, conformational and environmental changes of the biomolecule affect the protonation state of these groups. With few exceptions, conventional force field-based molecular dynamics (MD) simulations neither account for these effects nor do they allow for coupling to a pH buffer. Here, we present design decisions and applications of a rigorous Hamiltonian interpolation λ-dynamics constant pH method in GROMACS, which rests on GPU-accelerated Fast Multipole Method (FMM) electrostatics. Our implementation supports both CHARMM36m and Amber99sb*-ILDN force fields and is largely automated to enable seamless switching from regular MD to constant pH MD, involving minimal changes to the input files. Here, the first of two companion papers describes the underlying constant pH protocol and sample applications to several prototypical benchmark systems such as cardiotoxin V, lysozyme, and staphylococcal nuclease. Enhanced convergence is achieved through a new dynamic barrier height optimization method, and high pK_a accuracy is demonstrated. We use Functional Mode Analysis (FMA) and Mutual Information (MI) to explore the complex intra- and intermolecular couplings between the protonation states of titratable groups as well as those between protonation states and conformational dynamics. We identify striking conformation-dependent pK_a variations and unexpected inter-residue couplings. Conformation-protonation coupling is identified as a primary cause of the slow protonation convergence notorious to constant pH simulations involving multiple titratable groups, suggesting enhanced sampling methods to accelerate convergence.

This study focuses on numerical methods to compute charge carrier mobility in disordered materials, such as organic light-emitting diodes (OLEDs), based solely on molecular structures. The approach involves developing an ab initio method for calculating charge carrier mobility in organic materials using kinetic Monte Carlo (KMC) simulations. These simulations utilize Marcus rates derived from precise calculations of transfer integrals and site energies specific to the material's morphology. Going beyond the current approach to tackle a multicharge model system presents computational challenges, particularly in calculating site energies. To address this issue, a novel lattice mapping method was developed to efficiently determine transfer integrals and site energies from realistic morphologies while keeping computational costs manageable. Validation of the method was conducted by comparing mobility values computed using the KMC method with experimental data, showing a good agreement. Further insights into charge transport dynamics were gained through the analysis of charge carrier behavior using residence time calculations. Additionally, the model's applicability to multicharge systems was demonstrated by simulating exciton formation. In conclusion, the model has the potential to effectively and accurately simulate charge carrier trajectories in multicharge, multilayer models with minimal loss of information from realistic morphology, making it a valuable tool for designing and optimizing organic electronic devices.

We developed a k-point sampling real-time TDDFT (RT-TDDFT) program within the Gaussian and plane waves (GPW) framework of the CP2K software suite. In addition to standard real-time propagation of time-dependent Kohn-Sham orbitals, we make use of symmetry-based k-point reduction and k-point parallelization schemes so that our RT-TDDFT program in the GPW framework is feasible for practical large-scale calculations. We also implemented DFT + U as a relevant extension for real-time simulations of systems with strong electron correlations. In particular, we extended the "tensorial" subspace representation approach for DFT + U, following the formulation in [Chai, Z., et al. J. Chem. Theory Comput., 2024, 20, 8984], to k-point sampling RT-TDDFT. Our extension, which is, to our knowledge, the first application of the "tensorial" subspace representation approach to k-point sampling RT-TDDFT, is found to be robust and efficient with small additional costs owing to the locality of Gaussian basis functions, indicating that it is a promising approach to RT-TDDFT + U for solids. We show details of our implementation in CP2K and the results of our benchmark calculations.

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or their complexes, are strongly influenced by protonation changes of their typically many titratable groups, which explains their pH sensitivity. Conversely, conformational and environmental changes in the biomolecule affect the protonation state of these groups. With a few exceptions, conventional force field-based molecular dynamics (MD) simulations do not account for these effects, nor do they allow for coupling to a pH buffer. The λ-dynamics method implements this coupling and thus allows for MD simulations at constant pH. It uses separate Hamiltonians for the protonated and deprotonated states of each titratable group, with a dynamic λ variable that continuously interpolates between them. However, rigorous implementations of Hamiltonian Interpolation (HI) λ-dynamics are prohibitively slow for typical numbers of sites when used with particle mesh Ewald (PME). To circumvent this problem, it has recently been proposed to interpolate the charges (QI) instead of the Hamiltonians. Here, in the second of two companion papers, we propose a rigorous yet efficient Multipole-Accelerated Hamiltonian Interpolation (MAHI) method to perform λ-dynamics in GROMACS. Starting from a charge-scaled Hamiltonian, precomputed with the Fast Multipole Method (FMM), the correct HI forces are calculated with negligible computational overhead. However, other electrostatic solvers, such as PME, can also be used for the precomputation. We compare Hamiltonian interpolation with charge interpolation and show that HI leads to more frequent transitions between protonation states, resulting in better sampling and accuracy. Our accuracy and performance benchmarks show that introducing, e.g., 512 titratable sites to a one million atom MD system increases runtime by less than 20% compared to a regular FMM-based simulation. We have integrated the scheme into our GPU-accelerated FMM code for the simulation software GROMACS, allowing easy and effortless transitions from standard force field simulations to constant pH simulations.

This study extends the Grassmann extrapolation (G-Ext) method, which was introduced for Born-Oppenheimer molecular dynamics, to the context of geometry optimization. Using density matrices from previous optimization steps, the G-Ext approach applies a nonlinear, structure-preserving mapping onto the Grassmann manifold to provide an initial guess which accelerates the convergence of the self-consistent field (SCF) procedure. Using the optimal parameters identified by employing various descriptors and computational strategies across a diverse set of molecules, G-Ext shows excellent performance improvements, particularly with large molecular systems.

The dynamics of metal centers are challenging to describe due to the vast variety of ligands, metals, and coordination spheres, hampering the existence of general databases of transferable force field parameters for classical molecular dynamics simulations. Here, we present easyPARM, a Python-based tool that can calculate force field parameters for a wide range of metal complexes from routine frequency calculations with electronic structure methods. The approach is based on a unique labeling strategy, in which each ligand atom that coordinates the metal receives a unique atom type. This design prevents parameter shortage, labeling duplication, and the necessity to post-process output files, even for very complicated coordination spheres, whose parametrization process remain automatic. The program requires the Cartesian Hessian matrix, the geometry xyz file, and the atomic charges to provide reliable force-field parameters extensively benchmarked against density functional theory dynamics in both the gas and condensed phases. The procedure allows the classical description of metal complexes at a low computational cost with an accuracy as good as the quality of the Hessian matrix obtained by quantum chemistry methods. easyPARM v2.00 reads vibrational frequencies and charges in Gaussian (version 09 or 16) or ORCA (version 5 or 6) format and provides refined force-field parameters in Amber format. These can be directly used in Amber and NAMD molecular dynamics engines or converted to other formats. The tool is available free of charge in the GitHub platform (https://github.com/Abdelazim-Abdelgawwad/easyPARM.git).

This study explores the oxygen-binding mechanism and the potential peroxo-to-bis-μ-oxo isomerization in hemocyanin (Hc) using a quantum mechanics/molecular mechanics (QM/MM) approach at the multireference NEVPT2 level of theory (QM[NEVPT2]/MM). Our results support the previously proposed mechanism for Hc oxygen binding, involving two nearly simultaneous electron-transfer (ET) steps and a triplet-singlet intersystem crossing (ISC). However, we find that the first ET step occurs prior to ISC, resulting in the formation of a stable singlet superoxide intermediate through a low-energy barrier. The second ET leads to the formation of a singlet oxy-hemocyanin species featuring the characteristic peroxo-Cu₂O₂ "butterfly" core. Moreover, QM[NEVPT2]/MM simulations reveal a lower-energy barrier for the peroxo-to-bis-μ-oxo isomerization compared with density functional theory (DFT), although the peroxo form remains energetically favored within the protein environment. These findings offer new insights into the behavior of the hemocyanin active site, highlighting the importance of considering both the electronic correlation and the protein environment in accurately modeling copper-oxygen interactions in biological systems.

Graph neural network (GNN) architectures have emerged as promising force field models, exhibiting high accuracy in predicting complex energies and forces based on atomic identities and Cartesian coordinates. To expand the applicability of GNNs, and machine learning force fields more broadly, optimizing their computational efficiency is critical, especially for large biomolecular systems in classical molecular dynamics simulations. In this study, we address key challenges in existing GNN benchmarks by introducing a dataset, DISPEF, which comprises large, biologically relevant proteins. DISPEF includes 207,454 proteins with sizes up to 12,499 atoms and features diverse chemical environments, spanning folded and disordered regions. The implicit solvation free energies, used as training targets, represent a particularly challenging case due to their many-body nature, providing a stringent test for evaluating the expressiveness of machine learning models. We benchmark the performance of seven GNNs on DISPEF, emphasizing the importance of directly accounting for long-range interactions to enhance model transferability. Additionally, we present a novel multiscale architecture, termed Schake, which delivers transferable and computationally efficient energy and force predictions for large proteins. Our findings offer valuable insights and tools for advancing GNNs in protein modeling applications.