Implementation of Girsanov Reweighting in CP2K

TL;DR

Integrating this reweighting framework into the CP2K electronic-structure and molecular-dynamics software package delivers a robust and widely applicable tool for metadynamics, uncertainty quantification, and force-field optimisation based on both \textit{ab initio} and classical MD simulations.

Abstract

Dynamical reweighting of path measures is a powerful approach for accurately evaluating slow molecular processes using modified potential energy surfaces used in enhanced sampling methods. Integrating this reweighting framework into the CP2K electronic-structure and molecular-dynamics (MD) software package delivers a robust and widely applicable tool for metadynamics, uncertainty quantification, and force-field optimisation based on both \textit{ab initio} and classical MD simulations. Based on the Girsanov theorem for stochastic dynamical systems, the method is adapted to the Bussi-Donadio-Parrinello velocity-rescaling scheme. This scheme is accessible through the CSVR thermostat and can be interpreted as a Langevin OVRVO/OBABO update scheme requiring two random numbers per integration step. Comprehensive implementation details are provided, including a complete overview of the modified CP2K modules involved in the MD cycle and the computation of reweighting factors. The framework supports multiple sources of bias potentials, enabling reweighting via the PLUMED interface or CP2K-native external potentials and restraints. The feasibility of the implementation is demonstrated through rerun benchmarks, dynamical reweighting of Markov state models (MSMs) and computation of transport properties based on CP2K trajectories, showing excellent agreement with reference simulations.

Figures (12)

• Figure 1: a) Overlay of congruent update graphs from $(q_{i,k}, p_{i,k})$ to $(q_{i,k+1}, p_{i,k+1})$ for the $\mathcal{O'V'RV'O'}$ integrator, comparing the simulation potential $V(\bm{q}_k)$ with the target potential $\widetilde{V}(\bm{q}_k)$;X001M076 b) Impact of the location of the support point on $\mathcal{N}(0,1)$ for the calculation of $\widetilde{\eta}$: Depending on the sign of $\eta$, a one-step transition can become more or less probable upon reweighting.
• Figure 2: a) The two newly introduced Fortran types that travel through the MD cycle to gather the necessary data for calculating the reweighing factor. b) New Girsanov section as MD subsection. Each PRINT section represents a standard CP2K PRINT accepting keywords such as FILENAME and the EACH subsection for defining the output rate.
• Figure 3: Call diagramm of all CP2K modules involved in Girsanov reweighting using the CSVR thermostat. Each box refers to the file name and MODULE that can be found under the source path given in the upper line. SUBROUTINES are listed below the MODULE name and arrows indicate the CALL command. Light blue boxes correspond to existing modules that have been changed, whereas dark blue boxes correspond to new modules. Calls for setting up, configuring, and releasing the Girsanov environment are marked in orange. Runtime calls updating the Girsanov environment within the MD loop are shown in yellow.
• Figure 4: Simple Rerun Benchmarks: Results of the rerun benchmark of a Lennard-Jones particle forced onto a circular trajectory using CP2K and PLUMED. a) Superposition of the trajectories generated with CP2K applying the bias $-U(\bm{q})$ and the ones generated with the $\mathcal{O'V'RV'O'}$ script using the $\Delta\eta^{(z)}$ provided by the CP2K run; b) Correlation of $\eta^{(z)}$ and $\widetilde{\eta}^{(z)}=\eta^{(z)} + \Delta\eta^{(z)}$; c) Static (blue) and dynamic (orange) reweighting factor increments in their logarithmic form; d-f) Correlation of the positions in $x$, $y$, $z$, respectively, between the original run and the rerun. The color indicates the time step (black: starting point).
• Figure 5: Markov State Model on a 1D-periodic double basin potential: Results of a reweighting study using CP2K and PLUMED. a) 1D-profile of the double well showing the unbiased PES and the Boltzmann distribution as well as the $3RT$ percentile of the kinetic energy (top), 2D-image of the unbiased PES periodically extended into $x$, where white dashed lines mark the simulation cell (bottom); b) same as a, but after application of the bias potential (shown in yellow); c) implied time scales of the slowest process for the biased and unbiased runs as well as after Girsanov reweighting, with d) the corresponding stationary distributions, and e) the first eigenvector.