PaperlayerView exact revision on arXiv

Source-linked AI summary

CP2K: An Electronic Structure and Molecular Dynamics Software Package -- Quickstep: Efficient and Accurate Electronic Structure Calculations

Thomas D. Kühne, Marcella Iannuzzi, Mauro Del Ben, Vladimir V. Rybkin, Patrick Seewald, Frederick Stein, Teodoro Laino, Rustam Z. Khaliullin, Ole Schütt, Florian Schiffmann, Dorothea Golze, Jan Wilhelm, Sergey Chulkov, Mohammad Hossein Bani-Hashemian, Valéry Weber, Urban Borstnik, Mathieu Taillefumier, Alice Shoshana Jakobovits, Alfio Lazzaro, Hans Pabst, Tiziano Müller, Robert Schade, Manuel Guidon, Samuel Andermatt, Nico Holmberg, Gregory K. Schenter, Anna Hehn, Augustin Bussy, Fabian Belleflamme, Gloria Tabacchi, Andreas Glöß, Michael Lass, Iain Bethune, Christopher J. Mundy, Christian Plessl, Matt Watkins, Joost VandeVondele, Matthias Krack, Jürg Hutter

arXiv:2003.03868v2physics.chem-phcond-mat.mtrl-sciphysics.comp-ph

TL;DR

Atomistic simulations require efficient methods for challenging extended systems, motivating CP2K’s broad electronic-structure and molecular-dynamics capabilities. This review presents CP2K’s scalable methods and demonstrates efficient, accurate simulations across electronic-structure and ab-initio molecular-dynamics applications.

  • Problem

    Challenging atomistic simulations require computational methods that efficiently handle extended systems and provide broad methodological flexibility.

    Source: §I · p. 2

  • Method

    The review surveys CP2K’s density-functional-theory and post-Hartree–Fock methods, including Gaussian-and-plane-wave and augmented all-electron approaches.

    Source: §I · p. 2

  • Results

    CP2K combines efficient algorithms and parallel scalability with broad capabilities, including early-offset linear scaling and substantially accelerated molecular-dynamics workflows.

    Source: §I · p. 2, B. Linear scaling methods based on… · p. 24, B. Second-generation Car-Parrinello… · p. 27

  • Takeaways & Limitations

    CP2K supports flexible, efficient atomistic simulations across extended condensed-phase systems and diverse computational workflows.

    Source: §I · p. 2, E. CP2K workflows · p. 43

    The IC-QM/MM approach is appropriate when molecule–molecule interactions are primary and the surface mainly acts as a template.

    Source: §3 · p. 35

Abstract

from arXiv · show

CP2K is an open source electronic structure and molecular dynamics software package to perform atomistic simulations of solid-state, liquid, molecular and biological systems. It is especially aimed at massively-parallel and linear-scaling electronic structure methods and state-of-the-art ab-initio molecular dynamics simulations. Excellent performance for electronic structure calculations is achieved using novel algorithms implemented for modern high-performance computing systems. This review revisits the main capabilities of CP2k to perform efficient and accurate electronic structure simulations. The emphasis is put on density functional theory and multiple post-Hartree-Fock methods using the Gaussian and plane wave approach and its augmented all-electron extension.

I. INTRODUCTION · II. GAUSSIAN AND PLANE WAVES METHOD · A. Kohn–Sham energy, forces, and stress tensor

CP2K provides scalable atomistic simulation through broad electronic-structure and molecular-dynamics capabilities, centered on the Gaussian-and-plane-waves framework. Quickstep combines Gaussian orbitals with plane-wave density representations to support efficient energies, forces, and stress calculations.

I. INTRODUCTION

  • CP2K combines efficient algorithms and strong parallel scalability to support large-scale ab-initio molecular dynamics and flexible method selection for extended condensed-phase systems.Its simulations provide atomistic insight into chemical, physical, biological, and materials systems that may be difficult or impossible to investigate experimentally.Source: §I · pp. 1–2

II. GAUSSIAN AND PLANE WAVES METHOD

  • Quickstep supports semi-empirical, tight-binding, orbital-free, Kohn–Sham DFT, MP2, dRPA, and GW methods through shared integral and optimization infrastructure.The review emphasizes CP2K’s Gaussian-and-plane-waves methodology as a distinguishing approach among electronic-structure programs.Source: §II · p. 2
  • The GPW method represents orbitals with contracted Gaussian basis functions while expanding the electron density on an equidistant grid or equivalently in plane waves.Gaussian Fourier transforms and exponential grid convergence enable this dual representation; expansion accuracy is controlled by grid spacing or the plane-wave cutoff.Source: §II · p. 2
  • Reciprocal lattice vectors define the plane-wave basis for periodic systems, allowing real-space functions on grid points to be transformed into reciprocal-space representations.The reciprocal vectors derive from the computational-cell matrix and lattice vectors.Source: §II · p. 2
  • The density-matrix-to-plane-wave mapping uses multigrid methods, real-space screening, and Cartesian-Gaussian separability to achieve linear scaling with a small prefactor.This dual density representation facilitates efficient and accurate electrostatic, exchange, and correlation calculations.Source: §II · p. 3

A. Kohn–Sham energy, forces, and stress tensor

  • Within the GPW framework, the Kohn–Sham energy includes kinetic, ionic, electrostatic, and exchange-correlation terms, with optional Brillouin-zone k-point sampling.Electrostatics are evaluated using an Ewald method, with Gaussian compensation charges accounting for self, overlap, and related interactions.Source: A. Kohn–Sham energy, forces, and… · p. 3
  • GPW computes nuclear forces consistently, including Pulay contributions, while Fourier-transform methods and virial expressions provide the internal stress tensor.Special treatment is required for cell-dependent grid integration in GGA exchange-correlation contributions; the same mapping machinery also supports kinetic-energy-density functionals.Source: A. Kohn–Sham energy, forces, and… · pp. 3–4

B. Dual-space pseudopotentials … E. Gaussian-augmented plane waves approach

The GPW framework combines analytic dual-space pseudopotentials, transferable Gaussian basis sets, and local density fitting to improve accuracy, conditioning, and scaling. GAPW further separates smooth and atom-centered density components, enabling all-electron-sensitive properties while introducing basis-dependent approximation limits.

B. Dual-space pseudopotentials

  • Dual-space pseudopotentials use analytic Gaussian projectors that are transferable and norm-conserving, avoiding the rapidly increasing plane-wave cutoffs associated with Gaussian exponents.Their parameters are optimized against relativistic all-electron atomic calculations, with parameter databases available for different exchange-correlation potentials.Source: B. Dual-space pseudopotentials · p. 4

C. Basis sets

  • The molopt basis sets target accuracy, systematic improvement, broad environmental transferability, well-conditioned overlaps, and efficient treatment of diffuse functions.They use generally contracted functions with shared exponents and globally optimized parameters, supporting molecular, liquid, and dense-solid calculations.Source: C. Basis sets · p. 4
  • Reducing Gaussian primitives and removing very diffuse functions yields a 10-fold time reduction for mapping routines in liquid-water calculations with default settings, while retaining good performance in many environments.The smallest oxygen exponent in the basis set is approximately 0.162.Source: C. Basis sets · pp. 4–5

D. Local density fitting approach

  • LRIGPW replaces pair-based grid mapping with local auxiliary-Gaussian density fitting, reducing the 64-water example from ≈200000 atom pairs to 192 atom terms.The method fits pair densities while preserving electron number, and CP2K accelerates its demanding overlap integrals with solid-harmonic Gaussian functions.Source: D. Local density fitting approach · p. 5
  • 99% of atom pairs in the water-box example can be treated as distant pairs through single-atom expansions, substantially reducing LRIGPW memory and computational cost.A distance criterion and switching function preserve a smooth potential energy.Source: D. Local density fitting approach · p. 5

E. Gaussian-augmented plane waves approach

  • GAPW separates smooth interatomic density from rapidly varying atom-centered components, evaluating plane-wave contributions with FFTs and local terms on atom-centered grids.The decomposition is systematically improvable with larger local basis sets, but its accuracy is tied to the primitive Gaussians in the primary basis.Source: E. Gaussian-augmented plane waves… · pp. 5–6
  • GAPW enables properties depending on core electrons, including liquid-water X-ray scattering, X-ray absorption spectra, and nuclear and electronic magnetic properties.The approach is presented as an alternative to pseudopotentials or as a route to smaller cutoffs in pseudopotential calculations.Source: E. Gaussian-augmented plane waves… · pp. 5–6

III. HARTREE-FOCK AND HYBRID DENSITY FUNCTIONAL THEORY METHODS … 3. Implementation

CP2K extends DFT with Hartree-Fock exchange and hybrid functionals, using screening, parallel in-core integration, multiple-time-step evaluation, and ADMM to enable large condensed-phase simulations. It also implements correlated post-Hartree-Fock methods, including MP2 variants with RI, Laplace, sparse, massively parallel, and GPU-accelerated algorithms.

III. HARTREE-Fック AND HYBRID DENSITY FUNCTIONAL THEORY METHODS

  • Hybrid functionals replace part of DFT exchange with exact Hartree-Fock exchange, enabling condensed-phase AIMD, electrochemistry, and correlated methods such as MP2, RPA, and GW.HFX was introduced in CP2K for relatively large, disordered systems and provides a foundation for subsequent correlated wavefunction methods.Source: §III · p. 6

III. HARTREE-FOCK AND HYBRID DENSITY FUNCTIONAL THEORY METHODS

  • HFX screening reduces computational scaling from O(N 4) to O(N 2), while MPI/OMP-parallel in-core integral storage, compression, and multiple-time-step evaluation reduce practical AIMD costs.The multiple-time-step approach evaluates HFX energies only every few time steps, assuming the GGA-to-hybrid potential-energy difference varies slowly.Source: §III · p. 7
  • ADMM approximates the density matrix in an auxiliary basis and corrects the resulting exchange energy with a GGA-style functional, often bringing HFX within a few times GGA cost.The method enabled over 300 ps of AIMD for 160 water molecules and spin-density calculations for a solvated metallo-protein of approximately 3000 atoms.Source: §III · p. 7

IV. BEYOND HARTREE-FOCK METHODS

  • CP2K also provides post-Hartree-Fock methods that are more accurate than hybrid DFT, including MP2, RI-MP2, Laplace-transformed MP2, and SOS-MP2.These methods are available within CP2K as correlated wavefunction approaches beyond Hartree-Fock.Source: §IV · p. 7, §3 · p. 8

1. Theory

  • MP2 captures most dynamic electron correlation, while RI-MP2 retains O(N 5) scaling with a lower prefactor by reducing integral-computation cost.RI replaces transformed electron-repulsion integrals with approximated integrals using auxiliary basis functions and two-center quantities.Source: §1 · pp. 7–8

2. Scaled opposite-spin MP2

  • SOS-MP2 neglects same-spin correlation, rescales the opposite-spin contribution, and uses Laplace quadrature requiring approximately 7 points for micro-Hartree accuracy.Its most expensive matrix-element calculation scales as O(N 4).Source: §2 · p. 8

3. Implementation

  • The massively parallel MP2 implementation exploits sparse algebra and GPU-accelerated matrix multiplication; gradients for hundreds of atoms can run on thousands of CPU cores.Gradient evaluation costs 4–5 times the energy calculation, while open-shell MP2 is typically 3–4 times more expensive than closed-shell MP2.Source: §3 · p. 8

4. Applications … C. Ionization potentials and electron affinities from GW

CP2K implements scalable RPA, SOS-MP2, and GW methods for large systems, while applications span aqueous simulations, periodic water benchmarks, and molecular spectroscopy. Its GW approaches provide accurate ionization potentials and electron affinities and support efficient level-alignment calculations for molecules on metal surfaces.

4. Applications

  • CP2K’s RI-MP2 enabled aqueous simulations that reproduced liquid-water density, structure, and IR spectra using 64-water-molecule simulation cells.Applications also included ice XV structure refinement, bulk hydrated electrons in 47-water-molecule cells, and a condensed-phase radical simulation.Source: §4 · p. 8

B. Random Phase Approximation Correlation Energy Method

  • RPA correlation energies in CP2K use a dielectric-matrix formulation with Gaussian RI auxiliary bases, minimax frequency integration, O(N 4Nq) cost, and O(N 3) storage.The conventional EXX/RPA formalism combines DFT, exact exchange, and RPA correlation evaluations using DFT orbitals and orbital energies.Source: B. Random Phase Approximation… · pp. 8–9, X + ERPA · p. 9

1. Implementation of the quartic scaling RPA and SOS-MP2 methods

  • Quartic-scaling RPA and SOS-MP2 share the same computational cost for large systems when they use the same number of frequency or time quadrature points.Their common bottleneck is distributed matrix multiplication over independent quadrature points, enabling excellent parallel scalability to several thousand nodes.Source: §1 · p. 9

2. Cubic Scaling RPA and SOS-MP2 method

  • CP2K’s cubic-scaling RPA/SOS-MP2 scheme has effective O(N 3) or better scaling for all systems, extending RPA to systems containing thousands of atoms.Sparse tensor contractions and localized Gaussian basis functions reduce the dominant computational work, while large systems can benefit from sparse density matrices.Source: §2 · pp. 9–10
  • At 256 water molecules, cubic-scaling RPA is one order of magnitude faster than quartic RPA, with effective scaling of O(N 1.8).The cubic implementation becomes favorable near 90 water molecules and maintains relative RPA-energy errors below 0.01%.Source: §2 · p. 10, Figure 1 · p. 10

C. Ionization potentials and electron affinities from GW

  • GW has mean absolute deviations below 0.2 eV from CCSD(T) for IPHOMO and EA, while including vibrational effects reduces deviation from experiment to < 0.1 eV. Depending on the DFT functional, evGW improves HOMO-LUMO gaps relative to G0W0 by 0.1–0.3 eV.Source: C. Ionization potentials and electron… · pp. 11–12
  • The GW implementation routinely treats molecules with standard O(N 4) scaling, while its low-scaling variant effectively reaches O(N 2) to O(N 3) for molecules exceeding 1,000 atoms.Standard G0W0 targets valence orbital energies and uses analytic-continuation frequency integration through two-pole or Padé approximations.Source: C. Ionization potentials and electron… · p. 11
  • For molecules on metal surfaces, CP2K combines gas-phase evGW with an image-charge correction to compute HOMO-LUMO gaps while accounting for surface-induced level shifts.The image charge raises occupied states, lowers unoccupied states, and significantly reduces the molecular gap compared with the gas phase.Source: Figure 2 · p. 11, C. Ionization potentials and electron… · p. 12

V. DENSITY FUNCTIONAL PERTURBATION THEORY · A. Polarizability

DFPT obtains response properties as energy or density derivatives by perturbatively minimizing an extended energy functional, generally outperforming direct numerical perturbations. For polarizability, electric-field coupling is formulated through Berry-phase polarization in periodic systems and supports Raman-spectrum calculations.

V. DENSITY FUNCTIONAL PERTURBATION THEORY

  • DFPT computes response properties from perturbative energy and density derivatives, offering an alternative that performs better than directly perturbing systems and recomputing total energies.The perturbation is introduced through an external potential or functional, and the electronic response is evaluated variationally.Source: §V · p. 12
  • The energy, orbitals, and density are expanded in powers of λ, with the ground-state Kohn–Sham equations providing the zero-order solution and first-order orbitals determining second-order energy.The second-order energy is variational in the first-order wavefunction, while the energy kernel requires second-order exchange-correlation functional derivatives.Source: §V · p. 12
  • The first-order perturbation problem is solved directly with a preconditioned conjugate-gradient minimization, while orthonormality enforces conservation of total charge.The perturbation density retains dependence on the first-order orbitals, making the response equations self-consistent.Source: §V · pp. 12–13

A. Polarizability

  • For external electric fields, periodic systems use the modern theory of polarization and Berry phases because the position operator is ill-defined.The electric field couples to the electronic polarization rather than being represented directly by an ordinary Hamiltonian perturbation.Source: A. Polarizability · p. 13
  • Polarizability tensor elements are derivatives of electronic polarization with respect to electric field, measuring deformation of a molecule’s electron cloud.Changes in polarizability accompany Raman-active vibrations and define isotropic and anisotropic Raman transition polarizabilities.Source: A. Polarizability · p. 13
  • The arbitrary-functional formulation also supports orbital-specific perturbations, whose stationary condition produces an inhomogeneous nonlinear system for first-order Kohn–Sham orbitals.The simulation-cell matrix and electric-field coupling determine the perturbation and induced electric dipole.Source: A. Polarizability · p. 13
  • During AIMD, time-dependent polarizability enables depolarized Raman intensities from polarizability autocorrelation functions, alongside vibrational and infrared spectra from velocity and dipole autocorrelations.The relevant spectra are obtained by temporal Fourier transformation of the corresponding autocorrelation functions.Source: A. Polarizability · p. 13

B. Nuclear magnetic resonance and electron paramagnetic resonance spectroscopy

CP2K’s all-electron GAPW/DFPT linear-response framework computes NMR chemical shifts and, for spin-1/2 systems, EPR g-tensors from induced current densities. The implementation addresses periodic systems with localized Wannier-based operators and supports relativistic corrections and gauge choices relevant to accuracy and computational cost.

B. Nuclear magnetic resonance and electron paramagnetic resonance spectroscopy

  • The all-electron GAPW/DFPT treatment derives NMR chemical shifts from induced current densities and computes EPR g-tensors for systems with net electronic spin 1/2.The induced current density generated by a static magnetic perturbation determines the NMR chemical shift at each nucleus and the EPR g-tensor.Source: B. Nuclear magnetic resonance and… · p. 14
  • The spin-orbit correction is the leading relativistic contribution to EPR g-tensors, becomes more important for heavy elements, and is evaluated through spin-dependent current densities and effective-potential gradients.CP2K replaces computationally demanding two-electron integrals with an effective one-electron operator.Source: B. Nuclear magnetic resonance and… · p. 14
  • For periodic systems, CP2K uses maximally localized Wannier functions and localized position operators so the total induced current remains gauge-independent despite gauge-dependent diamagnetic and paramagnetic components.The induced current is separated into diamagnetic and paramagnetic terms, while the Wannier-based position operator avoids the invalid multiplicative position operator in periodic systems.Source: B. Nuclear magnetic resonance and… · p. 14
  • With real ground-state molecular orbitals, the unperturbed current and first-order density response vanish, allowing CP2K to skip the second-order energy kernel when solving the linear-response equations.The implementation uses orbital-dependent angular-momentum, momentum, and full-correction operators, with Wannier centers defining individual reference systems.Source: B. Nuclear magnetic resonance and… · pp. 14–15
  • Gauge choice strongly affects basis-set convergence: CSGT eliminates the diamagnetic current but requires rich basis sets near nuclei, whereas IGAIM places the gauge at a nearby nuclear center and enables large-scale paramagnetic-solid NMR analyses.Paramagnetic NMR contributions can be combined from hybrid-functional hyperfine couplings with gradient-corrected-functional g-tensors and orbital shieldings.Source: B. Nuclear magnetic resonance and… · p. 15

VI. TIME-DEPENDENT DENSITY FUNCTIONAL THEORY … A. Traditional diagonalization

CP2K supports time-dependent electronic-structure simulations through linear-response and real-time TDDFT, including Ehrenfest dynamics, while offering diagonalization-based eigensolvers for SCF calculations. These methods balance accuracy, scalability, and applicability across system types, basis sets, and occupation schemes.

VI. TIME-DEPENDENT DENSITY FUNCTIONAL THEORY

  • CP2K’s time-dependent DFT methods investigate many-body dynamics and properties under time-dependent potentials, including electric and magnetic fields.TDDFT provides the framework for studying these driven systems.Source: §VI · p. 15

A. Linear-response time-dependent density functional theory

  • Linear-response TDDFT computes vertical transition energies and oscillator strengths from the response to a weak electromagnetic perturbation.Its implementation uses the Tamm-Dancoff and adiabatic approximations, reducing the equations to a Hermitian eigenproblem.Source: A. Linear-response time-dependent… · p. 15
  • A block Davidson solver supports hybrid exchange, integral screening, truncated Coulomb operators, ADMM, separate ground-state and kernel functionals, and simultaneous response-vector computations.Simultaneous evaluation improves load balancing and reduces communication costs, enabling effective use of more CPU cores.Source: A. Linear-response time-dependent… · p. 16

1. Applications

  • LR-TDDFPT performance was comparable to ground-state calculations for systems of approximately 1000 atoms, sufficient to converge localized defect transitions in medium- to wide-band-gap materials.Applications included periodic imogolite nanotube, MgO, and HfO2 defect systems, while varying Hartree-Fock exchange and verifying ADMM accuracy.Source: §1 · p. 16

B. Real-time time-dependent density functional theory and Ehrenfest dynamics

  • Real-time TDDFT captures nonlinear effects and electron-driven dynamics, while Ehrenfest dynamics propagates electrons and nuclei simultaneously when electron–nuclear coupling matters.CP2K implements both cubic-scaling MO-RTP and linear-scaling density-matrix P-RTP forms.Source: B. Real-time time-dependent density… · p. 16
  • Pure linear-scaling Ehrenfest dynamics becomes cubic with system size because density-matrix decay can fail, but subsystem DFT can restore true linear scaling across subsystems.The Kim–Gordon approach minimizes subsystems independently and preserves block-diagonal overlap structure under P-RTP.Source: B. Real-time time-dependent density… · p. 17

VII. DIAGONALIZATION-BASED AND LOW-SCALING EIGENSOLVER

  • The listed methods are traditional diagonalization, pseudo diagonalization, orbital transformation, and purification; OT is most efficient and scalable for growing systems but requires fixed integer MO occupations.TD supports fractional MO occupations but its cubic scaling limits accessible system sizes; PD may provide significant speedups.Source: §VII · p. 17

A. Traditional diagonalization

  • CP2K solves its non-orthogonal Gaussian-basis eigenproblem using ScaLAPACK or ELPA, with Cholesky or symmetric Löwdin orthogonalization; Löwdin diagonalization detects and filters linear dependencies.Overlap-matrix eigenvalues smaller than 10−5 usually indicate significant linear dependencies.Source: A. Traditional diagonalization · p. 17

1. TD/DIIS … C. Orbital transformations

CP2K accelerates SCF convergence using TD combined with DIIS, while addressing metallic-system charge sloshing and offering pseudo-diagonalization and direct orbital-energy minimization as alternatives. These methods exploit occupied–unoccupied structure, variational formulations, and efficient orthogonality-preserving updates to reduce computational effort and improve robustness.

1. TD/DIIS

  • DIIS accelerates SCF convergence through the KS–density-matrix commutator error and is especially efficient from a sufficiently pre-converged density.Its iteration savings matter when Hamiltonian construction costs more than diagonalization.Source: §1 · p. 18

2. TD/Broyden and Kerker mixing

  • DIIS can struggle for metallic systems because imbalanced short- and long-range charge redistribution causes charge sloshing.Source: §2 · p. 18

B. Pseudo diagonalization

  • Pseudo-diagonalization transforms the KS matrix into the preceding step’s MO basis, where the matrix becomes diagonally dominant and occupied–unoccupied block structure enables wavefunction convergence.Only the occupied–unoccupied coupling block must be eliminated because the total energy depends on occupied MOs.Source: B. Pseudo diagonalization · p. 18
  • Jacobi rotations iteratively decouple the occupied and unoccupied blocks while preserving MO orthonormality and using computationally cheap BLAS level 1 operations.The relevant block contains only 10–20% of the KMO matrix, and its nonzero blocks decrease rapidly near convergence.Source: B. Pseudo diagonalization · p. 18

C. Orbital transformations

  • Direct minimization of the electronic energy functional is more robust in principle, can replace diagonalization with fewer matrix-matrix multiplications, and reduces time-to-solution for difficult large systems.Preconditioners are often used to promote faster and smoother convergence.Source: C. Orbital transformations · pp. 18–19
  • Orbital transformations formulate ground-state calculation as constrained minimization of E[C] over occupied orbitals satisfying C^TSC = 1.The minimized functional depends on the electronic-structure theory; in hybrid Hartree–Fock/DFT it uses core, Coulomb, exact-exchange, and exchange-correlation terms.Source: C. Orbital transformations · p. 19
  • Efficiently enforcing orbital orthogonality is a major challenge in direct minimization, motivating two approaches implemented in CP2K.Source: C. Orbital transformations · p. 19

1. Orthogonality constraints · 2. Minimizer · 3. Preconditioners

Quickstep enforces orbital orthogonality through orbital-transformation parametrizations, minimizes the energy with safeguarded iterative methods, and accelerates convergence using scalable preconditioners. These components balance robustness, accuracy, and computational cost across large simulations.

1. Orthogonality constraints

  • Orbital transformation methods convert nonlinear orbital orthogonality constraints into unconstrained or linearly constrained formulations that preserve C^TSC = 1.OT/Diag and OT/Taylor evaluate matrix functions by diagonalization or truncated Taylor expansion, while OT/IR uses refinement functions accurate to finite order.Source: §1 · p. 19
  • OT/IR approximates the constraint function recursively through iterative refinement expansions, extending accuracy to any finite order.The approximation is constructed around Z0 and uses refinement functions for an approximate inverse matrix factorization.Source: §1 · p. 19

2. Minimizer

  • CP2K safeguards DIIS by switching to steepest descent when a DIIS step points uphill, improving robustness while retaining DIIS’s typically fast convergence.DIIS minimizes the norm of a linear combination of gradient error vectors and reduces to steepest descent for the initial step.Source: §2 · pp. 19–20
  • Nonlinear conjugate gradients provide a robust, efficient, and numerically stable minimizer using Polak-Ribiere updates with restart and approximate line searches.Quadratic interpolation is generally faster than the more robust golden-section line search in Quickstep.Source: §2 · p. 20
  • Quickstep implements a quasi-Newton minimizer with a Broyden type 2 approximate-inverse-Hessian update and adaptive curvature estimation.Quasi-Newton methods avoid the expense of forming the full Hessian, whose Newton iteration can diverge far from the solution without line search or backtracking.Source: §2 · p. 20

3. Preconditioners

  • Preconditioning approximates the constrained energy curvature with a single positive-definite matrix, because orbital-specific preconditioners are impractical.FULL ALL shifts orbital-dependent eigenvalues and ensures positive definiteness, but requires diagonalization with O(N^3) scaling; lower-scaling alternatives use sparse S and T matrices.Source: §3 · pp. 20–21
  • FULL SINGLE INVERSE often offers the best quality-cost tradeoff by inverting occupied eigenvalues, emphasizing orbitals near the band gap without diagonalization.Its construction scales as O(NM^2) in the dense case, matching the complexity of the orbital-transformation algorithm.Source: §3 · p. 21
  • Sparse inversion and Hotelling iterations reduce preconditioner costs, but aggressive filtering can eventually destabilize iterations; previous inverses provide efficient initial guesses for MD and geometry optimization.The preconditioner is symmetric positive definite, and only a few iterations are needed after the initial inverse approximation when successive matrices change little.Source: §3 · p. 21

D. Purification methods … A. Localization of orthogonal and non-orthogonal molecular orbitals

CP2K provides linear-scaling density-matrix methods based on purification, sign-function iterations, and submatrix calculations. It also supports orthogonal and nonorthogonal localized molecular orbitals for electronic-structure analysis and large-system methods.

D. Purification methods

  • Purification maps KS-matrix eigenvalues to density-matrix occupations and supports linear-scaling calculations for sparse matrices without explicitly forming orbitals.CP2K implements TRS4, TC2, SIGN, and a PEXSI interface; canonical ensembles may require chemical-potential adjustments or extra iterations.Source: D. Purification methods · p. 21

E. Sign-Method

  • CP2K uses Newton-Schulz and higher-order Padé-based sign iterations for linear-scaling matrix inversions, inverse square roots, and density-matrix purification.Implementations reach seventh order, while the fifth-order iteration uses four matrix multiplications per iteration and truncates small elements to preserve sparsity.Source: E. Sign-Method · pp. 21–22
  • Sign-function calculations for the STMV virus show strong wall-time sensitivity to the matrix truncation threshold, motivating combinations with approximate computing.The benchmark uses periodic CP2K GFN2-xTB calculations on a system exceeding one million atoms.Source: E. Sign-Method · p. 22, Figure 3 · p. 22

F. Submatrix Method

  • The submatrix method obtains an approximate density matrix from small, relatively dense principal submatrices whose size is independent of overall system size, yielding linear scaling.Accuracy and computational cost depend on the SCF truncation threshold, and purification can use either sign iterations or direct eigendecomposition.Source: F. Submatrix Method · p. 22

VIII. LOCALIZED MOLECULAR ORBITALS

  • Localized molecular orbitals support visualization of bonding, electronic-structure interpretation, and local electronic-structure methods that reduce the cost of large-system calculations.They are also used in XAS modeling and dispersion-corrected exchange-correlation functionals based on atomic polarizabilities.Source: §VIII · pp. 22–23

A. Localization of orthogonal and non-orthogonal molecular orbitals

  • CP2K localizes orthogonal orbitals with Resta-Berghold and Pipek-Mezey functionals, with the former applicable to gas-phase and periodic systems but limited to Γ-point states in periodic cells.In gas-phase systems, the Resta-Berghold functional is equivalent to Boys-Foster localization; Pipek-Mezey preserves σ–π bond separation and is common for molecules.Source: A. Localization of orthogonal and… · p. 23
  • CP2K’s nonorthogonal localization minimizes orbital spread while penalizing linear dependence, producing nonorthogonal orbitals that are more localized than conventional orthogonal orbitals.The penalty is zero for orthogonal orbitals and diverges for linearly dependent nonorthogonal orbitals, converting the procedure into unconstrained optimization.Source: A. Localization of orthogonal and… · p. 23, Figure 5 · p. 23

B. Linear scaling methods based on localized one-electron orbitals … A. Born-Oppenheimer molecular dynamics

CP2K develops localized-orbital and machine-learning basis-set methods to reduce electronic-structure cost while controlling accuracy, and supports AIMD through on-the-fly KS-DFT dynamics. Its BOMD formulation emphasizes complete force terms because forces are more sensitive to electronic-density errors than total energies.

B. Linear scaling methods based on localized one-electron orbitals

  • Localized orbitals reduce variational degrees of freedom to occupied states, but explicit localization constraints are required because unconstrained orbitals delocalize and lose O(N) scaling.Each orbital is assigned a center and radius Rc and expanded only in localized Gaussian basis functions within that radius.Source: B. Linear scaling methods based on… · p. 24
  • CLMO accuracy depends critically on Rc because locality removes distant donor-acceptor interactions while retaining electrostatics, exchange, polarization, and dispersion.The approach is effective for gapped systems when Rc includes nearest and possibly next-nearest neighbors, but is impractical for metals and semimetals.Source: B. Linear scaling methods based on… · p. 24
  • CP2K addresses slow CLMO optimization with a two-stage SCF method for weakly interacting systems and an approximate-Hessian method for strongly covalent systems.The two-stage method first optimizes Rc = 0 orbitals and then relaxes neighboring delocalization, whereas the Hessian approach removes low-curvature occupied-occupied mixing modes.Source: B. Linear scaling methods based on… · p. 24
  • CLMO-based O(N) DFT enables early-offset linear scaling and substantial computational savings without compromising accuracy, including for condensed-phase systems and large basis sets.The localization radius Rc can be tuned to balance computational performance and accuracy.Source: B. Linear scaling methods based on… · p. 24

C. Polarized atomic orbitals from machine learning

  • PAO-ML predicts environment-adaptive polarized atomic-orbital transformations and preserves rotational invariance through neighboring-atom-anchored potentials, enabling analytic forces for AIMD.Training requires optimizing the transformation matrix and electronic density under the Pauli principle.Source: C. Polarized atomic orbitals from… · p. 25
  • Up to 200-fold lower computational cost and 4 orders of magnitude fewer flops were achieved for liquid-water AIMD with PAO-ML, while structural properties fairly matched basis-set-converged results.A very small training set was already satisfactory because the method is variationally robust.Source: C. Polarized atomic orbitals from… · p. 25

IX. AB-INITIO MOLECULAR DYNAMICS

  • AIMD evaluates thermal observables as temporal averages under ergodicity, with CP2K calculating the potential energy on the fly using KS-DFT.CP2K provides two AIMD flavors, including BOMD.Source: §IX · p. 26

A. Born-Oppenheimer molecular dynamics

  • In BOMD, orbitals are minimized at every step under orthonormality constraints, and accurate forces require Hellmann-Feynman, Pulay, and non-self-consistent terms.Neglecting constraint or non-self-consistent forces gives only perturbative generalized forces; forces depend linearly on charge-density errors, unlike energies.Source: A. Born-Oppenheimer molecular dynamics · p. 26

B. Second-generation Car-Parrinello molecular dynamics · C. Low-cost linear-scaling ab-initio molecular dynamics based on compact localized molecular orbitals · D. Multiple-time-step integrator

The paper presents second-generation CPMD, CLMO-based linear-scaling AIMD, and a reversible multiple-time-step integrator for efficient, accurate molecular dynamics. These methods avoid costly electronic iterations, support larger systems or time steps, and address force or sampling errors through damping and force splitting.

B. Second-generation Car-Parrinello molecular dynamics

  • Second-generation CPMD combines BOMD-sized time steps with CPMD efficiency, requiring only one preconditioned gradient calculation and yielding one-to-two-order-of-magnitude gains across diverse systems.The method avoids fictitious electronic masses and has been applied to molecules, solids, phase-change materials, aqueous systems, and interfaces.Source: B. Second-generation Car-Parrinello… · p. 27
  • The electronic dynamics bypasses SCF cycles and diagonalization while approaching the instantaneous ground state within one electronic-force calculation and remaining close to the Born–Oppenheimer surface.A predictor-corrector scheme preserves orthonormality and is time-reversible up to O(∆t2K−2), while allowing ∆t as large as in BOMD.Source: B. Second-generation Car-Parrinello… · pp. 27–28
  • Second-generation CPMD introduces slight nuclear dissipation, which can be corrected with a Langevin-type equation and an automatically bootstrapped or adaptively adjusted friction coefficient for canonical sampling.The friction coefficient can be tuned on-the-fly using a Berendsen-like algorithm or adjusted through adaptive Langevin dynamics.Source: B. Second-generation Car-Parrinello… · p. 28

C. Low-cost linear-scaling ab-initio molecular dynamics based on compact localized molecular orbitals

  • CLMO DFT enables low-overhead AIMD with computational complexity linear in the number of molecules, but nonvariational orbital errors can accumulate during long trajectories and impair sampling.CP2K compensates for neglected nonvariational force components using a modified Langevin equation tuned for stable dynamics with imperfect forces.Source: C. Low-cost linear-scaling ab-initio… · p. 28
  • CLMO AIMD reproduces multiple structural and dynamical properties of liquid water with Rc = 1.6 vdWR and ϵSCF as high as 10−2 a.u., while maintaining stable dynamics.These settings provide low computational overhead for routine medium-scale simulations and preserve linear scaling for larger molecular systems.Source: C. Low-cost linear-scaling ab-initio… · p. 29, Figure 7 · p. 29
  • The CLMO AIMD methodology cannot currently be combined with Harris-functional or perturbative CLMO methods, and its extension to strongly interacting atoms remains underway.Source: C. Low-cost linear-scaling ab-initio… · p. 29

D. Multiple-time-step integrator

  • AIMD-MTS uses reversible r-RESPA force splitting to evaluate accurate hybrid-DFT forces less frequently than cheaper local-functional forces while preserving time reversibility and energy conservation.The accurate-minus-approximate force is usually smaller and low-frequency, enabling distinct inner and outer integration loops.Source: D. Multiple-time-step integrator · pp. 29–30

X. ENERGY DECOMPOSITION AND SPECTROSCOPIC ANALYSIS METHODS … C. Mode selective vibrational analysis

CP2K provides EDA, NMA, and MSVA methods to analyze intermolecular bonding and vibrational spectra. Its MSVA implementation reduces force evaluations for selected modes while supporting large condensed-matter applications.

A. Energy decomposition analysis based on compact localized molecular orbitals

  • CP2K’s ALMO EDA separates total interaction energy into frozen-density, polarization, and charge-transfer terms for gas-phase and condensed-matter systems.The implementation supports closed-shell fragments, fractionally occupied ALMOs, spatially controlled charge transfer, and AIMD analyses with charge transfer switched off.Source: A. Energy decomposition analysis based… · pp. 30–31
  • ALMO EDA studies connect donor–acceptor interactions with X-ray absorption, infrared, and nuclear magnetic resonance spectra, including hydrogen bonding in water systems.AIMD extensions found that the insignificant covalent hydrogen-bond component influences liquid-water structure, dielectric constant, hydrogen-bond lifetime, self-diffusion, and viscosity.Source: A. Energy decomposition analysis based… · p. 31

B. Normal mode analysis of infrared spectroscopy

  • CP2K obtains infrared spectra through Born–Oppenheimer normal-mode analysis by diagonalizing a mass-weighted Hessian built from three-point force differences.For M particles, the central-difference Hessian requires 6M force evaluations and yields 3M normal-mode coordinates and frequencies.Source: B. Normal mode analysis of infrared… · p. 31

C. Mode selective vibrational analysis

  • MSVA reduces computational cost for selected vibrations by iteratively evaluating a Hessian subspace with the Davidson algorithm instead of constructing the full Hessian.Hessian actions are obtained from analytical forces and numerical force derivatives along collective displacement vectors.Source: C. Mode selective vibrational analysis · pp. 31–32
  • After at most 3M iterations, MSVA recovers the numerical-difference Hessian and therefore exact frequencies and normal modes within that approximation.The method requires fewer force evaluations before reaching this limit.Source: C. Mode selective vibrational analysis · p. 32
  • MSVA can target a single mode, frequency range, or atom-localized mode, but convergence to the desired mode depends on a suitable initial guess.Poor guesses can cause premature convergence to another mode or track the nearest mode rather than the intended one.Source: C. Mode selective vibrational analysis · p. 32
  • MSVA parallelizes through distributed force calculations or block Davidson environments, enabling vibrational analysis of large condensed-matter systems.Dipole derivatives must be computed with NMA to obtain IR intensities; dye–anatase simulations helped assign spectra and revealed protonation-dependent binding and self-assembly effects.Source: C. Mode selective vibrational analysis · p. 32

XI. EMBEDDING METHODS … 2. QM/MM for periodic systems

CP2K combines empirical, DFT-based, and quantum-chemical potential-energy methods through flexible embedding and QM/MM formulations. Its multi-grid GEEP implementation accelerates QM/MM electrostatics for isolated and periodic systems, while periodic treatments address long-range replica and QM-image interactions.

XI. EMBEDDING METHODS

  • CP2K supports flexible combinations of empirical force fields, DFT-based techniques, and quantum-chemical methods, including linear combinations of potentials and propagation of the lowest potential-energy surfaces.More sophisticated embedding methods are available alongside these simpler combinations.Source: §XI · pp. 32–33

A. QM/MM methods

  • CP2K’s QM/MM implementation uses an additive energy partition into quantum, molecular-mechanics, and QM/MM interaction terms, with Quickstep providing EQM and FIST providing EMM.The interaction term contains all non-bonded contributions between quantum and classical subsystems.Source: A. QM/MM methods · p. 33

1. QM/MM for isolated systems

  • For isolated systems, Gaussian MM charge distributions are assigned across grid levels, and the resulting multi-grid potential is interpolated from the coarsest to finest grid.The coarsest grid includes every atom’s smooth long-range GEEP component.Source: §1 · p. 33
  • The GEEP multi-grid decomposition lowers the electrostatic-potential prefactor from N_f*N_f*N_f to N_c*N_c*N_c, with other operations negligible for systems exceeding roughly 600–800 MM atoms.GEEP expands the electrostatic potential into Gaussian terms with different cutoffs and maps them onto separate real-space grid levels.Source: §1 · p. 33, A. QM/MM methods · p. 33
  • A four-level multi-grid implementation achieves an approximately 512-fold speed-up, making QM/MM electrostatic-potential evaluation two orders of magnitude faster than direct analytical grid evaluation.The implementation uses commensurate grids with N_f/N_c = 2^3(N_grid−1), and four grid levels.Source: §1 · p. 33

2. QM/MM for periodic systems

  • Periodic MM replicas affect only the long-range residual term, which is evaluated efficiently in reciprocal space using an Ewald-like treatment and mapped onto the coarsest grid.The periodic MM potential then uses the same interpolation and restriction operators as isolated systems.Source: §2 · p. 34
  • Periodic QM/MM calculations also decouple periodic QM images to restore correct QM periodicity, especially when the QM subsystem has net charge or a significant dipole moment.CP2K uses a technique proposed by Blöchl, with a comprehensive description reported in.Source: §2 · p. 34

3. Image charge augmented QM/MM … 2. Implementation

CP2K extends QM/MM and charge-fitting workflows for interfaces and periodic systems, while density functional embedding combines DFT environments with higher-level correlated treatments of embedded clusters. These implementations enforce physically consistent electrostatics or densities and support a broad range of system types and electronic-structure methods.

3. Image charge augmented QM/MM

  • The IC-QM/MM scheme treats adsorbates with KS-DFT and metallic substrates with MM, explicitly screening electrostatic interactions through image charges.Lennard-Jones-type potentials model dispersion and Pauli repulsion, while the image-charge formulation accounts for metallic polarization.Source: §3 · p. 34
  • Self-consistent Gaussian image charges impose a constant potential inside the metal, ensuring Ve(r) + Vm(r) = V0.The image-charge coefficients are determined from the constant-potential condition, correcting the electrostatics illustrated in Fig. 9.Source: §3 · p. 34, Figure 9 · p. 34
  • The IC-QM/MM implementation requires only metal-atom positions, has negligible current overhead, strengthens molecule–metal interactions, and partially captures molecular electronic rearrangement.It is intended for systems where molecule–molecule interactions dominate and the surface primarily acts as a template, including water/Pt(111), organic films on Au(100), and DNA at gold interfaces.Source: §3 · pp. 34–35

4. Partial atomic charges from electrostatic potential fitting

  • RESP charges fit quantum-mechanical electrostatic potentials by least squares, with grid-point selection tailored to interaction-relevant regions and slab surfaces.For slab-like systems, CP2K supports sampling that reproduces the potential above corrugated surfaces.Source: §4 · p. 35
  • For periodic systems, CP2K represents RESP charges as Gaussian functions on a real-space grid and uses GPW to compute periodic RESP potentials.The REPEAT implementation fits potential variance rather than absolute differences, stabilizing periodic fits and making it particularly suitable for non-slab systems.Source: §4 · p. 35
  • Periodic RESP and REPEAT charges support QM/MM simulations of surface systems and applications including hBN on Rh(111), supramolecular structures, metal-organic frameworks, and graphane.The hBN@Rh(111) charges were used for liquid-water-film QM/MM simulations.Source: §4 · p. 35

1. Theory

  • DFET embeds a correlated wavefunction treatment of a chemically relevant cluster in a DFT environment through a local potential vemb(r) that reconstructs the total DFT density.The embedding potential is optimized using the Wu–Yang functional until the embedded subsystem densities sum to the total-system density.Source: §1 · pp. 35–36

2. Implementation

  • The DFET workflow iterates subsystem DFT calculations to obtain vemb before applying the higher-level cluster calculation, with ASPC-based wavefunction extrapolation reducing computational cost.It supports open- and closed-shell variants, GPW pseudopotential calculations, hybrid DFT, MP2, RPA, embedded-cluster properties, differing subsystem basis sets, and real-space-grid or Gaussian representations of vemb.Source: §2 · p. 36

C. Implicit solvent techniques · D. Poisson solvers

CP2K provides implicit-solvent models that represent bulk solvent electrostatically while addressing interface smoothness and computational limitations. Its generalized Poisson solver extends supported boundary conditions, dielectric functions, convergence behavior, and force consistency for complex atomistic simulations.

C. Implicit solvent techniques

  • Implicit solvent methods reduce explicit-solvent cost by treating bulk solvent as a continuous dielectric medium while providing statistically averaged solvent effects.The first solvent layers can affect chemistry directly, whereas bulk solvent is described mainly through electrostatic interactions; SCRF in a sphere is the simplest implemented method.Source: C. Implicit solvent techniques · p. 36
  • COSMO and PCM represent solute-shaped cavities, but dielectric discontinuities can produce non-continuous forces that affect structural relaxation and AIMD energy conservation.Smoothed self-consistent dielectric functions based on the electronic density ρelec(r) were proposed to overcome these problems.Source: C. Implicit solvent techniques · p. 36
  • Density-dependent dielectric models require self-consistent polarization-charge iterations within each SCF step, creating potentially non-negligible computational overhead.The thermal-motion term Gtm and volume-change term P∆V are often ignored and are not yet available in CP2K.Source: C. Implicit solvent techniques · pp. 36–37

D. Poisson solvers

  • The generalized Poisson solver targets complex boundary configurations beyond existing methods, including nanoscale devices and other ab-initio applications.It was developed because conventional approaches, despite broad success, do not cover increasingly complex boundary setups.Source: D. Poisson solvers · p. 37
  • The generalized solver converges exponentially with respect to the density cutoff and combines advantages associated with plane-wave and real-space Poisson methods.Its formulation imposes the stated boundary setups through an equivalent constrained variational problem.Source: D. Poisson solvers · p. 37
  • The solver supports periodic or homogeneous von Neumann boundaries, arbitrarily shaped Dirichlet-potential regions, and sufficiently smooth position-dependent dielectric functions.These features model insulating interfaces and nanoscale-device contacts while allowing flexible dielectric descriptions.Source: D. Poisson solvers · p. 37
  • Consistent ionic forces make the generalized solver applicable to energy-conserving BOMD and EMD simulations.The implementation therefore extends its use to molecular-dynamics applications requiring consistent forces.Source: D. Poisson solvers · p. 38

XII. DBCSR LIBRARY … XIV. TECHNICAL AND COMMUNITY ASPECTS

The paper presents CP2K’s scalable infrastructure, emphasizing DBCSR’s distributed sparse-matrix operations, accelerator-aware execution, and interfaces for transport, plane-wave DFT, and external programs. It also reports quantitative validation for batched kernels, NEGF simulations, and CP2K+SIRIUS energies against QE.

XII. DBCSR LIBRARY

  • DBCSR provides distributed block-sparse and dense matrix operations across occupancies from 0.01% to dense, using blocked CSR storage on two-dimensional MPI grids.Its operations include matrix multiplication, sums, dot products, transposes, and traces, with block sizes determined by molecular basis-set properties.Source: §XII · p. 38

A. Message passing interface parallelization

  • General matrices use Cannon communication scaling as O(1/Tall-and-skinny matrices instead use an optimized algorithm with O(1) communicated data per process, while asynchronous MPI enables communication-computation overlap.Source: A. Message passing interface… · p. 38

B. Local multiplication

  • Local multiplication traverses small dense blocks cache-obliviously, batches at most 30,000 products, and executes them concurrently on CPUs and GPUs with minimized transfers.GPU execution uses coalesced and asynchronous operations, double buffering, and round-robin multi-GPU distribution across MPI ranks.Source: B. Local multiplication · pp. 38–39

C. Batched execution

  • 2–4x speedup over batched DGEMM in cuBLAS is achieved for {m, n, k} < 32, while larger sizes use cuBLAS or BLAS after performance saturation.LIBSMM ACC combines JIT compilation, autotuning, and machine-learning prediction across available kernel shapes.Source: C. Batched execution · p. 39

XIII. INTERFACES TO OTHER PROGRAMS

  • CP2K supports external-program workflows through library and CP2K-shell interfaces, connections to i-PI, PLUMED, and PyRETIS, and trajectory analysis with several tools.The library can also be embedded in external executables.Source: §XIII · p. 40

A. Non-equilibrium Green’s function formalism

  • Devices with realistic sizes and structural configurations can be routinely simulated by integrating CP2K’s DFT Hamiltonians with the OMEN NEGF solver.OMEN modifies the Hamiltonian for open boundaries and calculates charge and current densities using NEGF or wave-function/quantum-transmitting-boundary formalisms.Source: A. Non-equilibrium Green’s function… · p. 40

B. SIRIUS: Plane wave density functional theory support

  • 10−6 Ha agreement is obtained between CP2K + SIRIUS and Quantum ESPRESSO potential energies after removing a constant offset during Si7Ge AIMD.The total energy remains conserved up to 10−6 Ha, while the two codes show the same potential-energy evolution.Source: B. SIRIUS: Plane wave density… · p. 40, Figure 12 · p. 41

XIV. TECHNICAL AND COMMUNITY ASPECTS; D. Outlook

  • CP2K combines MPI, OpenMP, accelerator support, autogenerated and autotuned code, and external libraries to provide scalable performance and library-like extensibility.Its outlook identifies low arithmetic intensity and memory bandwidth as motivations for lower or mixed precision.Source: D. Outlook · pp. 39–40, §XIV · pp. 40–41

A. Hardware Acceleration … F. GitHub and general tooling

CP2K combines accelerator support, accuracy-preserving approximate computing, scalable and accurate electronic-structure methods, workflow automation, and open development tooling. Its implementations target GPUs, FPGAs, reproducible workflows, and broad software quality assurance.

A. Hardware Acceleration

  • CP2K supports mature GPU acceleration and emerging FPGA support through DBCSR offloading, GPU FFT libraries, and a dedicated FPGA FFT interface.The released FPGA library targets complex-valued single-precision 3D FFTs of sizes 32^3 to 128^3 on Intel Arria 10 and Stratix 10 devices.Source: A. Hardware Acceleration · p. 41

B. Approximate Computing

  • Approximate computing preserves exact ensemble-averaged expectations by rigorously compensating low-precision force errors within molecular-dynamics Langevin dynamics.The method treats inaccuracies as unbiased additive white noise and adaptively determines the damping coefficient to satisfy the fluctuation-dissipation condition [184, 396, 397].Source: B. Approximate Computing · pp. 41–42

C. Benchmarking

  • CP2K demonstrates strong scaling to thousands of CPU cores for GAPW Hartree-Fock, periodic RI-MP2, and linear-scaling DFT calculations.Figure 13 covers systems ranging from 216 lithium-hydride atoms to 2048 condensed-phase water molecules and up to 10240 CPU cores.Source: C. Benchmarking · p. 42, Figure 13 · p. 42
  • Up to one order of magnitude efficiency improvement is observed for million-atom STMV simulations using approximate computing, while force noise can be compensated by a modified Langevin equation.The approach applies to matrix-square-root and matrix-sign operations in periodic GFN2-xTB calculations, with energy and density-matrix deviations monitored.Source: C. Benchmarking · p. 42

D. MOLOPT basis set and deltatest

  • CP2K’s MOLOPT basis sets target diverse chemical environments and generally perform well in Δ-tests against WIEN2k and Abinit references.The comparison uses identical GTH pseudopotentials for CP2K and Abinit, with residual deviations attributed to the pseudization approach.Source: D. MOLOPT basis set and deltatest · pp. 42–43, Figure 15 · p. 43

E. CP2K workflows

  • CP2K interfaces with ASE and AiiDA to automate, reproduce, and share computational workflows.ASE supports rapid prototyping and common atomistic tasks, whereas AiiDA provides robust event-based high-throughput workflows with database-recorded provenance [402, 403, 404].Source: E. CP2K workflows · p. 43
  • AiiDA traces data dependencies and records provenance while handling rare failures across high-throughput workflows involving thousands of computations.Its workflows use functional building blocks with defined inputs and outputs and automatically screen candidate-material databases for target properties.Source: E. CP2K workflows · p. 43

F. GitHub and general tooling

  • GitHub-based development combines pull-request review, mandatory regression testing, broad supercomputer coverage, automated code analysis, and packaged installation support.Regression tests cover over 80% of CP2K code across MPI, OpenMP, CUDA/HIP, and FPGA variants, while tooling supports Fortran style, dependencies, and Linux, Spack, and EasyBuild installations.Source: F. GitHub and general tooling · pp. 43–44
Loading 2003.03868v2…