Browse Articles

Superionicity in ammonium polyhydrides at extreme pressures

Polyhydrides have been shown to form novel structures at high pressure, which may be found in the interiors of giant planets. With density functional molecular dynamics simulations, we studied the behavior of ammonium polyhydride compounds with stoichiometries of NH7, NH9, NH10, NH11, NH14, NH20, and NH24, which were predicted with crystal structure search methods to be metastable at 100–300 GPa. For every compound, we performed simulations at a range of temperatures (and for several compounds, pressures) covering the solid, superionic, and liquid phases. We show that when heated, high pressure ammonium polyhydride compounds exhibit hydrogen superionic diffusion. We demonstrate a number of metrics by which the solid-to-superionic and superionic-to-liquid transitions can be detected from simulation data, including changes in the internal energy and pressure, formation of new chemical species, and atomic diffusion rates. We find that both the solid-to-superionic and the superionic-to-liquid transitions decrease in temperature as proton fraction increases. These trends indicate that above a proton fraction of ∼0.97, ammonium hydride structures are likely to directly melt instead of first exhibiting a superionic phase. Our observed melting trend further indicates that at the extreme conditions of ice giant interiors, hydrogen rich ammonium hydrides such as those studied in this work would exist predominantly as liquids rather than exhibiting a superionic phase.

Erratum: “Interactions-driven structural morphology and motion in two-dimensional active Brownian systems” [J. Chem. Phys. 163, 194902 (2025)]

A patchy-particle three-dimensional octagonal quasicrystal

We devise an ideal three-dimensional octagonal quasicrystal that is based upon the two-dimensional Ammann–Beenker tiling and that is potentially suitable for realization with patchy particles. Based on an analysis of its local environments, we design a binary system of five- and eight-patch particles that in simulations assembles into a three-dimensional octagonal quasicrystal. The local structure is subtly different from the original ideal quasicrystal possessing a narrower coordination-number distribution; in fact, the eight-patch particles are not needed and a one-component system of the five-patch particles assembles into an essentially identical octagonal quasicrystal. We also consider a one-component system of the eight-patch particles; this assembles into a cluster with a number of crystalline domains, which, because of the coherent boundaries between the crystallites, has approximate eightfold order. We envisage that these systems could be realized using DNA origami or protein design.

Analytic polarizability and vibrational Raman spectra from constrained nuclear-electronic orbital density functional theory

Accurate and efficient simulation of vibrational Raman spectra for systems with strong anharmonicity and nuclear quantum effects remains challenging. Herein, we apply the recently developed constrained nuclear-electronic orbital (CNEO) framework to simulate Raman spectra. We implement analytic static polarizabilities within CNEO density functional theory (CNEO-DFT) and compute Raman spectra using both CNEO harmonic analysis and CNEO molecular dynamics (CNEO-MD). In the harmonic analysis approach, vibrational frequencies are obtained by diagonalizing the mass-weighted CNEO Hessian, and Raman intensities are derived from the projected derivatives of the polarizability with respect to normal-mode coordinates. In CNEO-MD, both frequency and intensity information are extracted from the Fourier transform of the projected polarizability time-derivative autocorrelation function. Applications to formic acid and the mixed water dimer/trimer system show that, compared to conventional DFT, CNEO-DFT achieves substantially improved accuracy, particularly for modes with significant hydrogen motion. Overall, the CNEO framework provides an accurate and efficient approach for simulating vibrational Raman spectra and is especially promising for systems in which hydrogen motion plays a critical role.

Second harmonic generation observed by dielectric relaxation at high electric fields: SHG without optics

We have measured second harmonic generation via the nonlinear dielectric permittivity of a polar glass-forming liquid, propylene glycol, with the inversion symmetry of the liquid broken by applying a dc bias field. For a given combined peak field, EB + E0, the highest second harmonic signals are obtained when the dc bias field and peak ac-field have the same amplitude, i.e., EB = E0. Second harmonic results measured in the static limit agree well with the theoretical prediction based upon the third-order nonlinear susceptibility term that connects polarization P with the static field E. Second harmonic signals are detectable at all frequencies at which dipole orientation contributes to permittivity, possibly occurring the instant a dc bias field is applied. The results imply that second harmonic generation as a signature of anisotropy, e.g., at interfaces, can be assessed by impedance spectroscopy at twice the fundamental frequency, not only by optical techniques.

Representation of electron–nucleus cusps in Slater and Gaussian basis sets

Prior to spherical averaging, the Coulomb cusp of exact electron densities is described by an equation involving two parameters: the nuclear charge Z and a vector quantifying the cusp anisotropy due to the presence of other nuclei. A similar three-parameter equation describes jump discontinuities of exact kinetic energy densities at nuclear positions. We demonstrate that these same equations, with basis-set-dependent parameters, also describe approximate electron and kinetic energy densities expanded in Slater-type orbitals (STOs) and even cuspless Gaussian-type orbitals (GTOs). Apart from the effective Z values, which are zero for GTOs, the cusp-shape parameters derived from STO and GTO wavefunctions are comparable in magnitude and converge to common complete-basis-set limits. These findings provide a rigorous justification for extending the concepts of effective nuclear charge and density anisotropy vectors to approximate wavefunctions within STO and GTO basis sets. Such extensions are not entailed by the properties of exact solutions to the Schrödinger equation and are not automatically applicable to all basis sets.

Interfacial ion hydration and electrostatics govern salt precipitation and crystal morphology

Salt precipitation, a phenomenon central to processes such as soil salinization, water treatment, and energy storage, is dictated by nanoscale hydration and electrostatic interactions at crystal interfaces. Understanding how their interplay gives rise to salt-specific growth behaviors and morphologies remains a critical challenge. Combining large-scale molecular dynamics simulations with free energy calculations, we investigate the precipitation behavior of three common salts, namely, NaCl, KCl, and Na2SO4, on their respective crystal surfaces. While all crystals grow steadily under controlled supersaturation, their surface morphologies differ markedly. KCl grows in a nearly ideal layer-by-layer mode, whereas NaCl and Na2SO4 develop increasingly rough, defect-rich surfaces. This roughening is accompanied by an early-stage charge imbalance at the crystal surface, characterized by preferential cation adsorption. The resulting local electrostatic environment modifies the structuring and dynamics of interfacial water, which in turn facilitates the premature growth of subsequent layers. We further find that ion adsorption is strongly defect-dependent: kinks provide stabilizing environments, whereas steps are energetically unfavorable due to their distinct hydration. Our study offers molecular-level insight into the coupled roles of ion-specific adsorption and interfacial hydration in shaping crystal morphology for different salts.

Stochastic representation of time-evolving neural network-based wavefunctions

Solving the time-dependent Schrödinger equation (TDSE) is pivotal for modeling nonadiabatic electron dynamics, a key process in ultrafast spectroscopy and laser–matter interactions. However, exact solutions to the TDSE remain computationally prohibitive for most realistic systems, as the Hilbert space expands exponentially with dimensionality. In this work, we propose an approach integrating the stochastic representation framework with a neural network wavefunction ansatz, a flexible model capable of approximating time-evolving quantum wavefunctions. We first validate the method on one-dimensional, single-electron systems, focusing on ionization dynamics under intense laser fields, a critical process in attosecond physics. Our results demonstrate that the approach accurately reproduces key features of quantum evolution, including the energy and dipole evolution during ionization. We further show the feasibility of extending this approach to three-dimensional systems. Due to the increased complexity of real-time simulations in higher dimensions, these results remain at an early stage and highlight the need for more advanced stabilization strategies.

Chirality, confinement, and dimensionality govern re-entrant transitions in active matter

The non-equilibrium dynamics of individual chiral active particles underpin the complex behavior of chiral active matter. Here, we present an exact analytical framework, supported by simulations, to characterize the steady states of two-dimensional chiral active Brownian particles and three-dimensional torque-driven counterparts in a harmonic trap. Using a Laplace-transform approach to the Fokker–Planck equation, we derive closed-form expressions for displacement moments and excess kurtosis, providing a precise probe of non-Gaussian statistics. Our analysis reveals three distinct regimes characterized by bimodal distributions with off-center peaks, Gaussian-like distributions, and weakly heavy-tailed distributions unique to two dimensions. We show that dimensionality plays a decisive role: in two dimensions, increasing chirality suppresses activity and restores Gaussian-like behavior, whereas in three dimensions, torque sustains anisotropic steady states and preserves non-Gaussianity even at high chirality. These behaviors are captured by simple active length-scale arguments that map the boundaries between Gaussian-like and non-Gaussian states. Our results offer concrete experimental signatures—including kurtosis crossovers, off-center peaks, and torque-induced anisotropy—that establish confinement as a powerful tool to probe and control chiral and torque-driven active matter.

Quantum one-dimensional materials at low temperature: Pre-critical effects, long-range correlations, and specific heats

In 1925, Ising showed that a one-dimensional (1D) system of particles that interact with short-range intermolecular forces cannot exhibit a phase transition at finite temperatures. However, he also found that the correlation length ξ diverges to infinity at T = 0 K and low density, the signature of a critical point. In this work, we report classical and quantum mechanical path integral Monte Carlo (PIMC) results for the static pair correlation function, the correlation length, and specific heats for a 1D system of spin 0 Lennard-Jones (LJ) molecules modeled on hydrogen, H2. Our PIMC results, which are for temperatures down to T* = kBT/ɛ = 0.01 and densities down to ρ* = ρσ = 0.1, exhibit strong long-range fluctuations out to 50 molecular diameters or more at low temperatures and low densities. In addition, the correlation length appears to diverge as the temperature approaches 0 K, thus strongly confirming Ising’s prediction of a critical point at T = 0 K. Our classical Monte Carlo results, by contrast, show smaller values of the correlation length, since it does not capture the quantum dispersion effect. The specific heat at constant length is also found to diverge on lowering the temperature toward absolute zero, a further signature of the approach to a critical point, but the third law of thermodynamics requires that it drop to zero at a temperature that is below the lowest value studied in this work. The findings reported here are relevant to recent advances in the synthesis of one-dimensional van der Waals materials.

On entropy driven spin crossover in pristine and fluorinated Fe(phenanthroline)2(NCS)2: Insight from DFT

An ab initio total energy calculation of the Fe(phenanthroline)2(NCS)2 (Fephen) molecule adsorbed on a copper surface revealed that substituting hydrogen with fluorine in the spin-crossover Fephen molecule drastically alters its electronic structure, leading to an inversion of its magnetic state. To understand the temperature-dependent physics of spin crossover, various vibrational modes were computed within the quasi-harmonic approximation using a finite-difference method. The analysis showed that the primary frequency differences between the high-spin (HS) and low-spin states are concentrated in the N–C–S and Fe–ligand stretching modes for both the pristine and fluorinated molecules. This observation aligns well with expectations, as these specific bonds are directly involved and most impacted by the spin-crossover phenomenon. Free energy calculations indicated that the pristine Fephen molecule undergoes a spin crossover at 167 K, which closely matches the experimental value of 175 K. In contrast, the fluorinated Fephen molecule does not exhibit any spin crossover, with its HS state remaining the most stable across all temperatures. These calculated frequencies were also instrumental in determining the evolution of the HS fraction with respect to temperature in the pristine molecule, utilizing the Slichter–Drickamer model. The results from this model showed good agreement with experimental observations, further validating this approach.

From soup to structure: Simulating hydrated semi-crystalline proton exchange membranes

A methodology is presented for the simulation of complex high molecular weight polymers, with a primary focus on reproducing the structure of Nafion® proton exchange membranes. Enhanced computational efficiency is achieved in comparison to commonly employed Monte Carlo techniques by implementing random insertions of long polymer chains followed by a ghost chain randomization process. Dissipative particle dynamics is applied to relax strongly overlapped configurations and to access the long time scales necessary to capture the crystallization process. The developed protocol first relaxes soft polymer chains in a stage analogous to thermal activation, which promotes crystallinity, and subsequently incorporates chain stiffening to reproduce crystal growth during cooling. The approach is validated against experimental SAXS measurements by computing intensity profiles from the electronic density of simulated Nafion. The results highlight the strong influence of side chain distribution on crystallinity, emphasizing its role in the formation of realistic semi-crystalline morphologies. Key processes, including chain ordering, local alignment, and molecular packing, are resolved, providing an improved understanding of structure–property relationships. The methodology accurately predicts crystallinity content as well as crystallite size and shape, thus demonstrating strong predictive capability for a wide range of PEMs and related polymer systems.

The liquidus curve of aqueous methanol mixtures. The high-alcohol concentration region at standard pressure

This work extends our previous simulation study on the water/methanol liquidus curve at low methanol concentrations [M. Martínez-Jiménez and A. L. Benavides, J. Chem. Phys. 157, 104502 (2022)]. We now explored methanol-rich mixtures that are located beyond the eutectic point, at p = 0.1 MPa. To our knowledge, the liquidus temperatures for these mixtures have not been obtained before with any simulation technique. In this case, the liquidus temperature curve characterizes the boundary between the liquid methanol/water mixture and a coexistence of a liquid mixture with the β-methanol solid. Direct-coexistence simulations were performed using the combination of the TIP4P/Ice model for water and the OPLS/2016 model for methanol. The obtained liquidus temperatures follow the experimental data trend, although a systematic overestimation is observed. Knowing that this deviation is possibly due to the OPLS/2016 model, which overestimates the methanol melting point, we hope that the information presented in this work will help refine this model, thereby improving the force field for mixtures. In addition, we present the methanol freezing-point depression by water that could be another target property to consider in the design of better models.

Erratum: “A generalized method for refining and selecting random crystal structures using graph theory” [J. Chem. Phys. 163, 094106 (2025)]

Water induced redistribution of polarons in rutile TiO2(011)-(2 × 1)

The interaction between water and TiO2 is a widely discussed topic due to its direct relevance to green hydrogen production. In addition to the impact of TiO2 substrate structure on the adsorption of water, adsorbates have also been found to affect the distribution of excess electrons, which form polaron-like bandgap states (BGS) in the substrate for rutile TiO2(110) and anatase TiO2(101). In order to further investigate the adsorbate induced redistribution of polarons in TiO2, we measured water coverage-dependent BGS and surface species on rutile TiO2(011) [R-TiO2(011)], another model metal oxide surface, using ultraviolet and x-ray photoelectron spectroscopy (UPS and XPS). We found the increase in BGS magnitude and surface hydroxyl quantity with increasing water coverage below 1 monolayer (ML). Ab initio molecular dynamics simulations suggested that excess electrons in the bulk of R-TiO2(011) are abstracted to the surface and induce water dissociation. Our results expand the applicable systems for the adsorbate induced polaron redistribution and indicate that it might be a general phenomenon. Knowledge about the adsorbate–polaron interaction in semiconductors may contribute to the understanding and optimization of surface processes involving molecule–substrate interaction, such as photocatalysis and photovoltaics.

Complex-energy second-order approximate coupled-cluster methods for electronic resonances

Electronic resonances are metastable states with finite lifetimes, encountered in processes such as photodetachment, electron transmission spectroscopy, and Auger decay. Resonances appear in Hermitian quantum mechanics as increased density of states in the continuum rather than as discrete energy levels. To describe resonances accurately, including their coupling to the continuum, methods based on non-Hermitian quantum mechanics can be used, which yield complex energies. In this study, we combine the complex absorbing potential and complex basis functions techniques with the RI-CC2 method. The second-order coupled cluster method (CC2) offers a good balance between accuracy and computational cost by approximating equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) theory, making it suitable for studying electronic resonances in larger molecules. The resolution-of-the-identity (RI) approximation further reduces computational demands without significant loss in accuracy. We investigate the numerical performance of the new complex-energy RI-CC2 methods, focusing on temporary anions. Negative electron affinities and decay widths can be computed using the electron-attachment (EA) variant of RI-CC2. For N2, C2H4, CH2O, and HCOOH, EA-CC2 yields electron affinities about 0.1–0.2 eV smaller than EOM-EA-CCSD, while deviations reach 0.5 eV for larger anions such as uracil, naphthalene, cyanonaphthalene, and pyrene. As a result of these trends, EA-CC2 is in better agreement with experiment for the negative electron affinities than EOM-EA-CCSD for all studied anions. The corresponding resonance widths from EA-CC2 calculations are about 0.05–0.25 eV smaller compared to EOM-EA-CCSD. Semi-empirical spin-scaling increases electron affinities by 0.3–0.5 eV and broadens resonance widths, improving the agreement with EOM-EA-CCSD but worsening the agreement with experiment.

A new reaction coordinate to study the translocation pathway of cell-penetrating peptides across lipid bilayers: The cases of transportan-10 and penetratin

Translocation pathway of cell-penetrating peptides remains elusive, as it is hard to observe by experimental and theoretical studies, which limits their effective use. Furthermore, lipid dynamics influence the translocation pathway, which is often overlooked due to its slow timescale. Current studies lack the effect of multiple peptides on the translocation process. Therefore, in this work, we employ the umbrella sampling technique with a preferential lipid–peptide interaction term in the reaction coordinate to explore the translocation activity of penetratin and transportan-10 (TP10) peptides in a heterogeneous membrane. In experiments, they follow different pathways according to their concentration, but the cause of this difference is unknown. We considered single and multiple (two and four) peptide translocation processes to understand the differences. Self-aggregation process is taken into account for multiple peptides. The interaction between peptides and peptides–lipids is found to be important for a proper overview of the translocation process. Peptide translocation was found to be related to the dynamics of the lipids, which change during the translocation process, making the system complex to study. In the case of multiple penetratin translocation, the anionic lipids were found to aggregate on the positive curvature of the upper leaflet, helping fold the membrane. Lipid composition of the TP10 multiple peptide case was found random. The increased mass and size of the solute in this case helped attain a radius more than the threshold value, leading to pore formation. Free energy barriers of single TP10 and penetratin are found to be 45.4 ± 2 and 33.7 ± 0.8 kJ mol−1, respectively.

State-to-state dynamics of P(2D) + H2(X1Σg+) reactions based on the new neural network potential energy surface

To investigate the state-to-state dynamics of the P(D2) + H2(XΣg+1)(v = 0, j = 0) reaction, we reconstructed the PH2(X2B1) potential energy surface (PES) using the permutation invariant polynomial neural network (PIP-NN) method based on 40 595 ab initio points. The aug-cc-pVQZ basis sets with Davidson correction were employed throughout the calculations. The reference wave function for the multi-reference configuration interaction calculations was constructed from a full valence complete-active-space self-consistent field wave function. To achieve higher PES accuracy, the double many-body expansion–scaled external correlation (DMBE–SEC) method was applied to extrapolate to the one-electron complete basis set limit, yielding a total root-mean square deviation of 2.4 meV for the final NN-PES. Based on the refined PH2(X2B1) NN-PES, geometries, energies, and harmonic frequencies of stationary points were obtained and analyzed in detail, showing excellent agreement with other theoretical results. Subsequently, quantum time-dependent wave packet (TDWP) and quasi-classical trajectory (QCT) methods were utilized to compute reaction probability, integral cross section (ICS), differential cross section (DCS), product rovibrational distribution, and rate constant on the developed NN-PES. The TDWP results reveal rich resonant structures, while the QCT calculations provide a qualitatively correct description of the coarse-grained reaction cross sections. These results also demonstrate a threshold-type microscopic reaction mechanism characterized by dominant forward–backward scattering, attributable to long-lived collision complex formation. Collectively, these findings provide fundamental mechanistic insights into the microscopic reaction mechanism and dynamics of phosphorus chemistry in interstellar environments.

Tutorial on quantifying and sampling biomolecular ensembles with ShapeGMM

Here, we present a detailed workflow for clustering and enhanced sampling of biomolecular conformations using the ShapeGMM methodology. This approach fits a probabilistic model of biomolecular conformations rooted in the idea that the free energy can be expressed in terms of local fluctuations in atomic positions around metastable states. We demonstrate using a single model system how to generate and fit equilibrium molecular dynamics simulation data. We then demonstrate how to use the resulting model to generate a reaction coordinate between two states, how to sample along that coordinate using metadynamics using our size-and-shape PLUMED module, and how to cluster those biased conformations to obtain a refined equilibrium ShapeGMM model.

Submersion of sodium clusters in helium nanodroplets and critical assessment of He–He pseudopotentials

The submersion of neutral sodium clusters Nan (n = 10, 20, 30) in HeN (N ≲ 90 000) is investigated by zero-point averaged dynamics (ZPAD) simulations, the quantum nature of helium atoms being taken into account through the use of He–He pseudopotentials. Nan clusters are described by tight-binding models; their initial structures are taken from the Quantum Cluster Database, locally optimized, and compared with structures predicted by popular density functionals in order to confirm their suitability. Helium nanodroplets (HNDs) are formed by successive coalescences of smaller droplets, a process that may excite a HND vibrational mode with no influence on the propensity of Nan for submersion in helium. Interestingly, Na20 is found to possibly rearrange at the HND surface at the early times of Na20He79989 dynamics, a behavior absent from other calculations. Three He–He pseudopotentials were mainly tested. Simulations based on the two shallowest pseudopotentials yield swift submersion of Nan (t < 100 ps) at all n investigated here, in disagreement with experiments that predict that submersion should be a slow process occurring when n ≳ 21. The submersion dynamics of Nan is slower with the third and deepest pseudopotential, but calculations up to 250 ps did not enable us to identify a critical size nc ≈ 20 for submersion. We believe that ZPAD can help model the submersion of alkali clusters in HNDs of experimental size in a reasonable amount of time, provided that some effort is made to build He–He pseudopotentials more specific to submersion and large HNDs.