Self-Interaction Correction

In the ADF2004.01 version the Perdew-Zunger (PZ) self-interaction energy correction (SIC) with the Krieger-Li-Iafrate (KLI) approximation to the self-interaction corrected optimized effective potential (OEP) is implemented[47-49]. The block key SICOEP should be used.

Note: The existing auxiliary fitting sets, employed in ADF, were optimized for the calculation of the Coulomb potential of the total electron density. These standard fitting sets are quite flexible in the valence region, but do not include functions of high angular momentum in the inner regions. It was found (Ref [47]) that self-consistent SIC calculations with the standard fitting sets result in very poor SCF convergence, and large fit incompleteness corrections, particularly for the orbitals with substantial p and d contributions. Once the auxiliary fit sets are augmented with additional functions of high angular momentum, and reoptimized, the SCF convergence problems largely disappear. The reason is that the fit set must be able to approximate densities of each individual localized orbital, both in direction, and shell structure.

The key

SINGULARFIT FAST

activates fast(er) treatment of linearly dependent fits. It precomputes and stores singular value decomposition of fit overlap integrals, making fitting itself faster. Because SIC needs to fit a lot of densities on each iteration, it can improve performance by an order of magnitude. Use of 'SINGULARFIT FRUGAL' requests the old treatment of linearly-dependent fits, recomputing SVD parameters on each iteration. This is always slower than 'SINGULARFIT FAST', but uses less disk space.

The localization transformation is essential for obtaining realistic total and atomization energies. For example, without localization, SIC-Vosko-Wilk-Nusair (VWN) predicts molecular oxygen to be unstable with respect to the dissociation to two neutral oxygen atoms. At the same time, when using the localized orbitals, SIC-VWN compares favorably to both VWN itself, and to sophisticated gradient-corrected functionals. In ADF the Foster-Boys localization procedure is used.

Remarks: Fractional occupation are treated right. GGA functionals are supported, even for open-shell. SIC works for Linear Transit calculations. Frozen core SIC potentials are computed (GGA). The code is not well suitable for treating d and f electrons within the valence shell. Ds and Fs are fine within the frozen core. The SIC code does not allow parallel calculations (yet). Geometry optimizations are not yet possible. SIC in combination with spin-orbit coupling is not possible yet.

Usage of the block key SICOEP:

SICOEP
 {IPRINT n}
 {NOLOCALIZE}
 {LOCALIZE {thrs}}
 {LDA name {Xa}}
 {GGA name}
 {SELFCONSISTENT n}
 {POSTSCF}
 {SKIPCYCLES every {start}}
 {STABLE frac}
 {NHOMO Nalpha {Nbeta}}
 {CORE cor}
 {DENSITY mode}
 {SHIPV filename}
 {READV filename}
 {READMOS filename {FREEZE}}
 {WRITEMOS filename}
 {NAILCANONICALS {eps}}
 {NPFITS n}
end

IPRINT n

-2 = no printing except for fatal errors
-1 = condensed printing
0 = normal printing
1 = verbose printing
5 = debug printing
10 = exorcism printing

NOLOCALIZE

Calculate SIC energies for canonical MOs. Except for atoms and certain small molecules, this is guaranteed to produce non-optimal SIC energies, and/or lead to severe convergence problems.

LOCALIZE {thrs}

Calculate SIC energies for localized MOs. In this version, Boys-Foster procedure is used to obtain localized orbitals. By default, all occupied orbitals will participate in the Boys-Foster procedure. This can be changed by supplying thrs parameter, which will only allow mixing of canonical orbitals within thrs eV of each other. Calculate SIC energies for canonical MOs. Except for atoms and certain small molecules, this is guaranteed.

LDA name {Xa}

Requests a specific local density functional in Perdew-Zunger energy correction and KLI contribution to the XC potential (If name is XALPHA, the default alpha of 0.7 can be changed by specifying the additional argument). Normally, SIC code will use LDA functional requested (explicitly or implicitly) in the XC keyblock. Specifying 'NONE' with suppress LDA contribution in SIC energy and potential. Because the Perdew-Zunger energy correction, and the corresponding term in the XC potential are trying to remove a non-physical self-interaction contribution from the Kohn-Sham E(XC)/v(XC) contribution, this key should only be used for debugging.

GGA name

Requests a specific gradient-corrected functional in the SIC part, overriding the default choice (same as in 'XC' keyblock). Specifying 'NONE' will suppress GGA contribution to the SIC energy and potential. Note that very little input checking is done for this key. It is possible to specify a non-existent functional here, which will lead to unpredictable results. Use this key only for debugging.

SELFCONSISTENT {n}

Includes KLI contribution in the XC potential. This is the default. In case of convergence difficulties, v(KLI) will be recomputed until n-th SCF cycle (90 cycles by default). All subsequent SCF cycles will use v(KLI) from the n-th cycle. If a SIC calculation runs into convergence difficulties, it is important to make sure that the SIC energy does not change significantly between the last cycle where potential was computed, and the final SCF cycle.

POSTSCF

Calculates Perdew-Zunger energy correction, using orbitals from standard Kohn-Sham calculation (possibly after localizing them - see LOCALIZE/NOLOCALIZE). Except in a few special cases, such as atoms or the hydrogen molecule, post-SCF SIC corrections are not reliable (see S. Patchkovskii and T. Ziegler, JCP 116, 7806, Ref.[48]).

SKIPCYCLES every {start}

Requests that v(KLI) is recomputed after a given number of SCF iterations. By default, v(KLI) is recomputed on each SCF iteration (every=1). Because v(KLI) evaluation is expensive, using 2 or 3 here may lead to a reduction in computation time. If an optional second parameter is included, v(KLI) is omitted during the first few SCF cycles. This may reduce calculation time if starting guess is not very accurate. Because SIC energy expression is defined as a functional of orbitals, v(KLI) is always omitted in the first SCF cycle, regardless of the "start" setting.

STABLE frac

Specifies mixing coefficient for v(KLI) contribution to exchange-correlation potential. On each iteration where v(KLI) is recomputed, the it is updated as: v(n) = frac*v(KLI) + (1-frac)*v(n-1). The default for frac is 1.0 (i.e. no mixing). Supplying a value smaller than 1 may improve SCF convergence.

NHOMO Nalpha {Nbeta}

Chooses treatment of the free parameter in the KLI approximation (see S. Patchkovskii, J. Autschbach, and T. Ziegler, JCP 115, 26, Ref.[47]). The default is to choose the free parameter such that per-orbital potential shifts are non-negative. If this key is supplied, the per-orbital shifts of Nalpha/Nbeta highest occupied orbitals will be set to zero instead. Both choices are usually identical for converged solutions, but the default exhibits much better convergence behavior.

CORE cor

cor can be one of:
IGNORE: Ignores contributions of frozen core orbitals to SIC energy and potential. Only valence orbitals will appear in KLI potential mixing expression (eq. 14 of [49]).
RHO: Includes core density in the total density of eq.14, but ignore core orbitals otherwise.
POTENTIAL: Includes SIC potential of the core orbitals in v(KLI), but sets corresponding potential shifts to zero.
FULL: Full treatment of the frozen core - frozen core orbitals participate in KLI potential equilibration (eq. 16) on the equal footing with the valence orbitals. This is the default.
For POTENTIAL and FULL, tesseral harmonics are used for the angular part of the core orbitals.

DENSITY mode

This keyword controls evaluation of per-orbital densities in a SIC calculation. mode can be one of:
EXACT: Evaluate densities from molecular orbitals. This is the default.
FIT: Evaluate densities from auxiliary fits. Using this option is not recommended - it is both slower, and less accurate than EXACT.

SHIPV filename

Write v(KLI) contribution to XC potential on grid, to an external file. When using this option, it is a good idea to write out the numerical integration grid as well, by adding "SAVE TAPE10" to the ADF input file.

READV filename

v(KLI) is taken from an external file, and add it to the XC potential. Current integration grid must match the grid used to calculate the potential. The only certain way to guarantee this is to save the integration grid on TAPE10, and pass it around, together with the v(KLI) potential. Specifying this keyword will deactivate add other processing in SIC code, including calculation of Perdew-Zunger energies and SIC potential updates.

READMOS filename {FREEZE}

For the first evaluation of the v(KLI) potential, loads localized orbitals from the specified file, instead of localizing canonical Kohn-Sham MOs. Unless FREEZE is specified, subsequent SCF cycles will use localized canonical MOs.

WRITEMOS filename

Stores localized orbitals to an external file. Reading these orbitals back with READMOS provides a rudimentary restart capability.

NAILCANONICALS {eps}

Stabilize orientation of degenerate canonical MOs prior to localization. The default is not to stabilize. eps is the degeneracy criterion in eV (0.001 by default). Supplying this key may improve convergence when high local symmetry is present in a molecule.

NPFITS n

Number of fit coefficient sets to be computed in a single pass. Large values will improve performance, but need more memory. The default of 15 is usually adequate.