Journal of Chemical Theory and Computation最新文献

The polarizable continuum model (PCM) is a computationally efficient way to incorporate dielectric boundary conditions into electronic structure calculations, via a boundary-element reformulation of Poisson's equation. This transformation is only rigorously valid for an isotropic dielectric medium. To simulate anisotropic solvation, as encountered at an interface or when parts of a system are solvent-exposed while other parts are in a nonpolar environment, ad hoc modifications to the PCM formalism have been suggested, in which a dielectric constant is assigned separately to each atomic sphere that contributes to the solute cavity. The accuracy of this "heterogeneous" PCM (HetPCM) method is tested here for the first time, by comparison to results from a generalized Poisson equation solver. The latter is a more expensive and cumbersome approach to incorporate arbitrary dielectric boundary conditions, but one that corresponds to a well-defined scalar permittivity function, ε(r). We examine simple model systems for which a function ε(r) can be constructed in a manner that maps reasonably well onto a dielectric constant for each atomic sphere, using a solvent-exposed dielectric constant ε_solv = 78 and a range of smaller values to represent hydrophobic environments. For nonpolar dielectric constants ε_nonp ≤ 2, differences between the HetPCM and Poisson solvation energies are large compared to the effect of anisotropy on the solvation energy. For ε_nonp = 4 and ε_nonp = 10, however, HetPCM and anisotropic Poisson solvation energies agree to within 2 kcal/mol in most cases. As a realistic use case, we apply the HetPCM method to predict solvation energies and pK_a values for blue copper proteins. The HetPCM method affords pK_a values that are more in line with experimental results as compared to either gas-phase calculations or homogeneous (isotropic) PCM results.

The weighted ensemble (WE) algorithm is gaining popularity as a rare event method for studying long timescale processes with molecular dynamics. WE is particularly useful for determining kinetic properties, such as rates of protein (un)folding and ligand (un)binding, where transition rates can be calculated from the flux of trajectories into a target basin of interest. However, this flux depends exponentially on the number of splitting events that a given trajectory experiences before reaching the target state and can vary by orders of magnitude between WE replicates. Markov state models (MSMs) are helpful tools to aggregate information across multiple WE simulations and have previously been shown to provide more accurate transition rates than WE alone. Discrete-time MSMs are models that coarsely describe the evolution of the system from one discrete state to the next using a discrete lag time, τ. When an MSM is built using conventional MD data, longer values of τ typically provide more accurate results. Combining WE simulations with Markov state modeling presents some additional challenges, especially when using a value of τ that exceeds the lag time between resampling steps in the WE algorithm, τ_WE. Here, we identify a source of bias that occurs when τ > τ_WE, which we refer to as "merging bias". We also propose an algorithm to eliminate the merging bias, which results in merging bias-corrected MSMs, or "MBC-MSMs". Using a simple model system, as well as a complex biomolecular example, we show that MBC-MSMs significantly outperform both τ = τ_WE MSMs and uncorrected MSMs at longer lag times.

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.

The MBX software provides an advanced platform for molecular dynamics simulations, leveraging state-of-the-art MB-pol and MB-nrg data-driven many-body potential energy functions. Developed over the past decade, these potential energy functions integrate physics-based and machine-learned many-body terms trained on electronic structure data calculated at the "gold standard" coupled-cluster level of theory. Recent advancements in MBX have focused on optimizing its performance, resulting in the release of MBX v1.2. While the inherently many-body nature of MB-pol and MB-nrg ensures high accuracy, it poses computational challenges. MBX v1.2 addresses these challenges with significant performance improvements, including enhanced parallelism that fully harnesses the power of modern multicore CPUs. These advancements enable simulations on nanosecond time scales for condensed-phase systems, significantly expanding the scope of high-accuracy, predictive simulations of complex molecular systems powered by data-driven many-body potential energy functions.

This work presents an efficient approach to optimizing force field parameters for sulfone molecules using a combination of genetic algorithms (GA) and Gaussian process regression (GPR). Sulfone-based electrolytes are of significant interest in energy storage applications, where accurate modeling of their structural and transport properties is essential. Traditional force field parametrization methods are often computationally expensive and require extensive manual intervention. By integrating GA and GPR, our active learning framework addresses these challenges by achieving optimized parameters in 12 iterations using only 300 data points, significantly outperforming previous attempts requiring thousands of iterations and parameters. We demonstrate the efficiency of our method through a comparison with state-of-the-art techniques, including Bayesian Optimization. The optimized GA-GPR force field was validated against experimental and reference data, including density, viscosity, diffusion coefficients, and surface tension. The results demonstrated excellent agreement between GA-GPR predictions and experimental values, outperforming the widely used OPLS force field. The GA-GPR model accurately captured both bulk and interfacial properties, effectively describing molecular mobility, caging effects, and interfacial arrangements. Furthermore, the transferability of the GA-GPR force field across different temperatures and sulfone structures underscores its robustness and versatility. Our study provides a reliable and transferable force field for sulfone molecules, significantly enhancing the accuracy and efficiency of molecular simulations. This work establishes a strong foundation for future machine learning-driven force field development, applicable to complex molecular systems.

Peptide-based coacervates are crucial for drug delivery due to their biocompatibility, versatility, high drug loading capacity, and cell penetration rates; however, their stability mechanism and phase behavior are not fully understood. Additionally, although Martini is one of the most famous force fields capable of describing coacervate formation with molecular details, a comprehensive benchmark of its accuracy has not been conducted. This research utilized the Martini 3.0 force field and machine learning algorithms to explore representative peptide-based coacervates, including those composed of polyaspartate (PAsp)/polyarginine (PArg), rmfp-1, sticker-and-spacer small molecules, and HBpep molecules. We identified key coacervate formation driving forces such as Coulomb, cation-π, and π-π interactions and established three criteria for determining coacervate formation in simulations. The results also indicate that while Martini 3.0 accurately captures coacervate formation trends, it tends to underestimate Coulomb interactions and overestimate π-π interactions. What is more, our study on drug encapsulation of HBpep and its derivative coacervates suggested that most loaded drugs were distributed on surfaces of HBpep clusters, awaiting experimental validation. This study employs simulation to enhance understanding of peptide-based coacervate phase behavior and stability mechanisms while also benchmarking Martini 3.0, thereby providing fundamental insights for future experimental and simulation investigations.

We investigate the convergence of quasiparticle energies for periodic systems to the thermodynamic limit using increasingly large simulation cells corresponding to increasingly dense integration meshes in reciprocal space. The quasiparticle energies are computed at the level of equation-of-motion coupled-cluster theory for ionization (IP-EOM-CC) and electron attachment processes (EA-EOM-CC). By introducing an electronic correlation structure factor, the expected asymptotic convergence rates for systems with different dimensionality are formally derived. We rigorously test these derivations through numerical simulations for trans-polyacetylene using IP/EA-EOM-CCSD and the G₀W₀@HF approximation, which confirm the predicted convergence behavior. Our findings provide a solid foundation for efficient schemes to correct finite-size errors in IP/EA-EOM-CCSD calculations.

We explore the melting of a lattice DNA in the presence of atmospheric disorder, which mimics the crowded environment inside the cell nucleus, using Monte Carlo simulations. The disorder is modeled by randomly retaining lattice sites with probability p while diluting the rest, rendering them unavailable to the DNA. By varying the disorder over a wide range from p = 1 (zero disorder) up to the percolation critical point p_c = 0.3116, we show the melting temperature (T_m) to increase nearly linearly with disorder up to p ≈ 0.6, while strong nonlinearity enters for p ≲ 0.6. Associated changes in the bubble statistics have been investigated, showing a substantial change in the bubble size exponent at corresponding melting points for p ≤ 0.5. Based on these findings, two distinct disorder regimes showing weak and strong effects on melting have been identified. For simulations, we use the pruned and enriched Rosenbluth method in conjunction with a depth-first implementation of the Leath algorithm to generate the underlying disorder.

The description of electron-phonon coupling in materials is complex, with varying definitions of coupling constants in the literature and different theoretical approaches available. This article analyzes different levels of theory to introduce and compute these coupling constants. Within the quasi-particle picture, we derive an effective linear-coupling Hamiltonian, describing the interaction of electronic quasi-particles with vibrations. This description allows a comparison between coupling constants computed using density functional theory and higher-level quasi-particle approaches by identifying the Kohn-Sham potential as an approximation to the frequency-independent part of the self-energy. We also investigate their dependence on the exchange-correlation (XC) functional. Despite significant deviations of the Kohn-Sham eigenvalues, which arise from different XC functionals, the resulting coupling constants are remarkably similar. A comparison to quasi-particle methods, such as the well-established G0W0 approach, reveals significant quasi-particle weight renormalization. Surprisingly, however, in nearly all the considered cases, the coupling constants computed in the DFT framework are excellent approximates of the ones in the quasi-particle framework, which is traced back to a significant cancellation of competing terms. Other quasi-particle methods, such as the Outer Valence Green's Function approach and the ΔSCF method, are also included in the comparison. Moreover, we investigate the coupling of vibrations to excitonic excitations and find, by comparison to time-dependent density functional theory and extended multiconfiguration quasi-degenerate second-order perturbation theory, that knowing the underlying electron- and hole-vibration couplings is sufficient to accurately determine the exciton-vibration coupling constants in the studied cases.