Implementation and Validation of Constrained Density Functional Theory Forces in the CP2K Package

Constrained density functional theory (CDFT) is a powerful tool for the prediction of electron transfer parameters in condensed phase simulations at a reasonable computational cost. In this work we present an extension to CDFT in the popular mixed Gaussian/plane wave electronic structure package CP2K, implementing the additional force terms arising from a constraint based on Hirshfeld charge partitioning. This improves upon the existing Becke partitioning scheme, which is prone to give unphysical atomic charges. We verify this implementation for a variety of systems: electron transfer in (H2O)2+ in a vacuum, electron tunnelling between oxygen vacancy centers in solid MgO, and electron self-exchange in aqueous Ru2+–Ru3+. We find good agreement with previous plane-wave CDFT results for the same systems, but at a significantly lower computational cost, and we discuss the general reliability of condensed phase CDFT calculations.

: Verification of analytical forces against forces calculated from centred finite differences of the total energy for a helium dimer, both with (right) and without (left) periodic boundary conditions. For the system under periodic boundary conditions, the helium dimer interacts at the same distance except with their periodic image and therefore the resultant force is the same. An isosurface of the weight function is shown on the bottom left of each figure.

CDFT geometry optimisation: MgO
For completeness, we provide all the electronic couplings and reorganisation energies tabulated below. Values calculated using Becke charge partitioning are also included.

CDFT geometry optimisation: 2D Pyrene COF
The usual process for calculating the reorganisation energy of organic molecules would be in the gas phase, however CDFT provides an ability to constrain the excess charge to a single unit in a periodic crystal and therefore account for the full outer-sphere reorganisation energy.
We show an example of hole transfer in a 2D pyrene based covalent organic framework (COF), where the subsequent geometry optimisation diverges and leads to an increase in energy of 54 eV. Several different constraint definitions have been tested: the charge difference between two units shown below, the charge difference between one unit and the rest of the system, and an absolute charge constraint over one unit. These, in addition to constraints defined including or excluding the acetylene linkers, all lead to a large IASD and divergent geometry optimisation. This can be attributed to the effect of a too strong constraint, that the polaron in these materials is band-like 3 and attempting to constrain the excess charge to a single unit is an inappropriate choice of constraint as the underlying functional is not able to correctly describe the resulting electronic state. Figure 3: CDFT weight function defined as the charge difference between two units in a 3x3 2D pyrene cof (left), increase in energy per atom as a function of geometry optimisation step (middle) and increase in Integrated Absolute Spin Density as a function of geometry optimisation step (right). The large initial IASD of 1.46 and subsequent increase is an indication of fractional charge transfer, which generally occurs when the DFT functional is unable to appropriately describe the charge localised state.

CDFT geometry optimisation: Pentacene crystal
Similar to the pyrene-cof, it would be useful to calculate the reorganisation energy for a pentacene crystal accounting for the full outer-sphere reorganisation energy. Different constraints have been tested, constraining either the absolute charge or the charge difference between two pentacene molecules. In all cases, the IASD is large and the geometry optimisation diverges with an increase in energy of 30 eV after 13 geometry optimisation steps. Figure 4: CDFT weight function defined as the charge difference between two pentacene molecules in a 3x2x1 pentacene cyrstal (left), increase in energy per atom as a function of geometry optimisation step (middle) and increase in Integrated Absolute Spin Density as a function of geometry optimisation step (right). The large initial IASD of 1.55 and subsequent increase is an indication of fractional charge transfer, which generally occurs when the DFT functional is unable to appropriately describe the charge localised state. Note that the repetition of the CDFT weight function (left) above and below the pentacene molecules is a visualisation artefact due to the non-cubic unit cell.

CDFT geometry optimisation: Pentacene in vacuum
Given the surprising results where the condensed phase CDFT geometry optimisation of pentacene fails, it is useful to confirm that CDFT works well for the pentacene dimer in vacuum. Here we confirm the exponential decay of the electronic coupling for both holes and electrons in a π-stacked pentacene dimer, despite their large IASD. This is consistent with the results published for the HAB11 dataset. 4-6  We further check the CDFT geometry optimisation of an excess hole for the π-stacked pentacene dimer in vacuum. For this example while the IASD increases, the total energy of the system decreases and the CDFT geometry optimisation succeeds with a reorganisation energy of 0.44 eV. This highlights a particular sensitivity of CDFT to condensed phase calculations, and suggests that the IASD is not always a reliable indicator for the breakdown of CDFT. It is possible however that CDFT only works for this system by fortuitous cancellation of errors, as described by the work of Van Voorhis and co-workers. 7 Figure 5: Decrease in total energy as a function of geometry optimisation step (left) and increase in Integrated Absolute Spin Density as a function of geometry optimisation step (right). Despite the large initial IASD of 1.34 and subsequent increase, the total energy decreases and the geometry optimisation converges.

CDFT-MD: H 2 O + -H 2 O in vacuum
Alternative Table 2 from main text, with added column providing bond lengths and angles from CDFT-MD simulation of water dimer in vacuum. While not statistically converged, these values agree well with the geometry optimised results.

Implementation of Hirshfeld CDFT
An implementation of CDFT forces using Hirshfeld partitioning of the electron density is now available in CP2K version 10.
Below, we present CDFT-MD energy conservation including data for a constraint convergence of 1 × 10 −2 e for condensed phase systems, and the time taken for an average CDFT-MD step relative to an equivalent DFT-MD step as a function of the constraint con- vergence. An important observation is that while hybrid DFT is more expensive than GGA DFT, the additional cost of CDFT is lower at the hybrid DFT level. This is likely a result of the under-binding of excess charge at the GGA level, which makes convergence of localised charges with CDFT more challenging and therefore the number of SCF steps is increased.