Stimulated by the renewed interest and recent developments in semi-empirical quantum chemical (SQC) methods for noncovalent interactions, we examine the properties of liquid water at ambient conditions by means of molecular dynamics (MD) simulations, both with the conventional NDDO-type (neglect of diatomic differential overlap) methods, e.g. AM1 and PM6, and with DFTB-type (density-functional tight-binding) methods, e.g. DFTB2 and GFN-xTB. Besides the original parameter sets, some specifically reparametrized SQC methods (denoted as AM1-W, PM6-fm, and DFTB2-iBi) targeting various smaller water systems ranging from molecular clusters to bulk are considered as well. The quality of these different SQC methods for describing liquid water properties at ambient conditions are assessed by comparison to well-established experimental data and also to BLYP-D3 density functional theory-based ab initio MD simulations. Our analyses reveal that static and dynamics properties of bulk water are poorly described by all considered SQC methods with the original parameters, regardless of the underlying theoretical models, with most of the methods suffering from too weak hydrogen bonds and hence predicting a far too fluid water with highly distorted hydrogen bond kinetics. On the other hand, the reparametrized force-matchcd PM6-fm method is shown to be able to quantitatively reproduce the static and dynamic features of liquid water, and thus can be used as a computationally efficient alternative to electronic structure-based MD simulations for liquid water that requires extended length and time scales. DFTB2-iBi predicts a slightly overstructured water with reduced fluidity, whereas AM1-W gives an amorphous ice-like structure for water at ambient conditions.

We present an implementation for the calculation of K-edge resonant inelastic X-ray scattering spectra based on time-dependent density functional theory in the CP2K package. The method evaluates the Kramers-Heisenberg cross section from transition dipole moments connecting the ground state, core-excited intermediate states, and valence-excited final states. The present implementation combines the existing linear response time-dependent density functional theory modules (XAS-TDP and TDDFT) within a unified framework. The resulting approach is computationally efficient and well suited for condensed-phase applications in combination with ab initio molecular dynamics. The implementation is validated against experiment and established reference calculations for gas-phase water and methanol at the oxygen K-edge. Its applicability to realistic condensed-phase systems is then demonstrated for crystalline kaolinite at the oxygen K-edge, and for aqueous ammonia at the nitrogen K-edge, where representative configurations are sampled from ab initio molecular dynamics trajectories.

CP2K is a versatile open-source software package for simulations across a wide range of atomistic systems, from isolated molecules in the gas phase to low-dimensional functional materials and interfaces, as well as highly symmetric crystalline solids, disordered amorphous glasses, and weakly interacting soft-matter systems in the liquid state and in solution. This review highlights CP2K's capabilities for computing both static and dynamical properties using quantum-mechanical and classical simulation methods. In contrast to the accompanying theory and code paper [J. Chem. Phys. 152, 194103 (2020)], the focus here is on the practical usage and applications of CP2K, with underlying theoretical concepts introduced only as needed.

Ab initio quantum Monte Carlo (QMC) methods are state-of-the-art electronic structure calculations based on highly parallelizable stochastic frameworks for accurate solutions of the many-body Schrödinger equation, suitable for modern many-core supercomputer architectures. Despite its potential, one of the major drawbacks that still hinders QMC applications, especially when targeting dynamical properties of large systems or extensive datasets, is the lack of an affordable method to compute atomic forces that are consistent with the corresponding potential energy surfaces (PESs), also known as unbiased atomic forces. Recently, one of the authors in the present paper proposed a way to obtain unbiased forces with the Jastrow-correlated Slater determinant Ansatz, where the determinant part is frozen to the values obtained by a mean-field method, such as density functional theory [K. Nakano, M. Casula, and G. Tenti, Phys. Rev. B 109, 205151 (2024)]. However, the proposed method has a significant drawback for its applications: for a system with N nuclei, one requires 6N additional density functional theory (DFT) calculations to get unbiased forces, which is not negligible as the system size increases. This paper presents a way to replace the 6N DFT calculations with a single coupled-perturbed Kohn-Sham calculation, following the so-called Lagrangian technique established in quantum chemistry. This improves the computational cost and scalability of the method. We also demonstrate that the developed unbiased variational Monte Carlo (VMC) force calculation improves not only the consistency with PESs but also its accuracy, by investigating three molecules from the rMD17 benchmark set, and comparing the unbiased VMC forces with those obtained by the coupled-cluster singles and doubles with perturbative triples [CCSD(T)] calculations. We found that the bare VMC forces are biased from the CCSD(T) ones, while the unbiased ones give values closer to those of the CCSD(T) ones. Our benchmark test also reveals that the unbiased VMC forces yield very consistent values with hybrid and meta generalized gradient approximations (e.g., ωB97X-D3BJ and ωB97M-D3BJ), but do not necessarily yield values that are very close to those of CCSD(T). Our finding paves the way to generate machine learning interatomic potentials based on VMC forces more efficiently and accurately.

One of the distinguishing aspects of CP2K is its seamless integration of diverse structural and transition-state optimization techniques with advanced sampling approaches including Monte Carlo, molecular dynamics, and metadynamics, enabling the efficient exploration of complex potential- and free-energy landscapes, including rare events. These capabilities are combined with a broad hierarchy of energy and force evaluation methods, ranging from classical and machine-learned interaction potentials and mixed quantum-classical multiscale and semiempirical schemes, to highly accurate quantum-mechanical electronic-structure approaches. At the heart of the latter lies the Gaussian and plane-wave framework, along with its augmented all-electron generalization, which have been described in detail in our previous code review [T. D. Kühne et al., J. Chem. Phys. 152, 194103 (2020)]. Building on this foundation, the present work revisits the methods within CP2K that turn electronic structure into dynamics, transport, and spectroscopic response. Particular emphasis is placed on the coupling between static response calculations and nuclear motion: spectra may be evaluated at optimized structures, averaged over thermally sampled configurations, obtained from time-correlation functions along ab-initio or path integral molecular trajectories, or followed in real time together with electronic and nuclear dynamics. The same modular structure also enables equilibrium and biased transport simulations, from Kubo-type linear response to open-boundary approaches under external potentials, highlighting CP2K's unique capability to unify quantum chemistry with quantum and statistical mechanics within a versatile, holistic simulation environment.

Machine-learned exchange–correlation (XC) functionals offer a route to improve Kohn–Sham density-functional theory without incurring the cost of explicitly correlated electronic-structure methods. Their use in production simulation codes, however, requires a well-defined mapping between the learned model and the host-code density representation. We formulate and implement a Skala-1.1 interface in CP2K through the external GauXC library. CP2K supplies the geometry, Gaussian basis, spin-resolved atomic-orbital density matrix, and communicator, while GauXC evaluates the XC energy, atomic-orbital potential matrix, and available nuclear derivatives. The interface accepts both all-electron and valence-only density matrices. The latter may arise from separable dual-space pseudopotentials or molecular effective-core potentials. Implementation errors are isolated from functional differences by comparing the Perdew–Burke–Ernzerhof (PBE) functional evaluated through GauXC with native CP2K PBE. The resulting interface gives consistent energies, forces validated against finite-difference total-energy checks, and force-based molecular-virial diagnostics for representative molecular cases. The dietGMTKN55 benchmark suite is evaluated with an all-electron Gaussian augmented plane-wave treatment for elements up to bromine and def2 effective-core potentials for the heavier elements. The resulting aggregate mean absolute deviation of 1.255 kcal/mol is within 0.020 kcal/mol of the corresponding Skala reference value of 1.235 kcal/mol. This work establishes a validated molecular implementation of Skala in CP2K through GauXC.

Kohn-Sham density functional theory is the workhorse of computational materials science and chemistry, but the optimal choice of exchange-correlation approximation for a given system is often unclear. Traditional benchmark studies typically compare selected density functional approximations (DFAs) against high-level quantum chemistry results for small structures and static energies, rather than against experimental thermodynamic or response properties. Here, we show how machine learning interatomic potentials (MLIPs) allow for molecular dynamics simulations of these properties at scale, thereby providing a thermodynamic benchmark for DFAs. Focusing on the wellstudied but challenging case of bulk water at ambient conditions, we benchmark 11 DFAs spanning the top four rungs of Jacob's ladder (GGA, meta-GGA, hybrid, and double hybrid) with and without empirical dispersion corrections. Using MLIPs trained on each DFA, we compute infrared spectra, heat capacities, radial distribution functions, and diffusivities. We find a general improvement up Jacob's ladder towards double hybrids, though there is substantial variation between DFAs within the same rung and across different properties. Dispersion corrections yield only a modest influence on the computed properties. We elucidate the similarities and hierarchy among DFAs, both based on correlating the DFA energies and forces for the same training configurations, and based on analyzing the finite-temperature observables. Our work thus demonstrates the power of MLIPs for finite-temperature benchmarking of electronic structure methods.

Metadynamics enables the reconstruction of free-energy surfaces from molecular dynamics; however, its application is often challenging because numerous methodological choices can influence the results. Here, we focus on how geometric confinement reshapes reactant and product basins and thereby contributes to catalytic effects. Our analysis is based on two prototypical reactions: a symmetric SN2 fluorine-substitution reaction, CH3F + F- ⇌ F- + CH3F and a Diels-Alder cycloaddition. We compare several metadynamics approaches: well-tempered metadynamics, on-the-fly probability-enhanced sampling (OPES), OPES-Explore, and a hybrid OPES + OPES-Explore scheme, and we examine two types of confinements: First, we analyze the influence of methodological confinement introduced by restraining walls, a commonly used but seldom assessed component of enhanced sampling simulations. Second, we examine chemical confinement inside carbon nanotubes of varying diameters. We find that the restraining wall does not alter the activation free energy of the studied reactions but significantly affects the widths of the reactant and product basins. In contrast, confinement inside carbon nanotubes changes the barrier height by enforcing axial alignment of the reactants in the SN2 reaction and by restricting the transition-state and product geometries in the Diels-Alder case. Concerning the metadynamics methods, we found different convergence behavior and sampling quality, even for these simple reactions. The results highlight how both physical and methodological confinement can influence the outcome of enhanced-sampling simulations, and they underscore the need for careful choice of metadynamics methods for reactions in restricted environments.

Reliable density-functional simulations require numerical settings whose residual errors are smaller than the chemical and materials trends being interpreted. In CP2K/QUICKSTEP, this requirement is complicated by the joint use of atom-centered Gaussian basis sets and norm-conserving pseudopotentials: a code-to-code discrepancy usually contains both contributions. We present the UZH protocol, a closed-loop CP2K workflow that calibrates molecularly optimized Gaussian basis sets on small molecules, validates the resulting settings in unary-crystal equation-of-state benchmarks, and identifies whether the limiting approximation is the Gaussian basis or the pseudopotential. The diagnosis is then used to revise the parameter files. The central diagnostic is a three-way comparison between production CP2K-GTH-UZH calculations, SIRIUS calculations using the same Goedecker-Teter-Hutter pseudopotential in a systematic plane-wave representation, and all-electron full-potential linearized augmented-plane-wave SIRIUS references. This construction decomposes the practical CP2K error into a Gaussian-basis component and a pseudopotential component. The protocol distinguishes basis-limited noble-gas and heavy-element cases from pseudopotential-limited transition-metal cases, guides targeted revisions with the CP2K basis and pseudopotential optimizers, and produces improved MOLOPT basis sets and GTH pseudopotentials as explicit outputs of the workflow. The UZH protocol is therefore constructive: it does not merely measure or reduce errors a posteriori but allows turning verification outliers into validated CP2K parameter files for simulations across molecules and condensed phases. The resulting MOLOPT/GTH parameter files are proposed as a practical default for future CP2K calculations.

We report the static and dynamical properties of liquid water at second-order Møller-Plesset perturbation theory level (MP2) with classical and quantum dynamics simulations using a neural network potential. We examined the temperature-dependent radial distribution function, diffusion and vibrational dynamics. MP2 theory predicts an over-structured liquid water at ambient conditions, which may be attributed to the incomplete basis set. The excellent agreement with experimental structural properties as well as the diffusion constant is observed at an elevated temperature of 340K.

Covalent organic frameworks (COFs) are materials of growing interest for electronic applications due to their tunable structures, chemical stability, and layered architectures that support extended π-systems and directional charge transport. While their electronic properties are strongly influenced by the choice of molecular building blocks and the stacking arrangement, experimental control over these features remains limited, and the number of well-characterized COFs is still relatively small. Here, we explore two alternative strategies, hydrostatic pressure and metal intercalation, to tune the electronic structure of COFs. Using periodic density functional theory (DFT) calculations, we show that the band gap of pristine COF-1 decreases by ∼1 eV under compression up to 10 GPa. Metal intercalation induces an even greater reduction, in some cases leading to metallic behavior. We demonstrate that pressure and intercalation offer effective, continuous control over COF electronic properties, providing powerful means to complement and extend conventional design approaches.

Augmented plane wave methods enable an efficient description of atom-centered or localized features of the electronic density, circumventing high energy cutoffs and thus prohibitive computational costs of pure plane wave formulations. To complement existing implementations for ground-state properties and excitation energies, we present the extension of the Gaussian and augmented plane wave method to excited-state nuclear gradients within the CP2K program package. Benchmarks for a test set of 35 small molecules demonstrate that maximum errors in the nuclear forces for excited states of singlet and triplet spin multiplicity are smaller than 0.1 eV/Å. The method is furthermore applied to the calculation of the zero-phonon line of defective hexagonal boron nitride. This spectral feature is reproduced with an error of 0.2 eV in comparison to GW-Bethe-Salpeter reference computations and of 0.4 eV in comparison to experimental measurements. Accuracy assessments and applications thus demonstrate the potential use of the outlined developments for large-scale applications on excited-state properties of extended systems.

Nonadiabatic molecular dynamics simulations provide a theoretical understanding of various excited-state processes in photochemistry, offering access to band widths, radiative or nonradiative relaxation and corresponding lifetimes, excited-state energies, and charge transfer. The range of method developments within the framework of time-dependent density functional theory is exceedingly large for molecular quantum chemistry. Still, it shrinks significantly when aiming to treat periodic boundary conditions. To address this gap and complement existing software packages for solid-state nonadiabatic molecular dynamics, we present an interface between the CP2K electronic structure and the NEWTON-X surface hopping codes. The interface features the generation of initial conditions, as well as adiabatic and nonadiabatic molecular dynamics, based on phenomenological or numerical time-derivative couplings. Setups are validated on gas-phase pyrazine, with electronic absorption spectra and excited-state populations for transitions between the lowest singlet states being in agreement with established molecular quantum chemistry methods. Extending the system size to crystalline pyrazine, limitations of approximate couplings are discussed, and the efficiency and applicability of the interface are demonstrated by computing broad spectra over several eV and 100 fs trajectories, considering couplings between all 80th lowest excited states, at low computational cost with a mixed semiempirical density functional theory setup.

We give an overview of the role of “quantum-chemoinformatics” in drug development. Quantum- chemoinformatics is a data-driven chemistry using descriptors on the basis of theoretical chemistry, especially quantum chemistry (QC) and ab initio molecular dynamics (MD) simulations. We focus especially on quantum-chemoinformatics for chemical reaction design and prediction, which is one of the important processes in basic research of drug development. We start with a brief historical overview and then introduces two projects of quantum-cheminformatics. The RMap project uses QC-based chemical reaction route networks for discovery and design of new molecules and reactions. The other project is related to environmental pollution by drug molecules, a property which should be taken into account in drug design and evaluation. The last section describes our recent attempt to accelerate QC-data acquisition by utilizing a limited amount of experimental data and machine learning (ML) technology.

We developed a general framework for hybrid quantum-classical computing of molecular and periodic embedding approaches based on an orbital space separation of the fragment and environment degrees of freedom. We demonstrate its potential by presenting a specific implementation of periodic range-separated DFT coupled to a quantum circuit ansatz, whereby the variational quantum eigensolver and the quantum equation-of-motion algorithm are used to obtain the low-lying spectrum of the embedded fragment Hamiltonian. The application of this scheme to study localized electronic states in materials is showcased through the accurate prediction of the optical properties of the neutral oxygen vacancy in magnesium oxide (MgO). Despite some discrepancies in the position of the main absorption band, the method demonstrates competitive performance compared to state-of-the-art ab initio approaches, particularly evidenced by the excellent agreement with the experimental photoluminescence emission peak.

The Random-Phase approximation (RPA) provides an appealing framework for semi-local density functional theory. In its Resolution-of-the-Identity (RI) approach, it is a very accurate and more cost-effective method than most other wavefunction-based correlation methods. For widespread applications, efficient implementations of nuclear gradients for structure optimizations and data sampling of machine learning approaches are required. We report a well scaling implementation of RI-RPA nuclear gradients on massively parallel computers. The approach is applied to two polymorphs of the benzene crystal obtaining very good cohesive and relative energies. Different correction and extrapolation schemes are investigated for further improvement of the results and estimations of error bars.

Isostructural metal-organic frameworks (MOFs), namely MFU-4 and MFU-4-Br, in which the pore apertures are defined by anionic side ligands (Cl− and Br−, respectively), were synthesized and loaded with noble gases. By selecting the type of side ligand, one can fine-tune the pore aperture size, allowing for precise regulation of the entry and release of gas guests. In this study, we conducted experiments to examine gas loading and release using krypton and xenon as model gases, and we complemented our findings with computational modeling. Remarkably, the loaded gas guests remained trapped inside the pores even after being exposed to air under ambient conditions for extended periods, in some cases for up to several weeks. Therefore, we focused on determining the energy barrier preventing gas release using both theoretical and experimental methods. The results were compared in relation to the types of hosts and guests, providing valuable insights into the gas trapping process in MOFs, as well as programmed gas release in air under ambient conditions. Furthermore, the crystal structure of MFU-4-Br was elucidated using the three-dimensional electron diffraction (3DED) technique, and the bulk purity of the sample was subsequently verified through Rietveld refinement.

We have studied polarized Au(100) and Au(111) electrodes immersed in electrolyte solution by implementing finite-field methods in density functional theory-based molecular dynamics simulations. This allows us to directly compute the Helmholtz capacitance of electric double layer by including both electronic and ionic degrees of freedom, and the results turn out to be in excellent agreement with experiments. It is found that the electronic response of Au electrode makes a crucial contribution to the high Helmholtz capacitance and the instantaneous adsorption of Cl can lead to a charge inversion on the anodic polarized Au(100) surface. These findings point out ways to improve popular semi-classical models for simulating electrified solid-liquid interfaces and to identify the nature of surface charges therein which are difficult to access in experiments.

Ab initio molecular dynamics (AIMD) is an important simulation method applied to understand systems encountered in chemistry and biochemistry, as well as in physics and materials science. It provides means to study complex systems undergoing chemical reactions or phase transitions, as occurring in systems under extreme conditions, or previously unknown compositions. We introduce the basic concepts and methods of AIMD when used together with density functional theory approaches. The connection between the currently most widely used methods, Born–Oppenheimer MD and 2nd Generation Car–Parrinello MD is derived. Important parameters for efficient and accurate simulation protocols are discussed with the focus on the implementations in the CP2K software. Step by step procedures to correctly set up simulations are provided. Sample applications from literature are presented with a special emphasis on the simulation system parameters and protocols used.

Isostructural metal-organic frameworks (MOFs), namely MFU-4 and MFU-4-Br, in which the pore apertures are defined by anionic side ligands (Cl- and Br-, respectively), were synthesized and loaded with noble gases. By selecting the type of side ligand, one can fine-tune the pore aperture size, allowing for precise regulation of the entry and release of gas guests. In this study, we conducted experiments to examine gas loading and release using krypton and xenon as model gases, and we complemented our findings with computational modeling. Remarkably, the loaded gas guests remained trapped inside the pores even after being exposed to air under ambient conditions for extended periods, in some cases for up to several weeks. Therefore, we focused on determining the energy barrier preventing gas release using both theoretical and experimental methods. The results were compared in relation to the types of hosts and guests, providing valuable insights into the gas trapping process in MOFs, as well as programmed gas release in air under ambient conditions. Furthermore, the crystal structure of MFU-4-Br was elucidated using the three-dimensional electron diffraction (3DED) technique, and the bulk purity of the sample was subsequently verified through Rietveld refinement.