External electric fields play a crucial role in chemistry, materials, and biology. However, accurately capturing external field effects remains a challenge for machine learning interatomic potentials (MLIPs). By revealing the limitations of conventional architectures in predicting field-dependent response properties, we introduce a physics-informed field-aware equivariant neural network framework that predicts response properties to external electric fields. External electric fields are incorporated through an equivariant embedding layer that concatenates field representations with atomic features, preserving rotational equivariance. The model leverages an energy-derivative approach to predict response properties in a single unified architecture. Furthermore, we embed a charge-equilibration electrostatic layer into the model to capture long-range electrostatic interactions. The close agreement between charged-based and energy derivative-based electric dipoles indicates that the model learns the correct underlying physics. The unified framework performs robustly for the prediction of response properties and simulation of IR and Raman spectra across molecular and periodic systems under different external field strengths. By incorporating physical principles into the model, we establish a general and transferable framework for modeling molecular and condensed-phase response under external electric fields, providing a robust routine for response properties prediction in external electric fields.

Singlet excited electronic states can be directly constructed with the ΔSCF method using a single electronic density by applying fractional occupation numbers. While the advantages and disadvantages of such a DFT-based ΔSCF method were demonstrated in several studies, the spurious energy shift of the singlet electronic states when constructed using hybrid DFT functionals remains unexplained. Here, we explain in detail the origin of the effect of losing the idempotence properties of density matrices on the DFT energy terms with hybrid-based DFT functionals, as well as a procedure to eliminate the artificial shift from ΔSCF singlet excited state energies when obtained using the hybrid-based DFT functionals.

We present a deterministic all-parameter optimization method based on analytic derivatives for vibrational wave functions in a separable product form represented in a distributed Gaussian basis. The optimization of the center and width parameters of the basis allows for compact wave function representations at a small basis size. We illustrate the approach by computing ground and excited vibrational states of anharmonic model systems and the 9-dimensional formic acid monomer.

Machine learning interatomic potentials (MLIPs) offer a computationally efficient and accurate approach for simulating complex chemical systems at the atomistic level. However, their predictive capabilities require further refinement to address broader chemical phenomena. This study introduces a charge-and spin-aware ability to equivariant graph neural network MLIPs that has the ability to predict energies, forces, and atomic charges and atomic spin moments in chemical systems. A novel spin-dependent charge equilibration (QEq) method is proposed, extending the MLIP’s applicability to spin-polarized systems. The spin-polarized QEq-enabled MACE models were benchmarked against density functional theory reference datasets, achieving superior predictive accuracy. Notably, the model accurately reproduces polaron distributions on the TiO2 (110) surface, demonstrating its robustness in capturing spin-related properties. This advancement significantly improves the versatility and precision of MLIPs, enabling more reliable atomistic simulations in chemistry and materials science.

Here we present a density matrix based KS inversion method formulated entirely within a Gaussian basis representation to optimize a KS potential matrix that reproduces a target electron density. Inverse Kohn-Sham (KS) density functional theory (DFT) aims to determine the effective local KS potential that reproduces a target electron density, and is important both for electronic structure analysis and for the development of orbital based correction methods. In finite Gaussian basis implementations, however, conventional inverse KS-DFT approaches such as the Zhao-Morrison-Parr (ZMP) method often become poorly constrained and inefficient, because the real space penalty potential is projected onto a limited number of Gaussian basis matrix elements, which can strongly coarse-grain its spatial variation. In the present method, the density matrix mismatch is defined in a Lowdin orthogonalized basis, which yields a penalty energy invariant under unitary rotations in that basis. The corresponding penalty potential contribution to the KS Hamiltonian is derived analytically in the original nonorthogonal Gaussian basis. Across a wide range of penalty strengths, the self consistent field (SCF) optimization remains robust and efficient for various open shell systems, while progressively tightening the penalty drives the electron density into accurate agreement with the target. Benchmarks on molecules and condensed phase systems show that the method achieves substantially smaller attainable density deviations than the conventional ZMP method. The method provides a fast and accurate route to KS inversion in finite Gaussian basis sets and may also be useful for future orbital based correction schemes.

Integration of Ru(bda) [bda = (2,2'-bipyridine)-6,6'-dicarboxylate] water oxidation catalysts together with dibenzofuran, dibenzothiophene, or carbazole heterocycles within macrocycles 1-3 affords structures with different conformational preferences. The latter are dictated by noncovalent interactions of the bda oxygens with NH (hydrogen bond) or S (chalcogen bond). The different orientations afford changes in the secondary coordination sphere, thereby influencing the rate and mechanism of photocatalytic water oxidation for the three mononuclear Ru(bda)-based catalysts, with the dibenzothiophene- and carbazole-based ones favoring the unimolecular pathway (first-order kinetics), typically related to water nucleophilic attack (WNA), and the dibenzofuran one operating via a bimolecular pathway (second-order kinetics), typically related to the interaction of two metal-oxo species (I2M). NMR and single crystal X-ray analyses provided insight into the conformational preferences of the precursor macrocycles in the Ru(II) state. Photocatalytic studies including the H/D kinetic isotope effect (KIE) in deuterated water together with theoretical studies on the orientation angle-dependent energy profile for macrocycles in the Ru(V) state afforded a structure-property relationship that explains the outcome of the water oxidation experiments.

We implemented ab initio Hubbard parameter calculation schemes in the k-point sampling real-time TDDFT (RT-TDDFT) program in CP2K. We propose a new linear-response-based calculation scheme for energy-dependent Hubbard parameters. Our scheme extends the minimum-tracking linear-response method proposed in [Moynihan et al., arXiv preprint arXiv:1704.08076(2017); E. B. Linscott et al., Phys. Rev. B 98, 235157 (2018)] to realize the calculation of energy-dependent Hubbard parameters that reflect the exchange-correlation (xc) effects included in the xc-functional. We discuss the properties of the minimum-tracking linear-response method in comparison to another promising scheme, ACBN0 [Agapito et al., Phys. Rev. X, 5, 011006 (2015)]. We show that, while neither clearly outperforms the other in the accuracy of static property calculations, each has a distinct dynamical application depending on its theoretical formulation.

Perovskite oxynitride LaTiO2N holds promise for visible-light-driven photocatalytic water splitting, yet its surface dynamics at the catalyst-water interface remain elusive. This study employs state-of-the-art density functional theory molecular dynamics (DFT-MD) and electrochemical measurements to unravel the intricate interplay of water arrangement and adsorption on the LaTiO2N(100) surface. By considering explicit solvent effects, we reveal a pronounced hydrophilic character, with spontaneous water dissociation forming hydroxyl groups at undercoordinated Ti sites at the surface, stabilizing the latter and potentially enhancing the OER activity. Our simulations identify the thermodynamically stable termination and demonstrate its role in fostering robust hydrogen-bond networks that facilitate proton-coupled electron transfer. The surface Pourbaix diagram underscores hydroxylated configurations across diverse electrochemical conditions, corroborated by DFT-MD insights. These findings highlight the critical role of treating an explicit solvent environment and its dynamics at a given temperature in the modeling of LaTiO2N, offering a blueprint for designing high-performance, sustainable water-splitting catalysts.

Harnessing functional groups in the outer coordination sphere to direct catalytic pathways is characteristic for enzymes, yet rather underexplored for man-made catalysts. Here, following our earlier work on Ru(bda) water oxidation catalysts (bda = 2,2'-bipyridine-6,6'-dicarboxylate), we demonstrate how a carboxyl-functionalized Ru(bda)-based macrocycle enables an oxide relay mechanism in light-driven water oxidation via an O─O bond formation between carboxylate and RuV = O units positioned on opposing sides. Through a combination of photocatalytic water oxidation studies, NMR and single crystal analysis, and 18O-labeling experiments, we provide evidence for the mechanistically distinct oxide relay pathway under photocatalytic conditions. Our findings underscore the role of second coordination sphere engineering in modulating reaction pathways and advancing molecular catalyst design.

Here we present a density matrix penalization method for finite basis Kohn-Sham (KS) matrix reconstruction in Gaussian basis representations. The method constructs a matrix-represented auxiliary KS Hamiltonian whose self-consistent density matrix and corresponding Gaussian-basis-representable real-space electron density reproduce a prescribed target, without assuming the recovery of a unique continuum local KS or exchange-correlation potential. This finite basis formulation is motivated by the numerical difficulties encountered by conventional inverse KS-DFT approaches, such as the Zhao-Morrison-Parr (ZMP) method, when a real-space penalty potential is projected onto a limited set of Gaussian basis matrix elements. Such a projection can strongly coarse-grain the constraining potential and lead to poorly constrained and inefficient self-consistent-field (SCF) optimizations. In the present method, the density matrix mismatch is defined in a Löwdin-orthogonalized basis, yielding a penalty functional that is invariant under basis rotations in that representation. The corresponding penalty Hamiltonian matrix contribution is derived analytically in the original nonorthogonal Gaussian basis. Across a wide range of penalty strengths, SCF optimization remains robust and efficient for various open-shell molecular and condensed-phase systems, while progressively tightening the penalty drives the density matrix and the associated real-space density into near machine precision agreement with the target for most systems. Benchmarks show that the method achieves substantially smaller attainable density deviations than conventional ZMP calculations in Gaussian bases. The method provides a stable, fast, and accurate route to finite basis KS matrix reconstruction and establishes a practical framework for density-matrix-based inverse reconstruction in Gaussian basis electronic structure calculations.

Accurate prediction of molecular electric dipole moments is crucial for determining molecular properties and understanding molecular reactivity. Recent developments in equivariant machine learning enable direct prediction of electric dipole moments as vector quantities. This raises the question as to whether it is still necessary to consider physics-informed charge-based principles in the model to achieve satisfactory performance. In this work, we systematically compare direct equivariant electric dipole prediction using MACE with charge-based approaches using adapted MACE which incorporates a variant charge equilibrium equation (QEq) into the MACE framework. Across diverse datasets including QM7b, QM9, as well as subsets of the SPICE and SN2 data sets, we demonstrate that both direct equivariant electric dipole prediction and physics-informed QEq approaches show good performance for short-to-medium range interaction datasets. The charge-based QEq model shows slightly better performance than direct prediction beyond small training data sizes. Notably, the charge-based QEq model outperforms direct electric dipole prediction for the systems with long range interactions. Our results highlight the importance of incorporating physics into the model for improved model interpretability and transferability.

We developed a k-point sampling real-time TDDFT (RT-TDDFT) program within the Gaussian and plane waves (GPW) framework of the CP2K software suite. In addition to standard real-time propagation of time-dependent Kohn-Sham orbitals, we make use of symmetry-based k-point reduction and k-point parallelization schemes so that our RT-TDDFT program in the GPW framework is feasible for practical large-scale calculations. We also implemented DFT + U as a relevant extension for real-time simulations of systems with strong electron correlations. In particular, we extended the "tensorial" subspace representation approach for DFT + U, following the formulation in [Chai, Z., et al. J. Chem. Theory Comput., 2024, 20, 8984], to k-point sampling RT-TDDFT. Our extension, which is, to our knowledge, the first application of the "tensorial" subspace representation approach to k-point sampling RT-TDDFT, is found to be robust and efficient with small additional costs owing to the locality of Gaussian basis functions, indicating that it is a promising approach to RT-TDDFT + U for solids. We show details of our implementation in CP2K and the results of our benchmark calculations.

We propose an anisotropic interfacial continuum solvation (AICS) model to simulate the distinct in-plane and out-of-plane dielectric constants of liquids near solid-liquid interfaces and their spatial variations along the surface normal direction. In low-electron-density regions, each dielectric function in the diagonal components of a dielectric tensor varies monotonically with distance from the solid surface along the surface normal direction; in high-electron-density regions near the surface, each dielectric function adopts the electron-density-based formulation proposed by Andreussi et al. (J. Chem. Phys. 2012, 136, 064102). The resulting dielectric tensor is continuously differentiable with respect to both electron density and spatial coordinates. We derived analytical expressions for electrostatic contributions to the Kohn-Sham potential and atomic forces and implemented the AICS model, including these analytical derivatives, into the CP2K software package. To solve the anisotropic Poisson equations, we developed a parallel finite-element anisotropic Poisson solver (FEAPS) based on the FEniCSx platform and its interface with CP2K. Analytical forces were validated against finite-difference calculations, while electrostatic potentials computed under vacuum and isotropic solvent conditions using AICS and FEAPS were benchmarked against standard vacuum DFT and SCCS results, respectively. In the anisotropic interfacial solvent environment characterized by the enhanced in-plane and reduced out-of-plane dielectric functions near the Ag(111) surface, we calculated the resulting work functions and electrostatic potentials and optimized the adsorption geometry for OH*. Compared to the isotropic case, we observed more pronounced work function shifts and spatially modulated electrostatic profiles across different charge states. Our results also showed that OH* tilted more toward the plane parallel to the surface under the anisotropic dielectric conditions.

With sunlight as the most abundant energy source on earth, solar water splitting has the potential to produce renewable hydrogen at a commercially competitive cost. Monoclinic BiVO4 is a promising n‐type semiconductor for photocatalytic and photoelectrochemical (PEC) water oxidation. One of the simplest and most energy‐efficient approaches for producing BiVO4 is hydrothermal synthesis. This method is carried out at moderate temperatures, while particle size, shape, and crystallinity are controlled by a wide range of synthesis parameters, which can be further expanded by using additives. In this work, these parameters systematically vary to study their influence on the hydrothermal synthesis of BiVO4, with a focus on KCl as an additive are systematically vary. By X‐ray diffraction, scanning electron microscopy, and transmission electron microscopy is shown that KCl acts as structure‐directing agent, leading to significant changes in morphology and crystallinity. Since the color and the optical spectra of BiVO4 powders indicate a redshift with increasing KCl concentration, an additional anionic substitution by Cl− takes place is proposed, a hypothesis supported by X‐ray photoelectron spectroscopy measurements and density functional theory calculations. The highest photocatalytic performance (1328 µmol g−1 h−1) is reached for 25 mmol L−1 KCl, while particle‐based photoelectrodes decorated with CoPi cocatalysts showed an improved photocurrent density (393 µA cm−2) at 1.23 V vs. reversible hydrogen electrode (RHE).

A thorough investigation of the delta self-consistent field (ΔSCF) method within the restricted-open Kohn-Sham (ROKS) formalism has been carried out. The ROKS-based ΔSCF targets unmixed singlet excited electronic states, avoiding the need for a spin purification procedure. A modification of the maximum overlap method to improve the convergence of the ΔSCF method is presented, in addition to a molecular orbital tracking and order-preserving algorithm. For benchmarking purposes a large-scale comparison with time-dependent density functional theory (TDDFT) was conducted: both single- and multi-reference singlet and triplet excited electronic states of various molecules were reproduced using ΔSCF and compared to the states provided by TDDFT. Besides the excitation energies, we also compared the electron densities and the transition dipole moments of the molecules between the two methods.

The CP2K software package provides a comprehensive suite of density functional theory-based methods for studying excited states and spectroscopic properties of molecular and periodic systems. In this review, we present recent developments and applications of several complementary approaches implemented in CP2K, including linear-response time-dependent (TD) and time-independent density functional perturbation theory (DFPT), delta self-consistent field (ΔSCF), and real-time TDDFT (RT-TDDFT). Nonadiabatic molecular dynamics (NAMD) capabilities are integrated with ΔSCF and TD-DFPT methods, in addition to Ehrenfest dynamics based on RT-TDDFT, enabling detailed investigations of photochemical processes and the excited-state dynamics in gas and condensed phase systems. Applications demonstrating the versatility of these methods include studies on solvated molecules, surface-bound photosensitizers, and two-dimensional materials. Spectroscopic methods encompass, e.g., ultraviolet-visible absorption, electronic circular dichroism, Raman (optical activity), infrared absorption, and vibrational circular dichroism spectra. We demonstrate that CP2K provides a unique and powerful toolkit for studying a wide range of excited-state phenomena in complex molecular and extended (periodic) systems.

Electrochemical water splitting is essential for reducing our dependence on fossil fuels through green hydrogen production and requires the design of new, low-cost, 3d transition-metal-based catalysts for the sluggish oxygen evolution reaction (OER). We report on the synthesis, characterization, and OER performance of new types of homo- and heterometallic iron(II)-based di(2-pyridyl)ketone cubanes, 1-[Fe4(dpy-C{OH}O)4(OAc)3(H2O)]ClO4 and 2-[Fe2Ni2(dpy-C{OH}O)4(OAc)3]ClO4 (dpy = di(2-pyridyl)diol), referred to as 1-{Fe4O4} and 2-{Fe2Ni2O4}. The heterometallic oxocluster is the first sought-after molecular cutout of the key active {H2O-Fe2Ni2(OR)2-OH2} motif, bridging molecular and heterogeneous OER catalysts as a model to understand the key catalytic synergisms in powerful NiFe oxide-based materials. The precise positions of Fe(II) and Ni(II) cations within the 2-{Fe2Ni2O4} oxocluster, along with the detailed structural and electronic features, were elucidated with a variety of methods, including extended advanced X-ray absorption spectroscopy and total neutron scattering techniques, accompanied by time-dependent density functional theory (TDDFT) calculations. (Spectro)electrochemical analyses showed that 2-{Fe2Ni2O4} exhibits synergisms between Fe(II) and Ni(II) centers, tuning both metals' redox potentials and the molecular stability of the mixed-valent {(FeIII)2Ni2O4} precatalyst species. 2-{Fe2Ni2O4} displays a maximum jcat of 1.75 mA/cm2 for the OER and a Faradaic efficiency of 84.8% under turnover conditions. Although the oxocluster experienced notable degradation during the OER, significant metal (oxy)hydroxide formation could be excluded. 2-{Fe2Ni2O4} is introduced as a new molecular platform for in-depth studies of NiFe synergisms, along with ligand engineering or polymer matrix strategies.

We present electric dipole polarizability calculations employing atomic-orbitals based linear response theory within the Kohn-Sham Density Functional Theory (KS-DFT) framework, considering both non-periodic and periodic boundary conditions. We adopt the optimization scheme introduced by T. Helgaker et al. in Chemical Physics Letters 327, 397 (2000) for the single-electron atomic-orbitals density matrix. We conduct a comparative analysis between the static polarizability computed using atomic orbitals-based and previously implemented molecular orbitals-based methods. In our calculations involving periodic boundary conditions, we implement polarizability calculation using velocity representation of the electric dipole operator in atomic orbitals-based algorithm, subsequently comparing the results with those computed using the Berry-phase formulation and velocity representation in molecular orbitals-based algorithm. We investigate 10 small and medium-sized molecules in the gas phase, analyze liquid-phase systems with up to 256 water molecules, and the solid-state structures of anatase TiO2 and bulk WO3. All polarizability results obtained from the AO-based solver exhibit good agreement with MO-based results. From our example calculations, we find that the AO-based solver exhibits better computational scaling and less memory demand than the MO-based solvers, which makes it better suited for very large systems.

In the self-consistent continuum solvation (SCCS) approach (J. Chem. Phys. 136, 064102 (2012)), the analytical expressions of the local solute-solvent interface functions determine the interface function and dielectric function values at a given real space position based solely on the electron density at that position, completely disregarding the surrounding electron density distribution. Therefore, the low electron density areas inside the solute will be identified by the algorithm as regions where implicit solvent exists, resulting in the emergence of non-physical implicit solvent regions within the solute and even potentially leading to the divergence catastrophe of Kohn- Sham SCF calculations. We present a new and efficient SCCS implementation based on the solvent-aware interface (J. Chem. Theory Comput. 15, 3, 1996-2009 (2019)) which addresses this issue by utilizing a solute-solvent interface function based on convolution of electron density in the CP2K software package, which is based on the mixed Gaussian and plane waves (GPW) approach. Starting with the foundational formulas of SCCS, we have rigorously derived the contributions of the newly defined electrostatic energy to the Kohn-Sham potential and the analytical forces. This comprehensive derivation, which to the best of our knowledge is not available in the current literature, utilizes the updated versions of the solute-solvent interface function and the dielectric function, tailored to align with the specifics of the GPW implementation. Our implementation has been tested to successfully eliminate non-physical implicit solvent regions within the solute and achieve good SCF convergence, as demonstrated by test results for both bulk and surface models, namely liquid H2O, titanium dioxide, and platinum.