7.21. Dynamic Correlation Dressed CAS¶

DCD-CAS(2) is a post-CASSCF MRPT method of the perturb-then-diagonalize kind, i.e. it can modify the CAS wavefunction compared to the previous CASSCF.[652] In cases where CASSCF already provides a good 0th order wavefunction, DCD-CAS(2) energies are comparable to NEVPT2.

7.21.1. Theory of Nonrelativistic DCD-CAS(2)¶

The DCD-CAS(2) method is based on solving the eigenvalue problem of an effective Hamiltonian of the form

\[H_{IJ}^{\text{DCD}, S} = \langle \Phi_I^{SS} | H | \Phi_J^{SS} \rangle - \sum_{K \in \text{FOIS} } \frac{\langle \Phi_I^{SS} | H | \tilde{\Phi}_K^{SS} \rangle \langle \tilde{\Phi}_K^{SS} | H | \Phi_J^{SS} \rangle }{E_K^S - E_0^S}\]

for each total spin \(S\) separately. The 0th order energies \(E_K^S\) of the perturbers \(|\tilde{\Phi}_K^{SS}\rangle\) are obtained by diagonalizing the Dyall’s Hamiltonian in the first-order interacting space (FOIS). The effective Hamiltonian has the form of a CASCI Hamiltonian that is dressed with the effect of dynamic correlation (dynamic correlation dressed, DCD), hence the name for the method. \(E_0^S\) is chosen to be the ground state CASSCF energy for the respective total spin \(S\). Since this choice is worse for excited states than for the ground state, excitation energies suffer from a “ground state bias”.

For the contribution coming from perturbers in which electrons are excited from two inactive (\(ij\)) to two virtual (\(ab\)) orbitals, we use (when writing the DCD Hamiltonian in a basis of CASCI states) the alternative expression

\[\langle \Psi_I^{SS} | H^\text{DCD}(ij \rightarrow ab) | \Psi_J^{SS} \rangle = - \delta_{IJ} E_\text{MP2}\]

\[E_\text{MP2} = \sum_{ijab} \frac{(ib|ja)^2 - (ib|ja)(ia|jb) + (ia|jb)^2}{\epsilon_a + \epsilon_b - \epsilon_i - \epsilon_j}\]

Since in this version the \(ij\rightarrow ab\) perturber class does not contribute at all to excitation energies (like it is assumed in the difference-dedicated configuration interaction method), we call this the difference-dedicated DCD-CAS(2) method. Since the \(ij\rightarrow ab\) class contributes the largest part of the dynamic correlation energy, this also removes the largest part of the ground state bias. This option is used as default in DCD-CAS(2) calculations. In order to also remove the ground state bias from the other perturber classes, we furthermore apply a perturbative correction to the final energies. At first order (which is chosen as default), it takes the form

\[\Delta E_I = - \Delta_I \sum_{K \in \text{FOIS} } \frac{\langle \tilde{\Psi}_I | H | \tilde{\Phi}_K\rangle \langle \tilde{\Phi}_K | H | \tilde{\Psi}_I \rangle }{(E_K - E_0)^2}\]

\[\Delta_I = \langle \tilde{\Psi}_I | H | \tilde{\Psi}_I \rangle - E_0\]

for the correction \(\Delta E_I\) to the total energy of the \(I\)th DCD-CAS(2) root \(|\tilde{\Psi}_I\rangle\).

7.21.2. Treatment of spin-dependent effects¶

The theory so far is valid for a nonrelativistic or scalar-relativistic Hamiltonian \(H\). If we modify it to a Hamiltonian \(H+V\), where \(V\) contains effects that are possibly spin-dependent, this leads us to a theory [491] which has a similar form as QDPT with all CAS roots included. The form of the spin-dependent DCD-CAS(2) effective Hamiltonian is

\[\langle \Phi_I^{SM} | H^\text{DCD} | \Phi_J^{S'M'}\rangle = \delta_{SS'} \delta_{MM'} H_{IJ}^{\text{DCD},S,\text{corr} } + \langle \Phi_I^{SM} | V | \Phi_J^{S'M'} \rangle.\]

\[\mathbf{H}^{\text{DCD}, S, \text{corr} } = \mathbf{C}^\text{DCD} \mathbf{E} (\mathbf{C}^\text{DCD})^T.\]

In order to construct it, we first need to solve the scalar-relativistic DCD-CAS(2) problem to construct the matrix \(\mathbf{H}^{\text{DCD},S,\text{corr} }\) from the bias corrected energies \(\mathbf{E}\) and DCD-CAS(2) CI coefficients \(\mathbf{C}\) and then calculate the matrix elements of the operators contributing to V in the basis of CSFs \(|\Phi_I^{SM}\rangle\).

Zero field splitting D tensors are extracted using the effective Hamiltonian technique, i.e. fitting the model Hamiltonian to a des-Cloiseaux effective Hamiltonian that is constructed from the relativistic states and energies by projection onto the nonrelativistic multiplet (see section Zero-Field Splitting and the reference [565]). There are limitations to this approach if spin orbit coupling becomes so strong that the relativistic states cannot uniquely be assigned to a single nonrelativistic spin multiplet.

Hyperfine A-matrices and Zeeman g-matrices for individual Kramers doublets consisting of states \(|\Phi\rangle, |\overline{\Phi}\rangle\) are extracted by comparing the spin Hamiltonians

\[H_\text{Zeeman} = \mu_B \vec B \cdot g \cdot \vec S\]

\[H_\text{HFC} = \sum_A \vec I^A \cdot A^A \cdot \vec S\]

to the matrix representation of the many-electron Zeeman and HFC operators in the basis of the Kramers doublet. This yields [491]

\[\begin{split}\begin{aligned} g_{k1} &= 2\Re \langle \overline{\Phi} | L_k + g_e S_k | \Phi \rangle \\ g_{k2} &= 2\Im \langle \overline{\Phi} | L_k + g_e S_k | \Phi \rangle \\ g_{k3} &= 2\langle \Phi | L_k + g_e S_k | \Phi \rangle\end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} A_{k1} &= -2\gamma_A \Re \langle \overline{\Phi} | B_k^\text{HFC}(\vec R_A) | \Phi \rangle \\ A_{k2} &= -2\gamma_A \Im \langle \overline{\Phi} | B_k^\text{HFC}(\vec R_A) | \Phi \rangle \\ A_{k3} &= -2\gamma_A \langle \Phi | B_k^\text{HFC}(\vec R_A) | \Phi \rangle\end{aligned}\end{split}\]

where \(B_k^\text{HFC}(\vec R_A)\) is the \(k\)th component of the magnetic hyperfine field vector at the position of nucleus \(A\) and \(\gamma_A\) is the gyromagnetic ratio.

7.21.3. List of keywords¶

The following keywords can be used in conjunction with the DCD-CAS(2) method:

%casscf
  dcdcas true                # Do a DCD-CAS(2) calculation
  diffded true               # Use difference-dedicated DCD-CAS(2) for the 
                             # ij->ab class
  corrorder 1                # Maximum order for the perturbative bias correction
                             # (lower orders come for free)
  dcd_ritrafo false          # Use RI approximation for electron repulsion integrals

  dcd_soc false              # Relativistic DCD-CAS(2) with spin orbit coupling
  dcd_ssc false              # Relativistic DCD-CAS(2) with direct electronic
                             # spin-spin coupling
  dcd_domagfield 0           # Number of user-specified finite magnetic fields
  dcd_dtensor false          # Calculate an effective Hamiltonian D-tensor
  dcd_nmultd 1               # The number of nonrelativistic multiplets for which the
                             # D-tensor is calculated
  dcd_gmatrix false          # Calculate an effective Kramers pair Zeeman g-matrix
  dcd_amatrix false          # Calculate an effective Kramers pair Hyperfine A-matrix
  dcd_kramerspairs 1         # The number of Kramers pairs for which g and/or A
                             # is calculated
  dcd_magnetization false    # Calculate the magnetization of the molecule in an
                             # external magnetic field

  dcd_cascimode false        # Run relativistic calculation in CASCI mode, i.e.
                             # without the dynamic correlation dressing
  dcd_natorbs false          # Calculate natural orbitals for each state and write
                             # them to disk
  dcd_populations false      # Perform population analysis on the DCD-CAS(2) states
end

Note that the calculation of SSC requires the definition of an auxiliary basis set, since it is only implemented in conjunction with RI integrals. If dcd_magnetization is requested, the values for magnetic flux density and temperature to be used can be specified via the keywords MAGTemperatureMIN, MAGTemperatureMAX, MAGTemperatureNPoints, MAGFieldMIN, MAGFieldMAX, MAGNpoints of the rel subblock of the %casscf block (see section Magnetization and Magnetic Susceptibility). If the keyword dcd_domagfield is set to a number different than 0, the magnetic fields can be entered as a matrix of xyz coordinates (in Gauss), e.g.

%casscf
  dcdcas true
  ...
  dcd_domagfield 2
  dcd_magneticfields[0] = 10000, 0, 0
  dcd_magneticfields[1] = 0, 10000, 0
end

Furthermore, there is the keyword DCD_EDIAG that when running the DCD-CAS(2) code in CASCI mode works analogously to the keyword EDiag in the soc subblock of the %mrci block (see section Zero-Field Splitting). The only difference is that the energies should be entered in atomic units, not in wavenumbers, e.g.

%casscf
  ...
  dcdcas true
  dcd_cascimode true
  dcd_soc true
  DCD_EDIAG[0] -2220.920028
  DCD_EDIAG[1] -2220.834377
  DCD_EDIAG[2] -2220.835871
  DCD_EDIAG[3] -2220.810290
  DCD_EDIAG[4] -2220.812293
  DCD_EDIAG[5] -2220.756732
end