i-PI 3.0: A flexible and efficient framework for advanced atomistic simulations.

Atomic-scale simulations have progressed tremendously over the past decade, largely thanks to the availability of machine-learning interatomic potentials. These potentials combine the accuracy of electronic structure calculations with the ability to reach extensive length and time scales. The i-PI package facilitates integrating the latest developments in this field with advanced modeling techniques thanks to a modular software architecture based on inter-process communication through a socket interface. The choice of Python for implementation facilitates rapid prototyping but can add computational overhead. In this new release, we carefully benchmarked and optimized i-PI for several common simulation scenarios, making such overhead negligible when i-PI is used to model systems up to tens of thousands of atoms using widely adopted machine learning interatomic potentials, such as Behler-Parinello, DeePMD, and MACE neural networks. We also present the implementation of several new features, including an efficient algorithm to model bosonic and fermionic exchange, a framework for uncertainty quantification to be used in conjunction with machine-learning potentials, a communication infrastructure that allows for deeper integration with electronic-driven simulations, and an approach to simulate coupled photon-nuclear dynamics in optical or plasmonic cavities.

Melting entropy of crystals determined by electron-beam-induced configurational disordering.

Upon melting, the molecules in a crystal explore numerous configurations, reflecting an increase in disorder. The molar entropy of disorder can be defined by Boltzmann's formula ΔSd = Rln(Wd), where Wd is the increase in the number of microscopic states, so far inaccessible experimentally. We found that the Arrhenius frequency factor A of the electron diffraction signal decay provides Wd through an experimental equation A = AINTWd, where AINT is an inelastic scattering cross section. The method connects Clausius and Boltzmann experimentally and supplements the Clausius approach, being applicable to a femtogram quantity of thermally unstable and biomolecular crystals. The data also showed that crystal disordering and crystallization of melt are reciprocal, both governed by the entropy change but manifesting in opposite directions.

The effect of ligands on the size distribution of copper nanoclusters: Insights from molecular dynamics simulations.

Controlling the size distribution in the nucleation of copper particles is crucial for achieving nanocrystals with desired physical and chemical properties. However, their synthesis involves a complex system of solvents, ligands, and copper precursors with intertwining effects on the size of the nanoclusters. We combine molecular dynamics simulations and density functional theory calculations to provide insights into the nucleation mechanism in the presence of a triphenyl phosphite ligand. We identify the crucial role of the strength of the metal-phosphine interaction in inhibiting the cluster's growth. We demonstrate computationally several practical routes to fine-tune the interaction strength by modifying the side groups of the additive. Our work provides molecular insights into the complex nucleation process of protected copper nanocrystals, which can assist in controlling their size distribution and, eventually, their morphology.

Short-Time Infrequent Metadynamics for Improved Kinetics Inference.

Infrequent Metadynamics is a popular method to obtain the rates of long time-scale processes from accelerated simulations. The inference procedure is based on rescaling the first-passage times of the Metadynamics trajectories using a bias-dependent acceleration factor. While useful in many cases, it is limited to Poisson kinetics, and a reliable estimation of the unbiased rate requires slow bias deposition and prior knowledge of efficient collective variables. Here, we propose an improved inference scheme, which is based on two key observations: (1) the time-independent rate of Poisson processes can be estimated using short trajectories only. (2) Short trajectories experience minimal bias, and their rescaled first-passage times follow the unbiased distribution even for relatively high deposition rates and suboptimal collective variables. Therefore, by basing the inference procedure on short time scales, we obtain an improved trade-off between speedup and accuracy at no additional computational cost, especially when employing suboptimal collective variables. We demonstrate the improved inference scheme for a model system and two molecular systems.

Combining stochastic resetting with Metadynamics to speed-up molecular dynamics simulations.

Metadynamics is a powerful method to accelerate molecular dynamics simulations, but its efficiency critically depends on the identification of collective variables that capture the slow modes of the process. Unfortunately, collective variables are usually not known a priori and finding them can be very challenging. We recently presented a collective variables-free approach to enhanced sampling using stochastic resetting. Here, we combine the two methods, showing that it can lead to greater acceleration than either of them separately. We also demonstrate that resetting Metadynamics simulations performed with suboptimal collective variables can lead to speedups comparable with those obtained with optimal collective variables. Therefore, applying stochastic resetting can be an alternative to the challenging task of improving suboptimal collective variables, at almost no additional computational cost. Finally, we propose a method to extract unbiased mean first-passage times from Metadynamics simulations with resetting, resulting in an improved tradeoff between speedup and accuracy. This work enables combining stochastic resetting with other enhanced sampling methods to accelerate a broad range of molecular simulations.

Fermionic physics from ab initio path integral Monte Carlo simulations of fictitious identical particles.

The ab initio path integral Monte Carlo (PIMC) method is one of the most successful methods in statistical physics, quantum chemistry and related fields, but its application to quantum degenerate Fermi systems is severely hampered by an exponential computational bottleneck: the notorious fermion sign problem. Very recently, Xiong and Xiong [J. Chem. Phys. 157, 094112 (2022)] have suggested to partially circumvent the sign problem by carrying out simulations of fictitious systems guided by an interpolating continuous variable ξ ∈ [-1, 1], with the physical Fermi- and Bose-statistics corresponding to ξ = -1 and ξ = 1. It has been proposed that information about the fermionic limit might be obtained by calculations within the bosonic sector ξ > 0 combined with an extrapolation throughout the fermionic sector ξ < 0, essentially bypassing the sign problem. Here, we show how the inclusion of the artificial parameter ξ can be interpreted as an effective penalty on the formation of permutation cycles in the PIMC simulation. We demonstrate that the proposed extrapolation method breaks down for moderate to high quantum degeneracy. Instead, the method constitutes a valuable tool for the description of large Fermi-systems of weak quantum degeneracy. This is demonstrated for electrons in a 2D harmonic trap and for the uniform electron gas (UEG), where we find excellent agreement (∼0.5%) with exact configuration PIMC results in the high-density regime while attaining a speed-up exceeding 11 orders of magnitude. Finally, we extend the idea beyond the energy and analyze the radial density distribution (2D trap), as well as the static structure factor and imaginary-time density-density correlation function (UEG).

Quadratic scaling bosonic path integral molecular dynamics.

Bosonic exchange symmetry leads to fascinating quantum phenomena, from exciton condensation in quantum materials to the superfluidity of liquid 4He. Unfortunately, path integral molecular dynamics (PIMD) simulations of bosons are computationally prohibitive beyond ∼100 particles, due to a cubic scaling with the system size. We present an algorithm that reduces the complexity from cubic to quadratic, allowing the first simulations of thousands of bosons using PIMD. Our method is orders of magnitude faster, with a speedup that scales linearly with the number of particles and the number of imaginary time slices (beads). Simulations that would have otherwise taken decades can now be done in days. In practice, the new algorithm eliminates most of the added computational cost of including bosonic exchange effects, making them almost as accessible as PIMD simulations of distinguishable particles.

Collective Variables for Conformational Polymorphism in Molecular Crystals.

Controlling polymorphism in molecular crystals is crucial in the pharmaceutical, dye, and pesticide industries. However, its theoretical description is extremely challenging, due to the associated long time scales (>1 µs). We present an efficient procedure for identifying collective variables that promote transitions between conformational polymorphs in molecular dynamics simulations. It involves applying a simple dimensionality reduction algorithm to data from short (∼ps) simulations of the isolated conformers that correspond to each polymorph. We demonstrate the utility of our method in the challenging case of the important energetic material, CL-20, which has three anhydrous conformational polymorphs at ambient pressure. Using these collective variables in Metadynamics simulations, we observe transitions between all solid polymorphs in the biased trajectories. We reconstruct the free energy surface and identify previously unknown defect and intermediate forms in the transition from one known polymorph to another. Our method provides insights into complex conformational polymorphic transitions of flexible molecular crystals.

Stochastic Resetting for Enhanced Sampling.

We present a method for enhanced sampling of molecular dynamics simulations using stochastic resetting. Various phenomena, ranging from crystal nucleation to protein folding, occur on time scales that are unreachable in standard simulations. They are often characterized by broad transition time distributions, in which extremely slow events have a non-negligible probability. Stochastic resetting, i.e., restarting simulations at random times, was recently shown to significantly expedite processes that follow such distributions. Here, we employ resetting for enhanced sampling of molecular simulations for the first time. We show that it accelerates long time scale processes by up to an order of magnitude in examples ranging from simple models to a molecular system. Most importantly, we recover the mean transition time without resetting, which is typically too long to be sampled directly, from accelerated simulations at a single restart rate. Stochastic resetting can be used as a standalone method or combined with other sampling algorithms to further accelerate simulations.

Prediction of a Supersolid Phase in High-Pressure Deuterium.

Supersolid is a mysterious and puzzling state of matter whose possible existence has stirred a vigorous debate among physicists for over 60 years. Its elusive nature stems from the coexistence of two seemingly contradicting properties, long-range order and superfluidity. We report computational evidence of a supersolid phase of deuterium under high pressure (p>800 GPa) and low temperature (T<1.0 K). In our simulations, that are based on bosonic path integral molecular dynamics, we observe a highly concerted exchange of atoms while the system preserves its crystalline order. The exchange processes are favored by the soft core interactions between deuterium atoms that form a densely packed metallic solid. At the zero temperature limit, Bose-Einstein condensation is observed as the permutation probability of N deuterium atoms approaches 1/N with a finite superfluid fraction. Our study provides concrete evidence for the existence of a supersolid phase in high-pressure deuterium and could provide insights on the future investigation of supersolid phases in real materials.

Attenuating the fermion sign problem in path integral Monte Carlo simulations using the Bogoliubov inequality and thermodynamic integration.

Accurate thermodynamic simulations of correlated fermions using path integral Monte Carlo (PIMC) methods are of paramount importance for many applications such as the description of ultracold atoms, electrons in quantum dots, and warm-dense matter. The main obstacle is the fermion sign problem (FSP), which leads to an exponential increase in computation time both with an increase in the system size and with a decrease in the temperature. Very recently, Hirshberg et al. [J. Chem. Phys. 152, 171102 (2020)] have proposed to alleviate the FSP based on the Bogoliubov inequality. In the present work, we extend this approach by adding a parameter that controls the perturbation, allowing for an extrapolation to the exact result. In this way, we can also use thermodynamic integration to obtain an improved estimate of the fermionic energy. As a test system, we choose electrons in 2D and 3D quantum dots and find in some cases a speed-up exceeding 106, as compared to standard PIMC, while retaining a relative accuracy of ∼0.1%. Our approach is quite general and can readily be adapted to other simulation methods.

Metadynamics of Paths.

We present a method to sample reactive pathways via biased molecular dynamics simulations in trajectory space. We show that the use of enhanced sampling techniques enables unconstrained exploration of multiple reaction routes. Time correlation functions are conveniently computed via reweighted averages along a single trajectory and kinetic rates are accessed at no additional cost. These abilities are illustrated analyzing a model potential and the umbrella inversion of NH_{3} in water. The algorithm allows a parallel implementation and promises to be a powerful tool for the study of rare events.

Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality.

We present a method for performing path integral molecular dynamics (PIMD) simulations for fermions and address its sign problem. PIMD simulations are widely used for studying many-body quantum systems at thermal equilibrium. However, they assume that the particles are distinguishable and neglect bosonic and fermionic exchange effects. Interacting fermions play a key role in many chemical and physical systems, such as electrons in quantum dots and ultracold trapped atoms. A direct sampling of the fermionic partition function is impossible using PIMD since its integrand is not positive definite. We show that PIMD simulations for fermions are feasible by employing our recently developed method for bosonic PIMD and reweighting the results to obtain fermionic expectation values. The approach is tested against path integral Monte Carlo (PIMC) simulations for up to seven electrons in a two-dimensional quantum dot for a range of interaction strengths. However, like PIMC, the method suffers from the sign problem at low temperatures. We propose a simple approach for alleviating it by simulating an auxiliary system with a larger average sign and obtaining an upper bound to the energy of the original system using the Bogoliubov inequality. This allows fermions to be studied at temperatures lower than would otherwise have been feasible using PIMD, as demonstrated in the case of a three-electron quantum dot. Our results extend the boundaries of PIMD simulations of fermions and will hopefully stimulate the development of new approaches for tackling the sign problem.

Path integral molecular dynamics for bosons.

Trapped bosons exhibit fundamental physical phenomena and are at the core of emerging quantum technologies. We present a method for simulating bosons using path integral molecular dynamics. The main difficulty in performing such simulations is enumerating all ring-polymer configurations, which arise due to permutations of identical particles. We show that the potential and forces at each time step can be evaluated by using a recurrence relation which avoids enumerating all permutations, while providing the correct thermal expectation values. The resulting algorithm scales cubically with system size. The method is tested and applied to bosons in a 2-dimensional (2D) trap and agrees with analytical results and numerical diagonalization of the many-body Hamiltonian. An analysis of the role of exchange effects at different temperatures, through the relative probability of different ring-polymer configurations, is also presented.

Hydrogenic Stretch Spectroscopy of Glycine-Water Complexes: Anharmonic Ab Initio Classical Separable Potential Calculations.

The anharmonic frequencies of O-H, C-H, and N-H stretching modes of hydrogen-bonded glycine-H2O complexes are calculated using ab initio classical separable potential approximation. In this approach, ab initio molecular dynamic simulations are used to determine an effective classical potential for each of the normal modes of the system. The frequencies are calculated by solving the time-independent Schrödinger equation for each mode using time-averaged potentials. Three complex structures are studied, which differ in the location of the water molecule on the amino acid. Significant differences are found between the spectra of the three structures, and signatures of individual complexes are established. It is demonstrated that anharmonic effects are essential in the discrimination between different structures, while frequency differences at the harmonic level are much smaller. Intensities are also computed and found to carry information on differences between structures, but the role of anharmonicity in this is small.

N2O5 at water surfaces: binding forces, charge separation, energy accommodation and atmospheric implications.

Interactions of N2O5 with water media are of great importance in atmospheric chemistry and have been the topic of extensive research for over two decades. Nevertheless, many physical and chemical properties of N2O5 at the surface or in bulk water are unknown or not microscopically understood. This study presents extensive new results on the physical properties of N2O5 in water and at the surface of water, with a focus on their microscopic basis. The main results are obtained using ab initio molecular dynamics and calculations of a potential of mean force. These include: (1) collisions of N2O5 with water at 300 K lead to trapping at the surface for at least 20 ps and with 95% probability. (2) During that time, there is no N2O5 hydrolysis, evaporation, or entry into the bulk. (3) Charge separation between the NO2 and NO3 groups of N2O5, fluctuates significantly with time. (4) Energy accommodation of the colliding N2O5 at the surface takes place within picoseconds. (5) The binding energy of N2O5 to a nanosize amorphous ice particle at 0 K is on the order of 15 kcal mol-1 for the main surface site. N2O5 binding to the cluster is due to one weak hydrogen bond and to interactions between partial charges on the N2O5 and on water. (6) The free-energy profile was calculated for transporting N2O5 from the gas phase through the interface and into bulk water. The corresponding concentration profile exhibits a propensity for N2O5 at the aqueous surface. The free energy barrier for entry from the surface into the bulk was determined to be 1.8 kcal mol-1. These findings are used to interpret recent experiments. We conclude with implications of this study for atmospheric chemistry.

Trapping and Structural Characterization of the XNO2·NO3- (X = Cl, Br, I) Exit Channel Complexes in the Water-Mediated X- + N2O5 Reactions with Cryogenic Vibrational Spectroscopy.

The heterogeneous reaction of N2O5 with sea spray aerosols yields the ClNO2 molecule, which is postulated to occur through water-mediated charge separation into NO3- and NO2+ ions followed by association with Cl-. Here we address an alternative mechanism where the attack by a halide ion can yield XNO2 by direct insertion in the presence of water. This was accomplished by reacting X-(D2O)n (X = Cl, Br, I) cluster ions with N2O5 to produce ions with stoichiometry [XN2O5]-. These species were cooled in a 20 K ion trap and structurally characterized by vibrational spectroscopy using the D2 messenger tagging technique. Analysis of the resulting band patterns with DFT calculations indicates that they all correspond to exit channel ion-molecule complexes based on the association of NO3- with XNO2, with the NO3- constituent increasingly perturbed in the order I > Br > Cl. These results establish that XNO2 can be generated even when more exoergic reaction pathways involving hydrolysis are available and demonstrate the role of the intermediate [XN2O5]- in the formation of XNO2.

Approximate Quantum Dynamics using Ab Initio Classical Separable Potentials: Spectroscopic Applications.

Algorithms for quantum molecular dynamics simulations that directly use ab initio methods have many potential applications. In this article, the ab initio classical separable potentials (AICSP) method is proposed as the basis for approximate algorithms of this type. The AICSP method assumes separability of the total time-dependent wave function of the nuclei and employs mean-field potentials that govern the dynamics of each degree of freedom. In the proposed approach, the mean-field potentials are determined by classical ab initio molecular dynamics simulations. The nuclear wave function can thus be propagated in time using the effective potentials generated "on the fly". As a test of the method for realistic systems, calculations of the stationary anharmonic frequencies of hydrogen stretching modes were carried out for several polyatomic systems, including three amino acids and the guanine-cytosine pair of nucleobases. Good agreement with experiments was found. The method scales very favorably with the number of vibrational modes and should be applicable for very large molecules, e.g., peptides. The method should also be applicable for properties such as vibrational line widths and line shapes. Work in these directions is underway.

Formation of Carbonic Acid in Impact of CO2 on Ice and Water.

A new mode of formation is proposed for carbonic acid in the atmosphere. It involves impact of vibrationally excited gas-phase CO2 molecules on water or ice particles. This is a first mechanism that supports formation on ice as well as on liquid water surfaces. Results of ab initio molecular dynamics simulations are presented on collisions of CO2 with (H2O)n clusters (n = 1, 4, 8, 12). Efficient formation of carbonic acid is seen with product lifetimes exceeding 100 ps. The reaction is feasible even for collision of CO2 with a single water molecule but in a different mechanism than for larger clusters. For clusters, the transition state shows charge separation into H3O(+)···HCO3(-), which transforms into neutral carbonic acid as the product, hydrated by the remaining waters. Possible atmospheric implications of the results are discussed.

Decomposition of condensed phase energetic materials: interplay between uni- and bimolecular mechanisms.

Activation energy for the decomposition of explosives is a crucial parameter of performance. The dramatic suppression of activation energy in condensed phase decomposition of nitroaromatic explosives has been an unresolved issue for over a decade. We rationalize the reduction in activation energy as a result of a mechanistic change from unimolecular decomposition in the gas phase to a series of radical bimolecular reactions in the condensed phase. This is in contrast to other classes of explosives, such as nitramines and nitrate esters, whose decomposition proceeds via unimolecular reactions both in the gas and in the condensed phase. The thermal decomposition of a model nitroaromatic explosive, 2,4,6-trinitrotoluene (TNT), is presented as a prime example. Electronic structure and reactive molecular dynamics (ReaxFF-lg) calculations enable to directly probe the condensed phase chemistry under extreme conditions of temperature and pressure, identifying the key bimolecular radical reactions responsible for the low activation route. This study elucidates the origin of the difference between the activation energies in the gas phase (~62 kcal/mol) and the condensed phase (~35 kcal/mol) of TNT and identifies the corresponding universal principle. On the basis of these findings, the different reactivities of nitro-based organic explosives are rationalized as an interplay between uni- and bimolecular processes.