The Perdew-Zunger Self-Interaction Correction (PZ-SIC)

The self-interaction corrected density functional with PZ-SIC [1] has the following form:

\[E^{PZ-SIC}[\{n_i\}] = E^{DFA}[n] - \beta \sum_{i=1}^{N} \left(E_H[n_i] + E^{DFA}_{xc}[n_i]\right)\]

here \(E^{DFA}\) - density functional approximation, \(E_H\) - Hartree energy and \(E^{DFA}_{xc}\) - exchange correlation part of density functional approximation. \(n\) is the total density and \(n_i = |\psi_i (\mathbf{r})|^{2}\) - orbital density. \(\beta\) is a scaling factor.

The SIC functional is not a unitary invariant functional and is dependent on orbital densities. Therefore, the fully variational approach [2], [3] is used to find the optimal orbitals which provide the ground state energy.

Example

The implementation support all three modes (PW, FD, and LCAO) using with direct minimization described here . Since the functional is not a unitary invariant functional, it is necessary to employ complex orbitals to find the lowest energy state. Here is an example using FD mode:

import numpy as np
from ase import Atoms
from gpaw import FD, GPAW
from gpaw.directmin.etdm_fdpw import FDPWETDM

# Water molecule:
d = 0.9575
t = np.pi / 180 * 104.51
H2O = Atoms('OH2',
            positions=[(0, 0, 0),
                       (d, 0, 0),
                       (d * np.cos(t), d * np.sin(t), 0)])
H2O.center(vacuum=5.0)

calc = GPAW(mode=FD(force_complex_dtype=True),
            xc='PBE',
            occupations={'name': 'fixed-uniform'},
            eigensolver=FDPWETDM(localizationtype='PM_PZ',
                                 functional={'name': 'PZ-SIC',
                                             'scaling_factor':
                                                 (0.5, 0.5)},
                                 grad_tol_pz_localization=1.0e-4),
            mixer={'backend': 'no-mixing'},
            symmetry='off'
            )

H2O.set_calculator(calc)
H2O.get_potential_energy()
H2O.get_forces()

To use PW mode, just import PW mode and replace FD with PW. While here is the example for LCAO mode:

import numpy as np
from ase import Atoms
from gpaw import GPAW, LCAO
from gpaw.directmin.etdm_lcao import LCAOETDM

# Water molecule:
d = 0.9575
t = np.pi / 180 * 104.51
H2O = Atoms('OH2',
            positions=[(0, 0, 0),
                       (d, 0, 0),
                       (d * np.cos(t), d * np.sin(t), 0)])
H2O.center(vacuum=5.0)

calc = GPAW(mode=LCAO(force_complex_dtype=True),
            xc='PBE',
            occupations={'name': 'fixed-uniform'},
            eigensolver=LCAOETDM(localizationtype='PM_PZ',
                                 functional={'name': 'PZ-SIC',
                                             'scaling_factor':
                                                 (0.5, 0.5)}),
            mixer={'backend': 'no-mixing'},
            nbands='nao',
            symmetry='off'
            )

H2O.calc = calc
H2O.get_potential_energy()
H2O.get_forces()

If you use this module, please refer to implementation papers Refs. [2], [3].

References