We present an open Molecular Crystal (MC) database of Machine-Learned Interatomic Potentials (MLIP) called MolCryst-MLIPs. The first release comprises fine-tuned MACE models for nine molecular crystal systems – Benzamide, Benzoic acid, Coumarin, Durene, Isonicotinamide, Niacinamide, Nicotinamide, Pyrazinamide, and Resorcinol – developed using the Automated Machine Learning Pipeline (AMLP), which streamlines the entire MLIP development workflow, from reference data generation to model training and validation, into a reproducible and user-friendly pipeline. Models are fine-tuned from the MACE-MH-1 foundation model (omol head), yielding a mean energy MAE of 0.141 kJ/mol/atom and a mean force MAE of 0.648 kJ/mol/Angstrom across all systems. Dynamical stability and structural integrity, as assessed through energy conservation, P2 orientational order parameters, and radial distribution functions, are evaluated using molecular dynamics simulations. The released models and datasets constitute a growing open database of validated MLIPs, ready for production MD simulations of molecular crystal polymorphism under different thermodynamic conditions.

Accurate solvation free energies from molecular dynamics simulations require efficient sampling of coupled slow variables, including solvent coordinates, solute conformational modes, and the alchemical coordinate λ. Here, we develop a λ-dynamics framework that combines mass scaling, on-the-fly probability enhanced sampling (OPES), and driven adiabatic free energy dynamics (d-AFED) to address these sampling challenges within a unified protocol. For rigid organic solutes, Hamiltonian replica exchange with mass scaling is first used to quantify the effect of octanol solvent relaxation. Reducing all octanol atomic masses by a factor of ten accelerates convergence by more than fivefold while preserving equilibrium solvation free energies. These calculations then provide reference benchmarks for λ-OPES, a dual-bias λ-dynamics strategy that combines the "standard" and "explore" variants of OPES to promote transitions along the alchemical coordinate. This approach reaches convergence on timescales comparable to replica exchange, but without predefined λ windows or multiple parallel simulations. For flexible N-acetyl amino-acid amide solutes, λ-OPES is coupled with d-AFED on selected backbone and side-chain dihedrals to enable simultaneous alchemical and conformational enhanced sampling. This combined strategy improves agreement with experimental octanol-water partition coefficients and reduces the mean absolute error from 0.75 log units with λ-OPES alone to 0.30 log units with λ-OPES-d-AFED. Overall, this work establishes an integrated enhanced sampling protocol for solvation free energy calculations across rigid organic solutes and flexible peptide-like solutes, and provides a foundation for the application of alchemical free energy methods to larger and more conformationally complex systems.

We present MXtalTools, a flexible Python package for the data-driven modeling of molecular crystals, facilitating machine learning studies of the molecular solid state. MXtalTools comprises several classes of utilities: (1) synthesis, collation, and curation of molecule and crystal data sets, (2) integrated workflows for model training and inference, (3) crystal parametrization and representation, (4) crystal structure sampling and optimization, (5) end-to-end differentiable crystal sampling, construction, and analysis. Our modular functions can be integrated into existing workflows or combined and used to build novel modeling pipelines. MXtalTools leverages CUDA acceleration to enable high-throughput crystal modeling. The Python code is available open-source on our GitHub page, with detailed documentation on ReadTheDocs.

Molecular crystals are a highly polymorphic class of materials, with a single molecule commonly crystallizing via multiple packing patterns, making structure and property prediction very challenging. Crystal structure prediction typically comprises the production of sets of promising candidate structures, each considered in isolation rather than as samples in a thermodynamic distribution. Likewise, modern generative approaches to this problem, despite naturally sampling distributions of crystals, lack a concrete formulation of the distributions being sampled. Two components are required to impart meaning to the distributions of crystals generated under such models: a canonical parameterization, and a loss function which equilibrates the generated samples to some target distribution. We develop such a parameterization, and train energy-based generative flow networks (GFlowNets) to approximate the Boltzmann distribution over crystal structures for target molecules and space groups. Combined, these components comprise our MXtalGFlow framework for molecular crystal modeling. Going beyond sampling disconnected sets of low-energy structures, MXtalGFlow yields a thermodynamic distribution over crystal structures. We sample and analyze distributions of crystals for two molecules, each under two energy functions, a Lennard-Jones potential and the Universal Model for Atoms. We characterize the local structural basins about the known polymorphs, and identify additional as-yet un-reported packing modes with competitive probabilities to the known experimental structures. With MXtalGFlow, we illustrate how to define and train a model to sample a thermodynamically meaningful distribution of molecular crystals, and analyze such a distribution to glean useful information.

We introduce a computationally simple algorithm for sampling open-chain distributions within the framework of imaginary-time Feynman path integration. The present method is based on the staging algorithm introduced by Pollock and Ceperley [Phys. Rev. B 30, 2555 (1984)] originally developed for computing position-dependent observables. Here, we sample off-diagonal elements of the density matrix, formulated as a distribution describing a linear polymer-like chain of beads, each connected via nearest-neighbor springs to calculate momentum-dependent quantities. This is achieved using a Monte Carlo scheme that ensures efficient and unbiased sampling of all beads along the chain from the free-particle distribution via a staging transformation; we refer to this approach as staging open path integral Monte Carlo (OPIMC). The proposed algorithm is straightforward to implement, as it only involves sampling Gaussian distributions through a transformation defined by a set of recursion relations, followed by a standard Metropolis acceptance/rejection step. The staging OPIMC method accurately reproduces end-to-end and momentum distributions for quantum systems ranging from coupled harmonic oscillators to liquid water.

Machine-learning interatomic potentials (MLIPs) have enabled molecular dynamics at near ab initio accuracy, yet remain limited to energies and forces by construction, leaving electronic observables such as dipole moments and polarizabilities inaccessible. We introduce DenSNet, a density-first approach to machine-learned electronic structure that learns the Hohenberg–Kohn map from nuclear configurations to the ground-state electron density. Our approach employs an SE(3)-equivariant neural network to predict density coefficients of a flexible atom-centered Gaussian basis, combined with a Δ-learning strategy that uses superposed atomic densities as a prior to accelerate training. A second equivariant network then maps the predicted density to the total energy, providing a unified framework for molecular dynamics and electronic structure. We validate DenSNet on ethanol, ethanethiol, and resorcinol, where infrared spectra from machine-learned trajectories show excellent agreement with experimental gas-phase measurements. To test scalability, we train on polythiophene oligomers with 1–6 monomers and extrapolate to chains of up to 12 monomers, generating stable long-time trajectories whose infrared spectra agree with reference density functional theory calculations. Here, we show that reinstating the electron density as the central learned quantity opens a practical route to transferable prediction of spectroscopic and electronic observables in large-scale molecular simulations.

Machine learning interatomic potentials (MLIPs) have become powerful tools to extend molecular simulations beyond the limits of quantum methods, offering near-quantum accuracy at much lower computational cost. Yet, developing reliable MLIPs remains difficult because it requires generating high-quality data sets, preprocessing atomic structures, and carefully training and validating models. In this work, we introduce an Automated Machine Learning Pipeline (AMLP) that unifies the entire workflow from data set creation to model validation. AMLP employs large-language-model agents to assist with electronic-structure code selection, input preparation, and output conversion, while its analysis suite (AMLP-Analysis) based on ASE supports a range of molecular simulations. The pipeline is built on the MACE architecture and validated on acridine polymorphs, where with a straightforward fine-tuning of a foundation model mean absolute errors of 1.7 meV/atom in energies and 7.0 meV/Å in forces are achieved. The fitted MLIP reproduces DFT geometries with sub-Å accuracy and demonstrates stability during molecular dynamics simulations in the microcanonical and canonical ensemble.

Structural, thermal, and dynamic properties of four deep eutectic solvents comprising choline chloride paired with ortho-phenolic derivative hydrogen-bond donors were probed using experiments and molecular simulations. The hydrogen-bond donors include phenol, catechol, o-chlorophenol, and o-cresol, in a 3:1 mixture with the hydrogen-bond acceptor choline chloride. Density, viscosity, and pulsed-field gradient NMR diffusivity measurements were conducted over a range of temperatures. Classical and ab initio molecular dynamics simulation results match experimental data reasonably well. The simulation results were then used to perform a more detailed analysis of the local structure and dynamics of these systems.

Candidate systems for next-generation battery electrolyte materials, such as deep eutectic solvents and ionic liquids, often suffer from the limitation of an inverse relation that exists between viscosity and conductivity, known as Walden’s rule, which can suppress rates of charge transport and limit their performance characteristics. A strategy for circumventing this problem draws its inspiration from the world of fuel-cell based ion exchange membranes and the types of charge transport processes operative in these systems. In this talk, I will describe a project aimed at leveraging machine learning and Feynman path-integral based quantum simulation strategies, in combination with experimental synthesis and characterization, to develop a novel class of battery electrolytes that demonstrates an ability to escape the limitations of Walden’s rule. In particular, I will describe how the charge transport processes in this new class electrolytes achieve breakthrough performance by harnessing their unusual quantum character to drive the structural diffusion or Grotthuss diffusion mechanism, a phenomenon discoverable due to the power of machine learning. I will discuss the selection of candidate chemical species for each of the component processes, protocols for combining these components into a high-performance electrolyte, current results, and next steps in the evolution of the project. This work will serve to illustrate both the power of modern computational and machine learning approaches in the design of electrochemical systems but also to broaden the perspective on what constitutes a “breakthrough” electrolyte. Figure 1

Collective variable (CV) and generalized ensemble-based enhanced sampling methods are widely used for accelerating barrier-crossing events and enhancing conformational sampling in molecular dynamics simulations. Temperature-accelerated molecular dynamics (TAMD)/driven-adiabatic free energy dynamics (d-AFED) uses extended variables thermostated at high temperature to achieve better exploration of conformational space. Replica exchange with solute tempering (REST2) achieves improved sampling by scaling the solute-solute and solute-solvent interaction energies of different replicas and swapping conformations between them. It has been observed that a combination of CV-based enhanced sampling and global tempering is needed to boost the conformational sampling of large biomolecular systems due to the presence of large entropic basins. In this work, we propose a method called "Solute Tempered d-AFED" or "STed-AFED" that combines both d-AFED/TAMD and REST2. We implemented this approach in the OpenMM-UFEDMM interface and demonstrated the efficiency of this method by studying the conformational landscapes of small peptides and proteins, in particular, chignolin, Trp-cage, and villin.

ABT-333 and ABT-072 are two potent non-nucleoside NS5B polymerase inhibitors designed for the treatment of the hepatitis C virus (HCV). These structural analogs differ only by a minor substituent change, which disrupts the planarity of the naphthyl group on the ABT-333 compound through the addition of a more flexible trans-olefin substituent. However, this minor change leads to significant differences in their conformational preferences and intermolecular interactions, resulting in a ripple effect with drug development implications, ranging from crystal polymorphism and low aqueous solubility to formulation development challenges. In this article, we demonstrate how a suite of molecular simulation approaches, including crystal structure prediction augmented with a new hydrate CSP algorithm, free-energy perturbation, molecular dynamics (MD) based solubility predictions, and topological assessment to evaluate surface re-crystallization tendencies, provide key atomistic-level insights into the differentiated performance of the two analogs. Through this study, we establish the importance of end-to-end physics-based modeling, which involves explicit considerations of 3-D structure and crystal packing interactions. This approach provides structural and energetic insights into the physicochemical properties and drug development challenges faced when designing best-in-class drug molecules.

Accurate modeling of the dynamic structures of ribonucleic acid (RNA) molecules is essential for understanding their biological roles. However, such modeling remains challenging due to limitations in current force fields. This study critically evaluates three RNA force fields, HB-CUFIX, AMBER-χOL3, and AMBER-ROC, comparing their performance against experimental nuclear magnetic resonance and small-angle x-ray scattering data for single-stranded oligonucleotides. Using enhanced sampling techniques, specifically Unified Free Energy Dynamics, we exhaustively sampled the conformational space of tetramer and hexamer RNA sequences, achieving a detailed and thermodynamically converged view of their structural dynamics. Our findings reveal that HB-CUFIX outperforms AMBER-χOL3 and AMBER-ROC, providing near-experimental accuracy in capturing sequence-dependent structural preferences. In particular, HB-CUFIX accurately predicts low energy states for the AAAA and CCCC sequences, favoring A-form helical conformations, while the UUUU sequence adopts an extended, heterogeneous structure. The mixed GACC sequence displays a predominantly A-form helix with flexible terminal residues. These results highlight the significant role of sequence in dictating RNA conformational spaces, which are driven by base stacking interactions and covalent geometry. We also emphasize the importance of enhanced sampling, particularly methods that can handle large numbers of collective variables, in evaluating RNA force fields, as traditional brute-force Molecular dynamics fails to capture the conformational diversity of flexible RNAs. Our study provides a reliable tool for RNA structure prediction and dynamic analysis, supporting future advancements in RNA-targeted research and therapeutic design.

Representations are a foundational component of any modeling protocol, including on molecules and molecular solids. For tasks that depend on knowledge of both molecular conformation and 3D orientation, such as the modeling of molecular dimers, clusters, or condensed phases, we desire a rotatable representation that is provably complete in the types and positions of atomic nuclei and roto-inversion equivariant with respect to the input point cloud. In this paper, we develop, train, and evaluate a new type of autoencoder, molecular O(3) encoding net (Mo3ENet), for multi-type point clouds, for which we propose a new reconstruction loss, capitalizing on a Gaussian mixture representation of the input and output point clouds. Mo3ENet is end-to-end equivariant, meaning the learned representation can be manipulated on O(3), a practical bonus. An appropriately trained Mo3ENet latent space comprises a universal embedding for scalar, vector, and tensorial molecule property prediction tasks, as well as other downstream tasks incorporating the 3D molecular pose, and we demonstrate its fitness on several such tasks.

An approach is introduced to reduce the computational cost associated with performing path integral simulations. This work is an extension of the ring-polymer contraction (RPC) approach of Markland and Manolopoulos [J. Chem. Phys. 129, 024105 (2008)], which was originally formulated for simulating closed ring-polymers via molecular dynamics using normal modes to transform the full ring into a contracted ring. This work considers several new contraction transformation schemes, specifically a transformation to the contracted ring employing staging variables, designated as staging RPC, and a transformation for open-chain path integrals used in the calculation of momentum-dependent properties and quantum time correlation functions, which we designate as staging open-polymer contraction (OPC). These approaches are shown to be useful for obtaining equilibrium and dynamical properties for systems ranging from one-dimensional model problems to three-dimensional gas-phase water clusters. Additionally, we show that the accuracy and efficiency of contraction approaches can be improved significantly via a simple reweighting scheme. The advantage of staging RPC/OPC over normal mode RPC is the ability to formulate all transformations as a set of a recursion relations rather than matrix multiplications or Fourier transforms, providing a computationally simpler approach that is as effective as normal mode based contraction.

Hydrogen bonded electrolytes that exhibit accelerated proton transport via sequential reactive hops have drawn interest for their promise in clean energy applications. Molecular dynamics simulations of these electrolytes offer the opportunity to uncover microscopic mechanistic details that could be used to design and tune the properties of candidate electrolyte technologies. However, accurately modeling the proton transfer reactions and transport properties that give rise to high charge conductivites in these electrolytes proves computationally challenging because of the need to perform lengthy condensed phase simulations, treating both the electronic and nuclear degrees of freedom quantum mechanically. In this paper, we demonstrate that such a modeling task can be efficiently achieved with the use of density functional theory (DFT)-trained machine learning potentials (MLP) to accelerate path integral molecular dynamics (PIMD) simulations. We highlight the practical utility of this approach by using it to benchmark how closely PIMD simulations employing different DFT exchange-correlation functionals reproduce the composition-dependent densities, diffusion coefficients, and electrical conductivities of mixtures consisting of imidazole and levulinic acid. Even with the speedup afforded by our MLPs, PIMD simulations remain quite expensive. In order to render PIMD more computationally tractable, we introduce and benchmark the accuracy of a ring polymer contraction approach that leverages a computationally efficient short-range MLP to accelerate our PIMD simulations by an additional factor of four.

Because of its importance in various aspects of everyday life, silica is a material that has been the subject of extensive research. Studies on its amorphous phase have particularly benefited from the contribution of atomistic simulations to understand the close relationships between its structure and properties. In this context, the main difficulty lies in the compromise that had to be made between the precision of the interactions that need to be computed at an ab initio level and the important statistics required to describe disorder. With the advent of machine learning approaches, it is now possible to couple accuracy and statistics by using interatomic potentials trained on ab initio databases. This opens up unprecedented prospects for studies where calculation accuracy, system size and trajectory length are critical. In this work, we propose a machine learning potential for silica obtained from a neural network trained on a database consisting of a few hundred configurations extracted from an ab initio molecular dynamics trajectory at the Density Functional Theory (DFT) level of a silica liquid at high temperature and under pressure. We show that this potential is sufficiently accurate to describe the liquid and amorphous phases of silica, and that it is also transferable to glasses under moderate pressure and, more surprisingly, to certain crystalline phases.

Dynamic or structurally induced ionization is a critical aspect of many physical, chemical, and biological processes. Molecular dynamics (MD) based simulation approaches, specifically constant pH MD methods, have been developed to simulate ionization states of molecules or proteins under experimentally or physiologically relevant conditions. While such approaches are now widely utilized to predict ionization sites of macromolecules or to study physical or biological phenomena, they are often computationally expensive and require long simulation times to converge. In this article, using the principles of adiabatic free energy dynamics, we introduce an efficient technique for performing constant pH MD simulations within the framework of the adiabatic free energy dynamics (AFED) approach. We call the new approach pH-AFED. We show that pH-AFED provides highly accurate predictions of protein residue pK a values, with a MUE of 0.5 pK a units when coupled with driven adiabatic free energy dynamics (d-AFED), while reducing the required simulation times by more than an order of magnitude. In addition, pH-AFED can be easily integrated into most constant pH MD codes or implementations and flexibly adapted to work in conjunction with enhanced sampling algorithms that target collective variables. We demonstrate that our approaches, with both pH-AFED standalone as well as pH-AFED combined with collective variable based enhanced sampling, provide promising predictive accuracy, with a MUE of 0.6 and 0.5 pK a units respectively, on a diverse range of proteins and enzymes, ranging up to 186 residues and 21 titratable sites. Lastly, we demonstrate how this approach can be utilized to understand the in vivo performance engineered antibodies for immunotherapy.

Organic molecular crystals constitute a class of materials of critical importance in numerous industries. Despite the ubiquity of these systems, our ability to predict molecular crystal structures starting only from a two-dimensional diagram of the constituent compound(s) remains a significant challenge. Most structure-prediction protocols require a customized interatomic interaction model on which the quality of the results can depend sensitively. To overcome this problem, we introduce a new topological approach to molecular crystal structure prediction. The approach posits that in a stable structure, molecules are oriented such that principal axes and normal ring plane vectors are aligned with specific crystallographic directions and that heavy atoms occupy positions that correspond to minima of a set of geometric order parameters. By minimizing an objective function that encodes these orientations and atomic positions, and filtering based on the vdW free volume and intermolecular close contact distributions derived from the Cambridge Structural Database, stable structures and polymorphs for a given crystal can be predicted entirely mathematically without reliance on an interaction model. Reliable prediction of the organic molecular crystal structures is critical across numerous industries yet remains a significant challenge. Here, the authors develop a mathematical workflow based on topological concepts that reduces solution time to mere hours.

A force field as accurate as quantum mechanics (QMs) and as fast as molecular mechanics (MMs), with which one can simulate a biomolecular system efficiently enough and meaningfully enough to get quantitative insights, is among the most ardent dreams of biophysicists—a dream, nevertheless, not to be fulfilled any time soon. Machine learning force fields (MLFFs) represent a meaningful endeavor in this direction, where differentiable neural functions are parametrized to fit ab initio energies and forces through automatic differentiation. We argue that, as of now, the utility of the MLFF models is no longer bottlenecked by accuracy but primarily by their speed, as well as stability and generalizability—many recent variants, on limited chemical spaces, have long surpassed the chemical accuracy of 1 kcal/mol—the empirical threshold beyond which realistic chemical predictions are possible—though still magnitudes slower than MM. Hoping to kindle exploration and design of faster, albeit perhaps slightly less accurate MLFFs, in this review, we focus our attention on the technical design space (the speed-accuracy trade-off) between MM and ML force fields. After a brief review of the building blocks (from a machine learning-centric point of view) of force fields of either kind, we discuss the desired properties and challenges now faced by the force field development community, survey the efforts to make MM force fields more accurate and ML force fields faster, and envision what the next generation of MLFF might look like.

The dynamic interplay between a copper electrode and a reactive ionic liquid is illustrated in the cover picture, showcasing key components in the electrochemical reduction of CO2. The role of the cation in covering the electrode surface and the in-situ generated hydrogen bond donor from CO2 binding to the ionic liquid in enhancing the kinetics were uncovered by in-situ spectroscopy, complemented by electroanalytical and computational methods. Multi-carbon products with reduced reaction energy were obtained, as reported by Burcu Gurkan et al. in their Research Article (e202312163). Cover image credit: Miguel Munoz.