In this paper, we propose a computational framework, based on the VASP and phono3py computer codes, to obtain the thermoelectric figure of merit from the electron-phonon and phonon-phonon interactions using finite displacements in supercells. Several numerical techniques are developed for efficiency. The method is applied to several thermoelectric materials. We found different behaviors for the lifetimes of the electrons in PbTe, PbSe, and in compounds of the half-Heusler and magnesium silicide family. This is traced back to the different frequencies of the phonons involved in the scattering around the Fermi level. They have a lower frequency in PbTe and PbSe. The magnitude of the thermoelectric figures of merit we computed compare well with experiments, but the agreement is far from perfect. The role of the defects, not explicitly considered in our calculations, but abundant in thermoelectric materials, is discussed as a possible explanation. It is also shown that the choice of the exchange-correlation functional can strongly impact the results.
Simulating water from first principles remains a significant computational challenge due to the slow dynamics of the underlying system. Although machine-learned interatomic potentials (MLPs) can accelerate these simulations, they often fail to achieve the required level of accuracy for reliable uncertainty quantification. In this study, we use MACE-an equivariant graph neural network architecture that has been trained using an extensive RPBE-D3 database-to predict density isobars, diffusion constants, radial distribution functions, and melting points. Although equivariant MACE models are computationally more expensive than simpler architectures, such as kernel-based potentials (KbPs), their significantly lower total energy errors allow for reliable thermodynamic reweighting with minimal bias. Our results are consistent with those of previous studies using KbPs; however, equivariant models can be validated against the ground-truth density functional theory (DFT) ensemble with significantly increased efficiency. These findings establish equivariant MLPs as robust and reliable tools for investigating the thermophysical properties of water with DFT-level accuracy.
Abstract Nonadiabatic effects arising from the breakdown of the Born–Oppenheimer approximation govern a wide range of photoinduced processes in molecules and materials. While trajectory surface hopping has become a standard tool for simulating excited-state dynamics in molecular systems, its extension to periodic solids remains comparatively underdeveloped. In this work, we present an interface between the trajectory surface hopping package Surface Hopping including Arbitrary Couplings (SHARC) and the plane-wave density functional theory code Vienna Ab initio Simulation Package (VASP) that enables nonadiabatic dynamics simulations in periodic solids based on first-principles electronic-structure calculations. This development extends the applicability of SHARC from finite molecular systems to extended materials while retaining its modularity and support for arbitrary electronic couplings. Further, it provides a broadly applicable platform for investigating excited-state dynamics in solids using concepts and workflows familiar from the molecular nonadiabatic dynamics community, thereby establishing a unified theoretical and computational framework for both molecular and extended materials. As a first application, the interface is tested on bulk silicon, a prototypical semiconductor whose photoinduced carrier relaxation dynamics is strongly influenced by electron–phonon coupling and has therefore been extensively investigated both experimentally and theoretically.
Predicting tensorial properties with machine learning models typically requires carefully designed tensorial descriptors. In this work, we introduce an alternative strategy for learning tensorial quantities based on scalar descriptors. We apply this approach to the Born effective charge tensor, showing that scalar (monopole) kernel models can successfully capture its tensorial nature by exploiting the definition of the Born effective charge tensor as the derivative of the polarization with respect to atomic displacements. We compare this method with tensorial (dipole) kernel models, as established in our previous work, in which the tensorial structure of the Born effective charge is encoded directly in the kernel and obtained via its derivative. Both approaches are then used for charge partitioning, enabling the separation of monopole and dipole contributions. Finally, we demonstrate the effectiveness of the framework by computing finite-temperature infrared spectra for a range of complex materials.
Variational excited-state density functional theory (DFT) enables the calculation of excited states at a cost comparable to ground-state calculations, but single-configuration approaches often suffer from spin contamination. We implement restricted open-shell Kohn-Sham (ROKS) DFT, which recovers spin-pure singlet excitation energies via the variational minimization of a weighted combination of mixed-spin and triplet configurations, within the plane-wave projector augmented-wave framework of VASP. The energy functional is optimized using a preconditioned conjugate-gradient or a direct inversion in the iterative subspace algorithm, and analytical atomic forces are derived. The implementation is validated for eight organic molecules by comparison to the Q-Chem quantum chemistry code, yielding mean deviations of approximately 30 eV. As a solid-state application, we investigate the three lowest lying excitations of MgO with a neutral oxygen vacancy. For a dielectric-dependent hybrid functional, vertical excitation energies from ROKS and time-dependent density functional theory (TDDFT) differ on average by about 0.21 eV. The Franck-Condon shifts deviate on average by 0.14 eV between the two methods and mass-weighted displacements between the excited states and the ground state by 0.12 amu^1/2 Ang. Additional calculations at the PBE level reveal that these properties depend less strongly on the DFT functional for ROKS than for TDDFT. These results demonstrate that ROKS provides excitation energies and excited-state forces with an accuracy similar to TDDFT while retaining the favorable scaling of ground-state DFT, making it a promising approach for affordable excited-state simulations in extended systems.
Electron-phonon coupling (EPC) is key for understanding many properties of materials such as superconductivity and electric resistivity. Although first principles density-functional-theory (DFT) based EPC calculations are used widely, their efficacy is limited by the accuracy and efficiency of the underlying exchange-correlation functionals. These limitations become exacerbated in complex d- and f-electron materials, where beyond-DFT approaches and empirical corrections, such as the Hubbard U, are commonly invoked. Here, using the examples of CoO and NiO, we show how the efficient r2scan density functional correctly captures strong EPC effects in transition-metal oxides without requiring the introduction of empirical parameters. We also demonstrate the ability of r2scan to accurately model phonon-mediated superconducting properties of the main group compounds (e.g., MgB_2), with improved electronic bands and phonon dispersions over those of traditional density functionals. Our study provides a pathway for extending the scope of accurate first principles modeling of electron-phonon interactions to encompass complex d-electron materials.
Implementing novel features and experimental algorithms into widely adopted density functional theory (DFT) codes is frequently hindered by complex legacy architectures and the use of compiled languages such as Fortran. These production codes, while optimised for high-performance computing clusters, present significant hurdles for software development and rapid prototyping, often requiring deep expertise in the code's internal structure to modify. To address this challenge, we present a Python plugin infrastructure for the Vienna ab-initio Simulation Package (VASP) that combines computational efficiency with the flexibility of high-level scripting. Our architecture uses a C++ intermediate layer and pybind11 to expose VASP data as NumPy arrays via shared memory buffers, ensuring high performance without data duplication. We implement two categories of plugins: those that modify quantities at the end of each converged self-consistent field (SCF) cycle, such as structure and force_and_stress, and those that operate during the SCF cycle, such as local_potential and occupancies. We demonstrate the utility of our implementation through three applications, structure relaxation using the scipy library, implementing an implicit solvent model, and adding the DFT-D4 dispersion corrections. This infrastructure effectively bridges the gap between high-performance electronic structure routines and the widespread scientific Python ecosystem.
We present a plane-wave (PW) implementation of the auxiliary-field quantum Monte Carlo (AFQMC) method within the projector augmented-wave (PAW) formalism in the Vienna ab initio Simulation Package (VASP). By employing an exact inversion of the PAW overlap operator, our approach maintains cubic scaling while naturally incorporating all excitations defined by the PW cutoff. We benchmark this framework by calculating the equilibrium lattice constants and bulk moduli of C, BN, BP, and Si. Our analysis demonstrates that AFQMC systematically corrects the lack of long-range screening in MP2 and the missing higher-order exchange in RPA. We identify RPA as the optimal reference method because of the rapid convergence of the remaining short-range correlations with respect to supercell size. The resulting lattice constants exhibit a mean absolute relative error of 0.14% relative to experiment, establishing the method as a rigorous benchmark tool for structural properties of prototypical semiconductors and insulators.
The atomic structure of the most stable reconstructed surface of cuprous oxide ( Cu 2 O ) ( 111 ) surface has been a longstanding topic of debate. In this study, we develop on-the-fly machine-learned force fields (MLFFs) to systematically investigate the various reconstructions of the Cu 2 O ( 111 ) surface under stoichiometric as well as O- and Cu-deficient or rich conditions, focusing on both ( 3 × 3 ) R 30 ∘ and ( 2 × 2 ) supercells. By utilizing parallel tempering simulations supported by MLFFs, we confirm that the previously described nanopyramidal and Cu-deficient ( 1 × 1 ) structures are the lowest energy structures from moderately to strongly oxidizing conditions. In addition, we identify two promising nanopyramidal reconstructions at highly reducing conditions, a stoichiometric one and a Cu-rich one. Surface energy calculations performed using spin-polarized PBE, PBE + U , r 2 SCAN , and HSE06 functionals show that the previously known Cu-deficient configuration and nanopyramidal configurations are at the convex hull (and, thus, equilibrium structures) for all functionals, whereas the stability of the other structures depends on the functional and is therefore uncertain. Our findings demonstrate that on-the-fly trained MLFFs provide a simple, efficient, and rapid approach to explore the complex surface reconstructions commonly encountered in experimental studies, and also enhance our understanding of the stability of Cu 2 O ( 111 ) surfaces.
Intriguing analogies between the nickelates and the cuprates provide a promising avenue for unraveling the microscopic mechanisms underlying high-Tc superconductivity. While electron correlation effects in the nickelates have been extensively studied, the role of electron-phonon coupling (EPC) remains highly controversial. Here, by taking pristine LaNiO2 as an exemplar nickelate, we present an in-depth study of EPC for both the nonmagnetic (NM) and the C-type antiferromagnetic (C-AFM) phases using advanced density functional theory methods without invoking U or other free parameters. The weak EPC strength lambda in the NM phase is found to be greatly enhanced ('4x) due to the presence of magnetism in the C-AFM phase. This enhancement arises from strong interactions between the flat Ni-3dz2 bands and the low-frequency phonon modes associated with Ni and La vibrations, rather than solely from the high-frequency oxygen breathing modes. The resulting phonon softening is shown to yield a distinctive kink in the electronic structure around 15 meV, which would provide an experimentally testable signature of our predictions. Our study highlights the critical role of local magnetic moments and EPC in the nickelate.
Density functional theory (DFT) is the standard approach for modeling MIL-101(Fe) and related Fe-based metal-organic frameworks, typically assuming a ferromagnetic high-spin configuration. However, this widely adopted approach overlooks a key electronic feature: Spin frustration in the triangular Fe 3 ( μ 3 ${\rm Fe}_{3}(\mu _{3}$ -O) nodes. Using flip-spin, broken-symmetry DFT, we identify the true ground state as an antiferromagnetic 2 S + 1 = 6 $2S+1=6$ state that standard DFT fails to capture. We demonstrate that neglecting spin frustration in MIL-101(Fe) leads to structural distortions, incorrect energetics, and misleading predictions of stability and reactivity. By explicitly accounting for spin frustration, we recover the correct structure and rationalize the temperature-dependent N 2 ${\rm N}_{2}$ and CO binding. Spin frustration enhances N 2 ${\rm N}_{2}$ fixation at room temperature, while its loss upon partial Fe III ${\rm Fe}^{\mathrm {III}}$ reduction suppresses this activity but promotes CO adsorption via π $\pi$ -backbonding. These findings challenge current computational conventions and highlight spin frustration as a critical electronic feature in these frameworks.
Reconstructive phase transitions involving breaking and reconstruction of primary chemical bonds are ubiquitous and important for many technological applications. In contrast to displacive phase transitions, the dynamics of reconstructive phase transitions are usually slow due to the large energy barrier. Nevertheless, the reconstructive phase transformation from β- to λ-Ti3O5 exhibits an ultrafast and reversible behavior. Despite extensive studies, the underlying microscopic mechanism remains unclear. Here, we discover a kinetically favorable in-plane nucleated layer-by-layer transformation mechanism through metadynamics and large-scale molecular dynamics simulations. This is enabled by developing an efficient machine learning potential with near first-principles accuracy through an on-the-fly active learning method and an advanced sampling technique. Our results reveal that the β-λ phase transformation initiates with the formation of two-dimensional nuclei in the ab-plane and then proceeds layer-by-layer through a multistep barrier-lowering kinetic process via intermediate metastable phases. Our work not only provides important insight into the ultrafast and reversible nature of the β-λ transition, but also presents useful strategies and methods for tackling other complex structural phase transitions.
(Received September 2025; accepted 2025; published 2025) We present a constrained random phase approximation (cRPA) method, termed spectral cRPA (s-cRPA), and compare it to established cRPA approaches for scandium and copper by varying the 3d shell filling. The s-cRPA method generally produces larger Hubbard U interaction values compared to conventional approaches. When applied to the realistic system CaFeO3, s-cRPA yields interaction parameters that align more closely with those required within DFT + U to reproduce the experimentally observed insulating state, addressing the metallic behavior predicted by standard density functionals. We examine the issue of negative interaction values encountered in the projector cRPA method for filled d shells. We show that s-cRPA provides improved numerical stability by preserving electron number conservation, a constraint that is violated in the projector cRPA method. The s-cRPA approach addresses some limitations of standard cRPA methods, particularly the tendency to underestimate U values, suggesting its potential utility for the community. Additionally, we have enhanced our implementation to include computation of multicentre interactions for analyzing spatial decay and developed an efficient low-scaling variant employing a compressed Matsubara grid to obtain full frequency-dependent interactions.
Density functional theory (DFT) is the standard approach for modeling MIL-101(Fe) and related Fe-based metal–organic frameworks, typically assuming a ferromagnetic high-spin configuration. This widely adopted approach overlooks a key electronic feature: spin-frustration in the triangular Fe₃O nodes. Using flip-spin DFT, we identify the true ground state as an antiferromagnetic 𝑀 = 6 state that standard DFT fails to capture. We show that using standard DFT for MIL-101(Fe) leads to structural distortions, incorrect energetics, and misleading predictions of stability and reactivity. By explicitly accounting for spin-frustration, we recover the correct structure and rationalize the temperature-dependent N₂ and CO binding: spin-frustration enhances N₂ fixation at room temperature, while its loss upon partial Fe(III) reduction suppresses this activity but promotes CO adsorption via 𝜋-back bonding. These findings challenge the current conventions and highlight spin-frustration as a critical electronic feature.
Density functional theory (DFT) calculations of charged molecules and surfaces are critical to applications in electro-catalysis, energy materials and related fields of materials science. DFT implementations such as the Vienna ab-initio Simulation Package (VASP) compute the electrostatic potential under 3D periodic boundary conditions, necessitating charge neutrality. In this work, we implement 0D and 2D periodic boundary conditions to facilitate DFT calculations of charged molecules and surfaces respectively. We implement these boundary conditions using the Coulomb kernel truncation method. Our implementation computes the potential under 0D and 2D boundary conditions by selectively subtracting unwanted long-range interactions in the potential computed under 3D boundary conditions. By combining the Coulomb kernel truncation method with a computationally efficient padding approach, we remove nonphysical potentials from vacuum in 0D and 2D systems. To illustrate the computational efficiency of our method, we perform large supercell calculations of the formation energy of a charged chlorine defect on a sodium chloride (001) surface and perform long time-scale molecular dynamics simulations on a stepped gold (211) | water electrode-electrolyte interface.
Electron-phonon coupling (EPC) is fundamental for understanding the behavior of molecules and crystals, influencing phenomena such as charge transport, energy transfer, phase transitions, and polaron formation. Accurate computational methods to calculate EPCs from first principles are essential, but their complexity has resulted in a variety of computational strategies, raising concerns about their mutual consistency. In this study, we provide a systematic benchmark of methods for EPC calculation by comparing two fundamentally different ab initio methodologies. We investigate Gaussian-type orbital methods based on the CP2K code and plane-wave-based projector-augmented-wave methods combined with maximally localized Wannier functions, as implemented in VASP and wannier90. In addition, we further distinguish between the derivative-of-Hamiltonian ( dH) and derivative-of-states ( d psi) approaches for obtaining EPC parameters. The comparison is conducted on a representative set of organic molecules, including pyrazine, pyridine, bithiophene, and quarterthiophene, varying significantly in size and flexibility. We find excellent agreement across implementations and basis sets when employing the same computational approach ( dH or d psi), demonstrating robust consistency between the numerical schemes. However, noticeable deviations occur when comparing the dH and d psi approaches within each code and for specific cases discussed in detail. Our findings emphasize the reliability of EPC computations using the dH method and caution against potential pitfalls associated with the d psi approach, providing guidance for future EPC calculations and model parameterizations.
While the periodic equation-of-motion coupled-cluster (EOM-CC) method promises systematic improvement of electronic band gap calculations in solids, its practical application at the singles and doubles level (EOMCCSD) is hindered by severe finite-size errors in feasible simulation cells. We present a hybrid approach combining EOM-CCSD with the computationally less demanding GW approximation to estimate thermodynamic limit band gaps for several insulators and semiconductors. Our method substantially reduces required cell sizes while maintaining accuracy. Comparisons with experimental gaps and self-consistent GW calculations reveal that deviations in EOM-CCSD predictions correlate with reduced single excitation character of the excited many-electron states. Our work not only provides a computationally tractable approach to EOM-CC calculations in solids but also reveals fundamental insights into the role of single excitations in electronic-structure theory.
Correction for ‘Absolute standard hydrogen electrode potential and redox potentials of atoms and molecules: machine learning aided first principles calculations’ by Ryosuke Jinnouchi et al., Chem. Sci., 2025, 16, 2335–2343, https://doi.org/10.1039/D4SC03378G.
We revisit the long-standing question of whether water molecules dissociate on the Ru(0001) surface through nanosecond-scale path-integral molecular dynamics simulations on a sizable supercell. This is made possible through the development of an efficient and reliable machine-learning potential with near first-principles accuracy, overcoming the limitations of previous ab initio studies. We show that the quantum delocalization associated with nuclear quantum effects enables rapid and frequent proton transfers between water molecules, thereby facilitating the water dissociation on Ru(0001). This work provides the direct theoretical evidence of water dissociation on Ru(0001), resolving the enduring issue in surface sciences and offering crucial atomistic insights into water-metal interfaces.
We investigate the relationship between the K-edge fine structure of isolated single-wall carbon nanotubes (SWCNTs) and the Van Hove singularities (VHSs) in the conduction-band density of states. To this end, we model x-ray absorption spectra of SWCNTs using the final-state approximation and the Bethe-Salpeter equation (BSE) method. Both methods can reproduce the experimental fine structure, where the BSE results improve on peak positions and amplitude rations compared to the final-state approximation. When the fine structure in the modeled spectra is related to the VHSs, significant differences are found. We suggest that these differences arise due to modifications of the core exciton wave functions induced by the confinement along the circumference. Additionally, we analyze the character of core excitons in SWCNTs, and we find that the first bright excitons are Frenkel excitons, while higher-lying excitons are charge resonance states. Finally, we suggest that the qualitative picture based on VHSs in the density of states holds when there is a large energy gap between successive VHSs.