3.28 Dynamical configuration interaction in FHI-aims

Dynamical configuration interaction (DCI) is a quantum embedding theory developed by Marc Dvorak and Patrick Rinke of Aalto University. The implementation in FHI-AIMS is considered stable but still somewhat experimental. Parties interested in doing their own DCI calculations should thoroughly understand the method before reporting any results! While the method has its strengths and the implementation has been quite thoroughly tested, both the theory and implementation are far from the stabile and mature status of something like DFT or $G W$ . Reach out to the developers if questions arise.

Briefly, DCI is a configuration interaction calculation of an energy-dependent effective Hamiltonian. The Hamiltonian is based on the exact Löwdin downfolding of the Schrödinger equation. The energy dependence is treated exactly (not linearized or approximated) by iterating the eigenvalue until it is self-consistent. In an active space (AS) of orbitals, matrix elements are computed using the bare many-body Hamiltonian with the bare Coulomb interaction. This is an exact treatment of the active space. Matrix elements for the remaining degrees of freedom, the bath, are calculated in a $G W$ -BSE-like approximation. Unlike the BSE, which is in a basis of only single excitations, the DCI bath includes all possible multiple excitations (singles, doubles, triples, etc.) up to an arbitrarily high excitation level. A BSE-like method is implemented for these multiple excitations. In this $G W$ -BSE space of the problem, a block diagonal approximation is adopted to dramatically lower the computational cost of the associated matrix inversion. This is where the computational savings lies $-$ many, many inversions of very small matrices instead of one inversion of the full, enormous matrix in the bath space. Even in this simplfied scheme, correlation is described at a $G W$ -BSE-like level of theory in the bath, and it is of course treated exactly in the active space.

What types of problems is DCI good for? DCI is a quantum embedding theory. It is a good choice when the relevant physics is described by a small number ( $\sim$ 2-10) of strongly interacting electrons surrounded by a host medium with weakly-interacting electrons. The strongly interacting electrons form the active space. Localized electrons are a logical choice of active space. Frontier orbitals (HOMO-2, HOMO-1, HOMO, LUMO, LUMO+1, LUMO+2) many times also form a reasonable active space. Defect or impurity systems are a good candidate for DCI. DCI has potential to treat strongly-correlated systems for which $G W$ -BSE fails. Finally, studies of double excitations (or even higher), complicated spin multiplets, or generally multireference systems may all benefit from DCI, which naturally includes all of these effects in the AS.

What can be calculated? Total energies and vertical excitation energies. Absolute values of total energies are not reliable $-$ energy differences (excitation energies or adsorption energies) are much preferred.

DCI is expensive. The largest system studied to date is a single, gas phase free-base porphyrin molecule, and it was challenging. DCI is suitable for diatomics, triatomics, small organic molecules, etc. $-$ small stuff. Active spaces of no more than 8 levels and a bath of no more than $\sim$ 300 orbitals are recommended.

How are calculations performed? The computational workflow is: Hartree-Fock, then RPA correlation energy, then $G W$ , then DCI. Output from these preliminary calculations is used as input for DCI. The workflow is done internally in AIMS and the entire calculation is controlled from control.in with several input variables. There are other input variables related to $G W$ or RPA calculations which also need to be specified in order to produce the matrix elements used as input for DCI.

Theory

The theory is covered in detail in the main DCI publications [83, 82], which should be cited with any DCI calculations you perform. The Löwdin downfolding reduces the full Schrödinger equation to an energy-dependent one existing in a chosen subspace.

H^{\mathrm{eff}}(E)=DHD+DHR\frac{1}{E-RHR}RHD

(3.102)

Here, $D$ and $R$ are projection operators that project onto the many-body Hilbert space and $E$ is the chosen eigenvalue of the Schrödinger equation, to be determined self-consistently. $D$ is the set of all configurations generated by active space excitations. The Hamiltonian is the interacting electronic Hamiltonian,

H=\sum_{ij}t_{ij}+\frac{1}{2}\sum_{ijkl}v_{ijkl}a_{i}^{\dagger}a_{j}^{\dagger}% a_{l}a_{k}\;.

(3.103)

The expense of solving Eq. 3.102 is in inverting the $R H R$ block of the Hamiltonian.

DCI rewrites this block of the Hamiltonian as a ground state energy ( $E_{0}^{R}$ ) plus excitation energies ( $\Omega^{R}$ ):

RHR\rightarrow H^{R}=E_{0}^{R}+\Omega^{R}

(3.104)

The ^R superscript on these quantities indicates that their calculation needs to be constrained within the $R$ subspace. The multiple excitations in $\Omega^{R}$ are treated with $G W$ and the screened interaction between quasiparticles, analogous to the Bethe-Salpeter equation. For a configuration of excitation level $m$ ( $m=1,2,3,...$ ), the excitation energy $\Omega^{R}$ contains $2m$ (factor of 2 for electrons and holes) $G W$ quasiparticles. Additionally, we include all possible pairwise interactions between quasiparticles via the screened Coulomb interaction (electron-hole, electron-electron, and hole-hole). To constrain the calculation to the $R$ subspace, as indicated by ^R in Eq. 3.104, we apply the constrained Random Phase Approximation (cRPA) everywhere that the polarizability enters the calculation. Active space transitions are omitted when calculating the polarizability within the cRPA. This constraint correctly satisfies the limits of the embedding: we recover a full CI calculation as $D\rightarrow\mathbb{I}$ and an ordinary many-body perturbation theory calculation as $R\rightarrow\mathbb{I}$ . Furthermore, if we restrict the basis to single excitations ( $m=1$ ) in the limit $R\rightarrow\mathbb{I}$ and select bare electron-hole exchange, we recover the ordinary Bethe-Salpeter equation. Very nice!

The scalar $E_{0}^{R}$ must also be calculated. The best approach is to calculate a correlation energy based on the cRPA. The calculation of RPA correlation energies is implemented in AIMS; it is a trivial modification to calculate the same quantity with a modified polarizability according to the cRPA. The cRPA correlation energy is added to the energy of the reference configuration (Hartree-Fock energy) to define $E_{0}^{R}$ . This RPA calculation is therefore a prerequisite before entering a DCI calculation, along with Hartree-Fock and $G W$ .

$\Omega^{R}$ is, in principle, a matrix, and an enormous one $-$ one the size of FCI. The major simplification of DCI is that this matrix is treated as block diagonal, greatly simplifying the computational expense of the inversion. Blocks of $\Omega^{R}$ are constructed in the spin subspaces for a given excitation to correctly capture the spin symmetry (for a single excitation, describe either a singlet or triplet state). The configuration basis of DCI is not an eigenbasis of $S^{2}$ , hence the need for the blocks in the spin subspaces to describe either singlet or triplet states.

Notes about implementation and running calculations

DCI can currently only run from a Hartree-Fock mean-field. DFT functionality is not enabled. A preliminary $G_{0}W_{0}$ calculation is always necessary, and a preliminary RPA calculation for the correlation energy with the constrained polarizability is necessary 99% of the time. It is strongly encouraged to use as many of the default DCI input variables as possible. RPA correlation is the default. Many of the DCI input variables do not need to be changed and, in fact, should not be changed for best results! One should, in principle, specify only the active space, the excitation level up to which to consruct AS configurations, and optionally parameters related to convergence and mixing strength. DCI can only be run for nonperiodic systems.

The memory requirements for DCI can be very high due to the need to store two different sets of 4-center integrals (bare and screened). To accommodate this, it is often necessary for a single process to use the entire memory available on a full node. The code uses an MPI shared window to accomplish this without wasting extra cores. All the cores on a single node access a single copy of the 4-center integrals on that node. This is taken care of automatically and should be transparent to the user. Improved memory management and performance is a topic of future work.

In addition to standard input variables controlling SCF convergence, etc., typical parameters for a DCI calculation are:

  xc                    hf

  qpe_calc              gw
  anacon_type           two-pole

  total_energy_method   rpa

  neutral_excitation    dci
  dci_tol               1.0d-5
  dci_states            1 2
  dci_exc_level         6
  dci_active_space      5 10 5 10
  dci_dynamic           1 20 1 20

  use_2d_corr           .false.

Note that use_2d_corr must be switched to .false.. Only the old, 1D distribution is currently compatible with DCI (upgrading this to support 2D distribution is a future goal). This is not as limiting to the overall calculation as it seems. The memory and CPU requirements for DCI are much larger than other parts of the calculation, and so restricting these portions to a 1D distribution is not a limiting factor.

This control.in will compute the ground state and first excitation energy with DCI (specified by dci_states). Occupied and empty states for DCI do not need to be specified separately. AIMS will find them based on the chemical potential. This active space includes all orbitals from 5 to 10, which is embedded in the set of orbitals from 1 to 20. Active space configurations are constructed up to excitation level 6, equivalent to FCI for this 6-level active space at half-filling (if HOMO is orbital 7). This is the minimal set of input variables needed for DCI.

The easiest way to choose the active space is by selecting frontier orbitals or by visual inspection of the wave functions. One can also run a BSE calculation and examine a specific eigenvector of interest there. One can then run a DCI calculation with the dominant BSE transitions placed in the DCI active space in order to calculate the same state with DCI. It is up to the user to choose a meaningful active space.

Tags for general section of control.in:

Tag: neutral_excitation(control.in)

Usage: neutral_excitation dci
Purpose: Turns on the dynamical configuration interaction (DCI) calculation. DCI runs as post-processing after the Hartree-Fock mean-field calculation and after $G W$ .

Tag: dci_tol(control.in)

Usage: dci_tol real
Purpose: Numerical tolerance on the self-consistency of the DCI total energy. The Hamiltonian is iterated until the output eigenvalue is the same as the energy at which the Hamiltonian is constructed, within the numerical tolerance specified by dci_tol.
real: real number in units of Ha. Default: 1.0d-5

Tag: dci_mix(control.in)

Usage: dci_mix real
Purpose: Mixing value for a simple linear mixing between old and new eigenvalues of the DCI Hamiltonian. The new solution takes the weight dci_mix while the old solution has weight $(1-$ dci_mix).
real: mixing parameter, real number. Default: 0.65

Tag: dci_exc_level(control.in)

Usage: dci_exc_level integer
Purpose: Maximum excitation level for configurations in the active space. It is up to the user to determine what that level should be if, for example, they want full configuration interaction. Exact diagonalization or FCI in the AS is not necessary, and an excitation level of no less than triples (dci_exc_level 3 is recommended. Maximum allowed is 8.
integer: maximum excitation level Default: 3

Tag: dci_active_space(control.in)

Usage: dci_active_space int1 int2 int3 int4
Purpose: Specify the active orbitals for DCI by an orbital range. In the first spin channel, all orbitals between int1 and int2 (inclusive) are included in the AS. In the second spin channel, all orbitals between int3 and int4 are included. Only one of dci_active_space or dci_active_select_XX must be specified.
int1 int2 int3 int4: 4 integers specifying active space Default: none, must be specified by user

Tag: dci_dynamic(control.in)

Usage: dci_dynamic int1 int2 int3 int4
Purpose: Specify the dynamic orbitals for DCI by an orbital range. In the first spin channel, all orbitals between int1 and int2 (inclusive) are included in the AS. In the second spin channel, all orbitals between int3 and int4 are included. dci_dynamic specifies a much larger range of orbitals, those which are correlated at the $G W$ /BSE-like level of theory. These are the bath orbitals of the DCI theory and should, in principle, include the full basis set, but the convergence with respect to bath orbitals specified by dci_dynamic can be tested. This range is not related at all to the preceding $G W$ or RPA calculations; these methods have their own input variables setting the orbitals used, and we recommend those calculations before DCI always be performed with the full basis set in all parts. It is reasonable to omit core states from dci_dynamic.
int1 int2 int3 int4: 4 integers specifying range of dynamic orbitals. Default: none, must be specified by user

Tag: dci_active_select_up(control.in)

Usage: dci_active_select_up integer array
Purpose: Specify the active orbitals for DCI by selecting individual states by their orbital number. dci_active_select_up is a list of space-separated integers (74 77 78 84), including both occupied and unoccupied states, that specifies orbitals in the up spin channel for the active space. If used, this input variable must be paired with dci_active_select_dn. Only one of dci_active_space or dci_active_select_XX must be specified.
integer array: list of integers, individually chosen active orbitals as int1 int2 int3 int4 ... Default: none, must be specified by user

Tag: dci_active_select_dn(control.in)

Usage: dci_active_select_dn integer array
Purpose: Specify the active orbitals for DCI by selecting individual states by their orbital number. dci_active_select_dn is a list of space-separated integers (74 77 78 84), including both occupied and unoccupied states, that specifies orbitals in the down spin channel for the active space. If used, this input variable must be paired with dci_active_select_up. Only one of dci_active_space or dci_active_select_XX must be specified.
integer array: list of integers, individually chosen active orbitals as int1 int2 int3 int4 ... Default: none, must be specified by user

Tag: dci_states(control.in)

Usage: dci_states integer array
Purpose: Specify the states to be computed with DCI. The ground state is 1 and must always be included first in the list. The first excited state is 2 and so on. Selected states need not be consecutive. For example, 1 3 4 5 8. Maximum number of states is 20.
integer array: States to compute with DCI, integer list. Default: none, must be specified by user

Tag: dci_energy_start(control.in)

Usage: dci_energy_start real array
Purpose: Specify the energies at which self-consistent DCI cycles are initialized. This functionality is essential for restart capability. If the calculation terminates, simply take the last computed energy from the DCI cycle and enter it here. The new calculation will start iterations from this value. The energy can be taken from the line Evaluation energy X for excitation Y in the output file. NOTE! The ground state energy is always required (specified as 1 in dci_states) and its value is entered in dci_energy_start as a total energy (Ha). All remaining entries for dci_energy_start are entered as excitation energies (Ha), defined as total energy differences between the state of interest and the internal DCI ground state.
real array: Starting energies for DCI iterations, real list. Default: none, must be specified by user

Tag: dci_write_els(control.in)

Usage: dci_write_els boolean
Purpose: Choose to write out to disk all the matrix elements needed to restart a DCI calculation. This includes bare and screened Coulomb matrix els, one-body matrix elements, $GW^{R}$ eigenvalues based on cRPA, Hartree-Fock eigenvalues, and a small number of additional parameters. Because of the $N^{4}$ scaling of the Coulomb matrix elements, this can consume a large amount of disk space. Similar to FCIDUMP.
boolean: .true. or .false. Default: .false.

Tag: dci_read_els(control.in)

Usage: dci_read_els boolean
Purpose: Choose to read from disk all the matrix elements needed to restart a DCI calculation. This includes bare and screened Coulomb matrix els, one-body matrix elements, $GW^{R}$ eigenvalues based on cRPA, Hartree-Fock eigenvalues, and a small number of additional parameters. There are no checks about the consistency of these elements with what is specified in control.in $-$ the user must be sure they are reading the correct elements!
boolean: .true. or .false. Default: .false.

The remaining input variables control more advanced features. Only use these if you are sure of what you are doing. It is recommended that most of them remain at their defaults. A complete list of DCI variables, including additional input variables for control.in not listed here and internal variables, is in the src file dci_vars.f90.

Tag: dci_eh_scr_ex(control.in)

Usage: dci_eh_scr_ex boolean
Purpose: Turn on or off the screened exchange between electron and hole when calculating the excitation energies $\Omega^{R}$ . This choice applies to the exchange interaction between all electron-hole pairs in a multiple excitation. .true. uses screened exchange and is the default choice. This was found to perform better even though the BSE uses bare e-h exchange.
boolean: .true. or .false. Default: .true.

Tag: dci_scissor(control.in)

Usage: dci_scissor real
Purpose: Scissor shift applied to the mean-field band structure immediately before starting the quasiparticle ( $G W$ ) calculation. All unoccupied states are shifted by this amount. You must enter a negative value to close the HOMO-LUMO gap. The shifted band structure is then used in the full DCI calculation: it enters the Green’s function, the polarizability, the $G W$ self-energy, and the calculation of the cRPA correlation energy. This does NOT affect the one-body matrix elements of the bare many-body Hamiltonian, only the quasiparticle portion of the DCI calculation.
real: scissor shift to unoccupied states in units of Ha, real number. Default: 0.0

Tag: dci_type(control.in)

Usage: dci_type integer
Purpose: Choose the DCI algorithm based on either the original theory published in [83] (dci_type 1), a ground state energy in $R$ based on a constrained random phase approximation (cRPA) correlation energy (dci_type 2), or an approximation based on the diagonal, bare many-body Hamiltonian (dci_type 3). Option 2 is by far the most preferred; option 1 should never be used; option 3 can be useful as a comparison with option 2. The difference between options 2 and 3 is the effect of screening in the $R$ subspace.
integer: Version of DCI theory. Default: 2

Tag: dci_opt(control.in)

Usage: dci_opt integer
Purpose: Different options for parallelism, memory management, and handling of block diagonal construction in a DCI calculation. dci_opt 6 is the preferred choice but has the highest memory demands. dci_opt 7 is a slower version but requires less memory. These two chocies are the only ones that should be considered.
integer: Choice of DCI computational algorithm. Default: 6

Tag: dci_dist(control.in)

Usage: dci_dist integer
Purpose: Different options for parallel distribution to compute matrix elements of DCI Hamiltonian. dci_dist 1 is the default. dci_dist 2 is a redistribution of matrix elements meant to reduce the memory requirements when combined with dci_opt 6. dci_dist 3 is also available. dci_dist 2 has noticeably lower memory requirements than dci_dist 1 when combined with dci_opt 6. The speed and memory performance of other combinations largely depends on the system.
integer: Distribution of DCI matrix elements. Default: 1