Journal of Chemical Theory and Computation最新文献

Free energy calculations underpin a wide spectrum of tasks in computational chemistry, from scanning free energy surfaces along collective variables to ranking affinities of ligands in computer-aided drug discovery. Their routine use, however, is constrained by the high cost of sampling an extensive set of thermodynamic states with replicas. We introduce here Integrated Boltzmann Sampling (IBS), a few-state framework that integrates multistate thermodynamic sampling into a small set of artificial ensembles. Trajectories generated from these ensembles are reweighted to recover the full thermodynamic information for dozens of alchemical and tempered states, reducing the formal sampling cost from K · S to (1 – ϵ + Kϵ) · S, ϵ ≪ 1, where K counts the number of thermodynamic states involved, and S represents the computational effort required to sample a single state. On the SAMPL6 host–guest benchmark and a 13-ligand Farnesoid X receptor panel, IBS achieves accuracy comparable to replica-based free energy methods while lowering wall time by approximately 50–60%. These results demonstrate that achieving chemically accurate free energy predictions does not require exhaustive replica sampling and that IBS offers an efficient, drop-in alternative for computational applications that rely on accurate free energy differences.

G Protein-Coupled Receptors (GPCRs) mediate signal transduction across cellular membranes and are major drug targets. Activation of these receptors upon binding of an extracellular ligand involves propagation of structural change across the transmembrane domain to the cytoplasmic G protein partner, a process generally thought to involve dynamic hydrogen(H)-bond networks. Here we present DNET, a graph-based tool and workflow that enables efficient computations of dynamic protein–water H-bond networks. DNET is a fully portable Python tool that reads simulation trajectories, computes graphs of the dynamic protein–water H-bond networks, and generates, for each H-bonding residue, a residue summary that includes water interactions, H-bond time series, histograms, potential of mean force estimates, and the number of conformations of the H-bond. To facilitate estimates of pK_a fluctuations within the H-bond network, DNET calls PROPKA and computes, for each titratable residue that is part of the H-bond network, time series and analyses of the pK_a estimates. To illustrate the usefulness of DNET we apply it to study the wild-type and two mutations of jumping spider rhodopsin 1, JSR-1, a visual rhodopsin GPCR activated by the photoisomerization of the covalently bound retinal chromophore. The UV–vis data we present here demonstrate that the mutated JSR-1 proteins express, but both have an altered electrostatic environment of the retinal Schiff base. The DNET analyses indicate a highly complex dynamics of the retinal H-bond network, with some H-bonds that have only one conformational mode, and other H-bonds with multiple conformational modes separated by small energy barriers, and pK_a fluctuations that associate with the H-bond dynamics. The mutations associate with an altered H-bond network of the retinal Schiff base.

Channelrhodopsin-2 (ChR2) is a light-gated ion channel widely used in optogenetics, a technique that enables precise control of neuronal activity by genetically engineering light-sensitive proteins into cell membranes. This protein exists in dimeric form, with each monomer containing a retinal Schiff base (RSB) moiety covalently bonded that undergoes trans–cis isomerization upon light absorption. However, the limited penetration depth of visible light in biological tissues motivates the use of multiphoton-absorption techniques, which enhance tissue penetration, improve focality, and reduce phototoxicity, thereby offering a promising alternative for optogenetic applications. In this paper, we present a fully atomistic multiscale methodology for computing the one-, two-, and three-photon absorption spectra of ChR2, where the protein, lipid bilayer, and solvent are explicitly considered throughout the workflow. This methodology integrates classical molecular mechanics (MM) molecular dynamics (MD), quantum mechanics/molecular mechanics (QM/MM)-MD, and fragment-based polarizable embedding (PE) to derive environment-specific PE potentials from the explicit protein–lipid-solvent environment. The final step in the methodology is to use these potentials to compute accurate spectra via PE-time-dependent density functional theory (PE-TD-DFT). Validation against experimental one-photon absorption spectra demonstrates excellent agreement. For the first time, we report the theoretical two- and three-photon absorption in ChR2, albeit without direct experimental comparison. We compare the multiphoton absorption (MPA) spectra where the two RSB moieties are sampled using classical MD and QM/MM-MD, respectively. The resulting spectral differences are attributed to variations in key structural parameters that we analyze and document.

We present a spin-permutation diabatization strategy that transforms ab initio spin-pure eigenstates into spin-localized diabatic states, enabling direct mapping to spin-effective Hamiltonians without projection or orbital localization. The method provides both a real-space decomposition of electronic states in terms of localized spins and a straightforward evaluation of exchange couplings. Applications to several representative systems, including ethylene torsion, prototypical diradicals (benzynes, xylylenes, methylene), trimethylenebenzene triradical, singlet–triplet excited states of organic chromophores, and triplet-pair states in a tetracene dimer, demonstrate that the approach provides magnetic couplings and affords a clear physical interpretation of interacting spins. This general and conceptually transparent framework bridges ab initio electronic structure theory and spin models, and is expected to be especially valuable for systems with nontrivial distributions of unpaired electrons, such as delocalized or strongly correlated molecular magnets and spin-active chromophores.

RNAs are key to understanding cellular mechanisms and a prime target for novel therapeutic interventions. However, RNA structural ensembles are notoriously difficult to study due to structural polymorphism and their highly dynamic nature, i.e., the fact that RNAs can adopt multiple structures and the transitions between them are fast. As a result, computational and experimental methods to study the RNA structure are limited in their usefulness, and often, a combination of multiple methods is required. Furthermore, RNA force fields are not yet developed to the same standard as those for proteins. Here, we demonstrate that energy landscape explorations via discrete path sampling enable a full mapping of structural ensembles and can capture mutational changes well. In this contribution, we show that the ensembles derived for the TAR stemloop and a derived stemloop (ES2) are sufficient to reproduce experimental observations without the need to introduce experimental data into the modeling beyond the force field parameters. Our modeling reveals significant complexity in both the structural ensemble and the activation pathway. Furthermore, we identify a transient binding pocket that emerges on the activation pathway.

We introduce a new interlayer potential (ILP) for simulating the adsorption and diffusion of planar organic molecules with partially charged heteroatoms on hexagonal boron nitride (hBN). Unlike previous models, this ILP incorporates all-atom electrostatic interactions alongside short-range repulsion and long-range attraction, enabling the accurate treatment of polar molecules. Parametrized against density functional theory data for pentacene and PTCDI, the ILP demonstrates transferability to related systems such as PTCDA. Comparative studies of nonpolar pentacene and polar PTCDA reveal distinct behaviors in single-molecule diffusion, cluster formation, and monolayer growth. PTCDA exhibits stronger binding due to electrostatic contributions, limiting diffusion to short-ranged hops, while pentacene undergoes long-range translocations facilitated by out-of-plane motions. At low coverage, PTCDA molecules lock into place via carbonyl-mediated hydrogen bonds, enabling only collective motion, whereas pentacene remains mobile. Monolayer simulations reproduce experimentally observed epitaxial morphologies: PTCDA forms a dense square lattice, while pentacene aligns parallel along its long axis. This ILP offers a computationally efficient and accurate alternative to ab initio and machine-learning methods, opening avenues for modeling polar organic molecules on hBN. Its utility extends to understanding layer formation and structural properties in hBN-encapsulated or -supported organic systems.

Retrosynthesis, the process of deconstructing a target molecule into simpler precursors, is crucial for computer-aided synthesis planning (CASP). Widely adopted tree-search methods often suffer from exponential computational complexity. In this work, we introduce FragmentRetro, a novel retrosynthetic method that leverages fragmentation algorithms, specifically BRICS and r-BRICS, combined with stock-aware exploration and pattern fingerprint screening to achieve quadratic complexity. FragmentRetro recursively combines molecular fragments and verifies their presence in a building block set, providing sets of fragment combinations as retrosynthetic solutions. We present the first formal computational analysis of retrosynthetic methods, showing that tree search exhibits exponential complexity $O\left(b^h\right)$, DirectMultiStep scales as $O\left(h^6\right)$, and FragmentRetro achieves $O\left(h^2\right)$, where h represents the number of heavy atoms in the target molecule and b is the branching factor for tree search. Evaluations on PaRoutes, USPTO-190, and natural products demonstrate that FragmentRetro achieves high solved rates with competitive runtime, including cases where tree search fails. The method benefits from fingerprint screening, which significantly reduces substructure matching complexity. While FragmentRetro focuses on efficiently identifying fragment-based solutions rather than full reaction pathways, its computational advantages and ability to generate strategic starting candidates establish it as a powerful foundational component for scalable and automated synthesis planning.

Atomistic models of water with increased dispersion interactions (e.g., OPC, TIP4P-D) have been proposed to produce extended conformational ensembles of intrinsically disordered proteins (IDPs). However, the role of the protein force field in obtaining accurate macromolecular ensembles of IDPs remains unclear. Isolating the influence of the protein and water models by comparison to an experimental measure (such as X-ray scattering) requires an atomic consideration of the hydration layer waters around a thermally fluctuating protein. To enable an atomically detailed scattering calculation around a thermally fluctuating solute, we have developed a new scattering model termed small and wide angle X-ray scattering for all molecular dynamics engines (SWAXS-AMDE). SWAXS-AMDE can handle the trajectory files from all of the popular molecular dynamics (MD) simulation softwares, thus facilitating a straightforward validation of force field improvements. SWAXS-AMDE computed scattering profiles for polyampholyte peptides show the AMBER ff19SB protein force field to play a crucial role in obtaining accurate macromolecular ensembles of both folded proteins and IDPs.

We have developed Si parameters for the orthogonalization and dispersion corrected semiempirical method, ODM3. The new parameter set is called ODM3.25. As an illustration of the new capability allowed by these new parameters, we simulated the nonadiabatic dynamics of the cubic hydrogen silsesquioxane (Si₈O₁₂H₈, HSQ) electron beam photoresist by using curvature-driven trajectory surface hopping with energy-based decoherence (κTSH-EDC) interfaced with ODM3.25. Our simulations involve four coupled singlet states and four coupled potential energy surfaces. The isomerization of the HSQ molecule involves an interconversion of kinetic energy and potential energy during the first 100 fs, and this is followed by Si–H bond breaking and O–H bond formation. We find that about a quarter of the product remains in an electronically excited state. We anticipate that the combination of the cost-effective ODM3.25 with curvature-driven algorithms for electronically nonadiabatic dynamics will allow simulations of nonadiabatic dynamics in a broad range of problems involving Si-containing nanoparticles.

Interpreting time-resolved magnetic resonance experiments, sensitive to slow motions in molecules, requires access to at least the microsecond time scale. Today, all-atom classical molecular dynamics simulations allow exploration of such a long time scale; however, this comes at the price of a considerable computational effort. Stochastic models, based on a hierarchical distinction of the coordinates into relevant (treated explicitly) and irrelevant (treated as generators of fluctuation and dissipation), offer a relatively low-cost solution to this problem. In the past, ad hoc but essentially phenomenological approaches based on Langevin or Fokker–Planck equations have been employed, which are good in catching relevant differences among (even complex) molecular systems, but lack of predictive power since a map between such parameters and atomistic details is not always clear or defined. Recently, a rigorous derivation of a stochastic description of the dynamics of a macromolecule from the complete equations of motion has been provided. In this paper, a computational strategy based on the solution of the Brownian dynamics equations associated with the original model is discussed for the calculation and interpretation of nuclear magnetic resonance relaxation data. The approach merges the ability of stochastic approaches to perform a targeted complexity reduction of the system with the flexibility of molecular dynamics simulations in describing at the atomistic level the time evolution of the system. By expressing the stochastic dynamics in the relevant natural internal coordinates and exploiting the acceleration power of GPU-based hardware, the proposed approach lays the foundations for an effective interpretation of long-time dynamics of generic semiflexible complex molecules.