Constrained DFT (cDFT)

Introduction

cDFT is a method for build charge/spin localized or diabatic states with a user-defined charge and/or spin state. As such cDFT, is a useful tool for widening the scope of ground-state DFT to excitation processes, correcting for self- interaction energy in current DFT functionals, excitation energy, and electron transfer as well as parametrizing model Hamiltonians, for example.

Theoretical Background

The general cDFT methodology is reviewed in [1] and the publication of the GPAW implementation is available in [2]. In short, the cDFT works by specifying an additional constraining term to the KS functional. The role of this constraint is to enforce a user specified charge/spin state (\(N_c\)) on the chosen regions of atoms. The constrained regions are specified by a weight function \(w(\mathbf{r})\) and the strength of the constraining potential acting on the specified region is \(V_c\). With these definitions the new energy functional with the constraint is

\[F[n(\mathbf{r}), V_c] = E^{KS}[n] + \sum_c V_c \sum_{\sigma}\left(\int d\mathbf{r}w_c^{\sigma}(\mathbf{r}) n^{\sigma}(\mathbf{r})-N_c\right)\]

where \(E^{KS}\) is the Kohn-Sham energy, \(c\) specifies the region, and \(\sigma\) is the spin variable. It is also seen that \(V_c\) is also a Lagrange multiplier. The modified energy functional leads to an additional external potential

\[v_{\rm eff}^{\sigma}=\dfrac{\delta E^{KS}[n(\mathbf{r})]} {\delta n(\mathbf{r})} + \sum_c V_cw_c^{\sigma}(\mathbf{r})\]

This is just a sum of the usual KS potential and the constraining potential which is also used in the self-consistent calculation. The constraint is further enforced by introducing the convergence criteria

\[C \geq \bigg\lvert \sum_{\sigma} \int {\rm d}\mathbf{r} w_c^{\sigma}(\mathbf{r})n^{\sigma}(\mathbf{r}) - N_c \bigg\rvert \, ,\forall\, c\]

The \(V_c\) is self-consistently optimized so that the specified constraints are satisfied. Formally, this is written as

\[F\left[V_{c}\right]= \min _{n} \max _{\left\{V_{c}\right\}} \left[E^{\mathrm{KS}}[n(\mathbf{r})]+ \sum_{c} V_{c}\left(\int \mathrm{d} \mathbf{r} w_{c}(\mathbf{r}) n(\mathbf{r})-N_{c}\right)\right]\]

\(V_c\) is obtained from

\[\frac{\mathrm{d} F}{\mathrm{d} V_{c}} = \int \mathrm{d} \mathbf{r} w_{c}(\mathbf{r}) n(\mathbf{r})-N_{c}=0\]

In the end, one ends up with a modified electron/spin density localized on the chosen atoms.

Notes on using cDFT

The weight function

In the GPAW implementation a Hirschfeld partition scheme with atom-centered Gaussian functions is used. These Gaussian have two tunable parameters: the cut-off \(R_c\) and the width \(\mu\). If the constrained region cuts a covalent bond, the results are sensitive to width parameter. There is no universal guide for choosing the width in such cases. A sensible choice is to compute match the Gaussian-Hirschfeld charges with e.g. Bader charges. The function get_number_of_electrons_on_atoms() helps in this process.

class gpaw.cdft.cdft.CDFT(calc, atoms, charge_regions=[], charges=None, spin_regions=[], spins=None, charge_coefs=None, spin_coefs=None, promolecular_constraint=False, txt='-', minimizer_options={'gtol': 0.01}, Rc={}, mu={'F': 0.7, 'Li': 0.5, 'O': 0.7, 'V': 0.5}, method='BFGS', forces='analytical', use_charge_difference=False, compute_forces=True, maxstep=100, tol=0.001, bounds=None, restart=False, hess=None)

Constrained DFT calculator.

calc: GPAW instance

DFT calculator object to be constrained.

charge_regions: list of list of int

Atom indices of atoms in the different charge_regions.

spin_regions: list of list of int

Atom indices of atoms in the different spin_regions.

charges: list of float

constrained charges in the different charge_regions.

spins: list of float

constrained spins in the different charge_regions. Value of 1 sets net magnetisation of one up/alpha electron

charge_coefs: list of float

Initial values for charge constraint coefficients (eV).

spin_coefs: list of float

Initial values for spin constraint coefficients (eV).

promolecular_constraint: bool

Define charge and/or spin constraints from promolecular densities, see: dx.doi.org/10.1021/cr200148b Eq. 29-31 If true, user specified charges/spins are overwritten! The atoms (of Atoms object) specifying the charge/spin regions need to contain have correct charge/spin state! (atoms.set_initial_charges([atomic_charges]) and atoms.set_initial_magnetic_moments([moments]))

txt: None or str or file descriptor

Log file. Default id ‘-’ meaning standard out. Use None for no output.

minimizer_options: dict

options for scipy optimizers, see:scipy.optimize.minimize

method: str

One of scipy optimizers, e.g., BFGS, CG

forces: str

cDFT weight function contribution to forces ‘fd’ for finite difference or ‘analytical’

difference: bool

If True, constrain charge difference between two regions Then charge_regions needs two regions and charges needs only one item which is the charge difference between the two regions, the first beign donor, the second acceptor

If False, each region is treated with the corresponding charge constraint

compute_forces: bool

Should the forces be computed?

restart: bool

starting from an old calculation from gpw

hess: ‘2-point’, ‘3-point’, ‘cs’

scipy hessian approximation

calculate(atoms, properties, system_changes)

Do the calculation.

properties: list of str

List of what needs to be calculated. Can be any combination of ‘energy’, ‘forces’, ‘stress’, ‘dipole’, ‘charges’, ‘magmom’ and ‘magmoms’.

system_changes: list of str

List of what has changed since last calculation. Can be any combination of these six: ‘positions’, ‘numbers’, ‘cell’, ‘pbc’, ‘initial_charges’ and ‘initial_magmoms’.

Subclasses need to implement this, but can ignore properties and system_changes if they want. Calculated properties should be inserted into results dictionary like shown in this dummy example:

self.results = {'energy': 0.0,
                'forces': np.zeros((len(atoms), 3)),
                'stress': np.zeros(6),
                'dipole': np.zeros(3),
                'charges': np.zeros(len(atoms)),
                'magmom': 0.0,
                'magmoms': np.zeros(len(atoms))}

The subclass implementation should first call this implementation to set the atoms attribute and create any missing directories.

get_number_of_electrons_on_atoms()

Return the number of electrons with each gaussian.

implemented_properties: List[str] = ['energy', 'forces']

Properties calculator can handle (energy, forces, …)

Optimizing \(V_c\)

Updating and optimizing the Lagrange multipliers \(V_c\) is achieved using Scipy optimization routines. As it is easy to compute the gradient of the energy wrt \(V_c\), gradient based optimizers are recommended. The best performing optimizer seems to be the L-BFGS-B. The accuracy of the optimization is controlled mainly by the minimizer_options={'gtol':0.01}) parameter which measurest the error between the set and computed charge/spin value. A typical gtol value is around 0.01-0.1.

Choosing the constraint values and regions

Both charge and spin constraints can be specified. The charge constraints are specified to neutral atoms: specifying charges = [1] constrains the first regions to have a net charge of +1. Note, that if the usual DFT calculation yields e.g. a Fe ion with charge +2.5, specifying the charge constraint +1 will result in +1, not an additional hole on Fe with a charge state +3.5!

A constrained regions may contain several atoms and also several constrained regions can be specified. However, converging more than two regions is usually very difficult.

Tips for converging cDFT calculations

Unfortunaly, the cDFT sometimes exhibits poor convergence.

1. Choose a meaningful constraints and regions

2. Try to provide a good initial guess for the \(V_c\). In the actual calculation this initial guess given by the parameter charge_coefs or spin_coefs.

3. Use L-BFGS-B to set bounds for \(V_c\) by using minimizer_options={'bounds':[min,max]}).

4. Converge the underlying DFT calculation well i.e. use tight convergence criteria.

Constructing diabatic Hamiltonians

One of the main uses for cDFT is constructing diabatic Hamiltonians which utilize cDFT states as the diabats. For instance, a 2x2 Hamiltonian matrix would have diagonal elements \(H_{11}\) and \(H_{22}\) as well as off-diagonals \(H_{12}\) and \(H_{21}\). \(H_{11}\) and \(H_{22}\) values are given directly by cDFT energies (not cDFT free energies). The off-diagonal coupling elements \(H_{12}\) and \(H_{21}\) are also needed and often utilized in e.g. Configuration Interaction-cDFT, in parametrizing model Hamiltonians or in computing non-adiabatic charge transfer rates. Note, that all parameters need to be computed at the same geometry.

The coupling elements are computed using the CouplingParameters class. There are several options and optional inputs. There are two main methods for computing coupling elements: the original cDFT approach from [1] and a more general overlap or Migliore method detailed in [3].

cDFT coupling constant

This method is the original approach in the context of cDFT. It works well for simple system. For complex system it becomes quite sensitive to small errors in the Lagrange multipliers and constraints, and the overlap method is recommended instead. The inputs for the calculation are two cDFT objects with different constraints. The coupling constant is computed using get_coupling_term.

Overlap coupling constant

This approach has been found to perform very well for a range of systems. However, it has one major limitation: it can only be used if the two diabatic states have different energies. In addition to the two cDFT states/objects also a normal DFT calculator without any constraints is needed. The coupling constant is computed using get_migliore_coupling.

Additional comments

In [2] the coupling constants were computed by computing all-electron wave functions on a grid. However, this is quite slow and much faster implementation utilizing only pseudo wave functions and atomic corrections has been added. For the tested cases both give equally good values for the coupling constant. Hence, it is recommended to use the pseudo wave functions which is set by AE=False.

The quantities needed for computing the coupling constants can be parallellized only over the grid.

Example of hole transfer in \(He_2^+\)

Both the cDFT calculation and extraction of the coupling element calculation are demonstrated for the simple \(He_2^+\) system.

from ase import Atoms
from gpaw import GPAW, FermiDirac, Davidson, Mixer
from gpaw.cdft.cdft import CDFT
from gpaw.cdft.cdft_coupling import CouplingParameters

# Set up the system
distance = 2.5
sys = Atoms('He2', positions=([0., 0., 0.], [0., 0., distance]))
sys.center(3)
sys.set_pbc(False)
sys.set_initial_magnetic_moments([0.5, 0.5])

# Calculator for the initial state
calc_a = GPAW(
    h=0.2,
    mode='fd',
    basis='dzp',
    charge=1,
    xc='PBE',
    symmetry='off',
    occupations=FermiDirac(0., fixmagmom=True),
    eigensolver=Davidson(3),
    spinpol=True,  # only spin-polarized calculations are supported
    nbands=4,
    mixer=Mixer(beta=0.25, nmaxold=3, weight=100.0),
    txt=f'He2+_initial_{distance:3.2f}.txt',
    convergence={
        'eigenstates': 1.0e-4,
        'density': 1.0e-1,
        'energy': 1e-1,
        'bands': 4})

# Set initial state cdft
cdft_a = CDFT(
    calc=calc_a,
    atoms=sys,
    charge_regions=[[0]],  # choose atom 0 as the constrained region
    charges=[1],  # constrain +1 charge
    charge_coefs=[2.7],  # initial guess for Vc
    method='L-BFGS-B',  # Vc optimization method
    txt=f'He2+_initial_{distance:3.2f}.cdft',  # cDFT output file
    minimizer_options={'gtol': 0.01})  # tolerance for cdft

# Get cdft energy
sys.calc = cdft_a
sys.get_potential_energy()

# the same for the final state
calc_b = GPAW(h=0.2,
              mode='fd',
              basis='dzp',
              charge=1,
              xc='PBE',
              symmetry='off',
              occupations=FermiDirac(0., fixmagmom=True),
              eigensolver=Davidson(3),
              spinpol=True,
              nbands=4,
              mixer=Mixer(beta=0.25, nmaxold=3, weight=100.0),
              txt=f'He2+_final_{distance:3.2f}.txt',
              convergence={
                  'eigenstates': 1.0e-4,
                  'density': 1.0e-1,
                  'energy': 1e-1,
                  'bands': 4})

cdft_b = CDFT(
    calc=calc_b,
    atoms=sys,
    charge_regions=[[1]],  # choose atom 1
    charges=[1],  # constrained charge +1
    charge_coefs=[2.7],
    method='L-BFGS-B',
    txt=f'He2+_final_{distance:3.2f}.cdft',
    minimizer_options={'gtol': 0.01})

sys.calc = cdft_b
sys.get_potential_energy()

# Now for the coupling parameter
coupling = CouplingParameters(cdft_a, cdft_b, AE=False)  # use pseudo orbitals
H12 = coupling.get_coupling_term()  # use original cDFT method

The most important cDFT results are found in the .cdft files. For instance, the errors, iterations, and cDFT parameters are shown. Also, the energy can be found in this file. The most relevant energy is the Final cDFT Energy (not the free energy).

References

[1] (1,2)

B. Kaduk, T. Kowalczyk, T. Van Voorhis, Constrained Density Functional Theory, Chem. Rev., 112 321–370 (2012)

[2] (1,2)

M. Melander, E. Jońsson, J.J. Mortensen, T. Vegge,J.M. Garcia-Lastra, Implementation of Constrained DFT for Computing Charge Transfer Rates within the Projector Augmented Wave Method, J. Chem. Theory Comput., 12, 5367−5378 (2016)

[3]

A. Migliore, Nonorthogonality Problem and Effective Electronic Coupling Calculation: Application to Charge Transfer in π-Stacks Relevant to Biochemistry and Molecular Electronics, J. Chem. Theory Comput., 7, 1712-1725 (2011)