3.12. Correlated Methods using Automatic Code Generation (AUTOCI)

The AUTOCI module hosts a number of methods, where their manual implementation would be tedious or practically impossible. The module works with all types of reference wave function available in ORCA, i.e., RHF, ROHF, UHF and CASSCF and offers CI and related methods.

All the methods are implemented in canonical orbital basis and storing all integrals on disk.

3.12.1. Introduction

All the theories are obtained by the means of automated programming within the ORCA-AGE (Automated Generator Environment for ORCA).[359, 360] The CI module reads in the SCF wavefunction and optimizes the coefficient of the CI expansion. Conceptually, the module is similar to orca_mdci, therefore the input and output do have a lot in common.

3.12.2. Basis Usage / List of Features

The following single-reference methods are currently implemented in the AUTOCI.

Reference	Correlation	Energy	Gradient
RHF	CID	✓	✓
RHF	CISD	✓	✓
RHF	CISDT	✓
RHF	CCD	✓	✓
RHF	CCSD	✓	✓
RHF	CCSD[T]	✓	✓
RHF	CCSD(T)	✓	✓
RHF	CCSDT	✓
RHF	CEPA(0)	✓	✓
RHF	CC2	✓
RHF	QCISD	✓
RHF	MP2	✓	✓
RHF	MP3	✓	✓
RHF	MP4	✓	✓
RHF	MP4(SDQ)	✓	✓
UHF	CID	✓	✓
UHF	CISD	✓	✓
UHF	CISDT	✓	✓
UHF	CCD	✓	✓
UHF	CCSD	✓	✓
UHF	CCSD[T]	✓	✓
UHF	CCSD(T)	✓	✓
UHF	CCSDT	✓	✓
UHF	CCSDTQ	✓
UHF	CCSDT-1	✓
UHF	CCSDT-2	✓
UHF	CCSDT-3	✓
UHF	CCSDT-4	✓
UHF	CEPA(0)	✓	✓
UHF	CC2	✓
UHF	CC3	✓
UHF	QCISD	✓
UHF	MP2	✓	✓
UHF	MP3	✓	✓
UHF	MP4	✓	✓
UHF	MP4(SDQ)	✓	✓
UHF	MP5	✓	✓
ROHF	CISD	✓

Any AUTOCI single-reference methods can be called from the simple keyword by prepending AUTOCI- to the correlation method, for instance

! AUTOCI-CCSD

The prefix ‘AUTOCI’ is added to dinstinguish from the handwritten implementations in orca_mdciand orca_mp2.

Warning

Higher order CC calculations like CCSDT and CCSDTQ have extreme memory requirements in the CC Iterations. For UHF CCSDTQ particularly (computational efford scales as \(N^{10}\) !) only systems with few electrons and basis functions are feasible depending on your achitecture.

For a CASSCF reference wavefunction, the prefix is omitted. The following multi-reference methods and their unrelaxed densities are available.

Reference	Correlation	Energy	Unrelaxed Density
CASSCF	FIC-MRCI	✓	✓
CASSCF	FIC-MRCEPA(0)	✓	✓
CASSCF	FIC-MRACPF	✓	✓
CASSCF	FIC-MRAQCC	✓	✓
CASSCF	FIC-MRCC	✓	Linearized
CASSCF	FIC-DDCI3	✓	✓
CASSCF	FIC-NEVPT2	✓	✓
CASSCF	FIC-NEVPT3	✓	✓
CASSCF	FIC-NEVPT4(SD)	✓	✓

Note

If a multi-reference method is specified as a simple keyword, all information about reference spaces, number of roots etc. is taken from the CASSCF module that is assumed to be run in advance. More refined settings require the %autoci block in the input.

3.12.3. Analytic Nuclear Gradients with AUTOCI

Obtaining accurate geometries is crucial to computing molecular properties accurately. In order to perform geometry optimisations, the nuclear gradient is necessary and while this can easily be obtained using numerical finite difference methods, it is also quite costly. More importantly, perhaps, is the fact that numeric derivatives tend to become unstable. Therefore, being able to evaluate analytic gradients is of vital importance. Using the AGE, a general framework has been built that supports arbitrary-order CI, CC and MPn nuclear gradients (and other derivatives).[360]

An example is shown below how to optimise a geometry using AUTOCI’s gradients at the CCSD level of theory

! RHF cc-pVTZ AUTOCI-CCSD VerytightSCF Opt

%maxcore 10000

*xyz 0 1
...
*

The analytic gradients can even be used to perform semi-numerical frequency calculations

! AUTOCI-CCSD NumFreq

Besides nuclear gradients, all other first-order properties available in ORCA are available for the respective methods, such as dipole/quadrupole moments, hyperfine couplings or quadrupole splittings. As discussed above, (un)relaxed densities can be requested via

%autoci
  density relaxed
end

For geometry optimisations, both the unrelaxed and relaxed densities are computed automatically and do not need to be requested explicitly.

3.12.4. AUTOCI Response Properties via Analytic Derivatives

For single-reference methods (CCSD and CCSD(T)), some response properties could be calculated via the analytic derivative scheme.

Response Property	Reference	Correlation	Orbital Response
Static Polarizabilities*	RHF	CCSD	Unrelaxed / Relaxed
	RHF	CCSD(T)	Unrelaxed / Relaxed
	UHF	CCSD	Unrelaxed / Relaxed
	UHF	CCSD(T)	Unrelaxed / Relaxed
NMR Shielding (GIAO)	RHF	CCSD	Unrelaxed / Relaxed
	RHF	CCSD(T)	Unrelaxed / Relaxed
	UHF	CCSD	Unrelaxed / Relaxed
	UHF	CCSD(T)	Unrelaxed / Relaxed
G-Tensor (GIAO)	UHF	CCSD	Unrelaxed / Relaxed
	UHF	CCSD(T)	Unrelaxed / Relaxed

(* dipole-dipole, dipole-quadrupole, quadrupole-quadrupole) Note that for NMR shielding and EPR g-Tensor the fully analytic GIAO approach (see Nuclear Magnetic Resonance (NMR) Parameters and Gauge origin treatment) is the default and the common gauge approach (less accurate) is not available.

To calculate these properties, the following quantities need to be computed:

• Unperturbed CC amplitudes

• Unperturbed CC \(\Lambda\) (Lagrange Multiplier)

• Unperturbed 1-RDM

• Unperturbed 2-RDM (\(\dagger\))

• Unperturbed Z-Vector (\(\dagger\))

• Perturbed CC amplitudes

• Perturbed CC \(\Lambda\) (Lagrange Multiplier)

• Perturbed 1-RDM

• Perturbed 2-RDM (\(\dagger\))

• Perturbed Z-Vector (\(\dagger\))

Tasks labelled with (\(\dagger\)) are optional - they are only required for the “relaxed” scheme of computing response properties. With the “relaxed” scheme, the orbitals are allowed to relax to external perturbations. (See: orbital relaxation

3.12.4.1. Unperturbed Parameters

From CC gradient theory and as stated in the previous MDCI section, the CC energy Lagrangian could be written as (denoting the similarity-transformed Hamiltonian as \(\bar{H}=e^{-T}H e^{T}\)):

\[ L _{\text{CC}} = \langle \Phi_{0}|\bar{H}|\Phi_{0}\rangle + \sum_{\eta}\lambda_{\eta}\langle \Phi_{\eta}|\bar{H}|\Phi_{0}\rangle + \sum_{ai}f _{ai}z _{ai} \]

Equivalently, we could also express this Lagrangian in terms of density matrices:

\[\begin{split} L _{\text{CC}} = \sum_{pq}h _{pq} \frac{\partial L_{\text{CC}}}{\partial h_{pq}} \\ + \sum_{pqrs}(pq|rs)\frac{\partial L_{\text{CC}}}{\partial (pq|rs)} \\ = \sum_{pq}h _{pq}D_{pq} + \sum_{pqrs}(pq|rs)d_{pqrs} \end{split}\]

in which the one- and two-body reduced density matrices (RDMs) could be decomposed into the amplitude and orbital contributions:

\[\begin{split} D_{pq} = \gamma_{pq} + z _{pq}\delta_{pa}\delta_{qi} \\ d_{pqrs} = \Gamma^{pr}_{qs} + z_{pq}\sum_{k}\left(\delta_{pa}\delta_{qi}\delta_{rk}\delta_{sk} -\delta_{pa}\delta_{qk}\delta_{rk}\delta_{si}\right) \end{split}\]

The above are known as the relaxed 1- and 2-RDMs, since the terms involving the \(z_{pq}\) arise due to orbital relaxation. The CC amplitude contributions (a.k.a. the unrelaxed 1- and 2-RDMs) are defined as:

\[\begin{split} \gamma_{pq} = \langle \Phi_{0}|(1+\Lambda)e^{-T}E^{p}_{q}e^{T}|\Phi_{0} \rangle \\ \Gamma^{pr}_{qs} = \langle \Phi_{0}|(1+\Lambda)e^{-T}\left(E^{p}_{q}E^{r}_{s}-E^{p}_{s}\delta_{rq}\right)e^{T}|\Phi_{0} \rangle \end{split}\]

These quantities depend on the CC amplitudes \(t_{\eta}\) and Lagrange multipliers \(\lambda_{\eta}\), for which we need to solve the amplitude- and \(\Lambda\)-equations:

\[\begin{split} \langle \Phi_{\eta}|\bar{H}| \Phi_{0}\rangle = 0 \\ \langle \Phi_{0}|(1+\Lambda)(\bar{H}-E_{\text{CC}})|\Phi_{\eta}\rangle = 0 \end{split}\]

3.12.4.2. Perturbed Parameters

The well-known Wigner’s \(2n+1\) rule of perturbation theory states: wavefunction parameters up to the \(n\)-th order are required for the energies at the order of \(2n+1\). For non-variational methods like CC, the Lagrangian approach is employed, resulting in a variational energy functional (albeit non-Rayleigh-Ritz, thus no guaranteed upper bound), which enables the analytic evaluation of derivatives. The \(2n+2\) rule of the Lagrange multiplier in such Lagrangian functional states: Lagrange multipliers up to the \(n\)-th order are required for the energies at the order of \(2n+2\).

In practice, and here in ORCA-AUTOCI implementation, these two rules (especially the \(2n+2\) rule) are not fully exploited. This is due to the consideration of computational efficiency for certain molecular propeties. For example, the general expression for the second derivative of the CC energy, subject to two external perturbations \(x\) and \(y\), could be written as:

\[\begin{split} \frac{\mathrm{d} ^{2}E}{\mathrm{d} x \mathrm{d} y} = \langle \Phi_{0}|(1+\Lambda)e^{-T}\frac{\partial ^{2}H}{\partial x \partial y}e^{T}|\Phi_{0}\rangle \\ + \langle \Phi_{0}|(1+\Lambda)\left[e^{-T}\frac{\partial H}{\partial x}e^{T},\frac{\mathrm{d} T}{\mathrm{d} y}\right]|\Phi_{0}\rangle \\ + \langle \Phi_{0}|(1+\Lambda)\left[e^{-T}\frac{\partial H}{\partial y}e^{T},\frac{\mathrm{d} T}{\mathrm{d} x}\right]|\Phi_{0}\rangle \\ + \langle \Phi_{0}|(1+\Lambda)\left[\left[e^{-T}H e^{T}, \frac{\mathrm{d} T}{\mathrm{d} x}\right],\frac{\mathrm{d} T}{\mathrm{d} y}\right]|\Phi_{0}\rangle \end{split}\]

which is known as the “symmetric” expression, satisfying both the \(2n+1\) and \(2n+2\) rules. However, the following “asymmetric” expression might be preferred for mixed derivatives:

\[\begin{split} \frac{\mathrm{d} ^{2}E}{\mathrm{d} x \mathrm{d} y} \\ = \langle \Phi_{0}|(1+\Lambda)e^{-T}\frac{\partial ^{2}H}{\partial x \partial y}e^{T}|\Phi_{0}\rangle \\ + \langle \Phi_{0}|(1+\Lambda)\left[e^{-T}\frac{\partial H}{\partial x}e^{T},\frac{\mathrm{d} T}{\mathrm{d} y}\right]|\Phi_{0}\rangle \\ + \langle \Phi_{0}|\frac{\mathrm{d} \Lambda}{\mathrm{d}y}e ^{-T}\frac{\partial H}{\partial x}e ^{T}|\Phi_{0}\rangle \end{split}\]

In this way, we solve the perturbed \(T\) and \(\Lambda\) equations subject to only a single perturbation. For NMR shielding this means we can solve these equations only for the magnetic-field perturbation and there are \(3\) equations for \(\mathrm{d} T/\mathrm{d} B_{i}\) and \(3\) for \(\mathrm{d} \Lambda/\mathrm{d} B_{i}\). If we were to use the first “symmetric” expression, although we could save some effort by not having to solve for any perturbed \(\Lambda\), we still need to solve for perturbed \(T\), subject to both the magnetic-field perturbation (\(3\) equations) and the nuclear magnetic moment perturbation (\(3N_{\text{atoms}}\) equations).

To obtain the equations for perturbed \(T\) and \(\Lambda\), we differentiate the unperturbed equations, w.r.t. external perturbation \(\chi\)

\[\begin{split} \frac{\partial}{\partial\chi} \langle \Phi_{\eta}|\bar{H}|\Phi_{0}\rangle = 0 \\ \frac{\partial}{\partial\chi} \langle \Phi_{0}|(1+\Lambda)(\bar{H}-E_{\text{CC}})|\Phi_{\eta}\rangle = 0 \end{split}\]

Then we build the perturbed 1- and 2-RDMs, by taking derivative of the unperturbed expressions:

\[\begin{split} \frac{\partial\gamma_{pq}}{\partial\chi} = \frac{\partial}{\partial\chi}\langle \Phi_{0}|(1+\Lambda)e^{-T}E^{p}_{q}e^{T}|\Phi_{0} \rangle \\ \frac{\partial\Gamma^{pr}_{qs}}{\partial\chi} = \frac{\partial}{\partial\chi}\langle \Phi_{0}|(1+\Lambda)e^{-T}\left(E^{p}_{q}E^{r}_{s}-E^{p}_{s}\delta_{rq}\right)e^{T}|\Phi_{0} \rangle \end{split}\]

The expressions and relevant modules are produced by the code generator. Note that for the currently available response properties - polarizabilities, NMR shielding, EPR g-tensors - which are technically one-particle properties, for which the 2-RDMs do not appear explicitly in the energy derivative expression, both the unperturbed 2-RDMs and the perturbed 2-RDMs are only needed when orbital relaxation needs to be computed.

3.12.4.3. Orbital Relaxation

By including the orbital relaxation effect into property calculations, we essentially allow the MO to respond to external perturbations. This leads to the following parameterization of MO coeffcients:

\[ \mathbf{C}(\chi) = \mathbf{C}(0) \mathbf{U}(\chi) \]

then MO derivatives could be written as:

\[ \mathbf{C}^{\chi}\equiv \frac{\partial \mathbf{C}(\chi)}{\partial \chi} = \mathbf{C}(0) \frac{\partial \mathbf{U}(\chi)}{\partial \chi} \equiv \mathbf{C} \mathbf{U}^{\chi} \]

where \(\mathbf{U}^{\chi}\) are solutions to the CPSCF equations. For unrelaxed calculations the MO coefficients are always \(\mathbf{C}(0)\). In relaxed calculations, the MO response enters relevant CC equations in two ways.

First, as described above, for second-order response properties (like polarizabilities, NMR shielding, etc.) we need to solve the first-order \(T\) and \(\Lambda\) equations. Within a relaxed scheme, any derivative integrals (Fock, ERI, etc.) entering the CC equations contain CPSCF contributions. For example:

\[\begin{split} \frac{\partial (pq|rs)}{\partial \chi} = \\ \sum_{\mu\nu\kappa\tau}C_{\mu p}C_{\nu q} \frac{\partial (\mu\nu|\kappa\tau)}{\partial \chi}C_{\kappa r}C_{\tau s} \\ + \sum_{\mu\nu\kappa\tau}C_{\mu p}^{\chi}C_{\nu q}(\mu\nu|\kappa\tau)C_{\kappa r}C_{\tau s} \\ + \sum_{\mu\nu\kappa\tau}C_{\mu p}C_{\nu q}^{\chi}(\mu\nu|\kappa\tau)C_{\kappa r}C_{\tau s} \\ + \sum_{\mu\nu\kappa\tau}C_{\mu p}C_{\nu q}(\mu\nu|\kappa\tau)C_{\kappa r}^{\chi}C_{\tau s} \\ + \sum_{\mu\nu\kappa\tau}C_{\mu p}C_{\nu q}(\mu\nu|\kappa\tau)C_{\kappa r}C_{\tau s}^{\chi} \end{split}\]

For unrelaxed calculations the last 4 terms are zero.

Second, to compute the relaxed 1- and 2-RDMs and their derivatives, the \(z\)-vector contribution needs to be included. This is a consequence of requiring the CC energy to be invariant to orbital rotation, in the presence of external perturbations. This nice trick by Schaefer and Handy allows gradient to be evaluated without having to solve \(3N_{atoms}\) equations for the orbitals. Instead, only the so-called “\(z\)-vector” equation needs to be solved. It thus opens a door towards efficient analytic response property evaluations. The \(z\)-vector equations are of the same form as the CPSCF equations:

\[ \mathbf{A}\mathbf{z} = \mathbf{X} \]

where the \(\mathbf{A}\) matrix is exactly the same as the CPSCF LHS, and the RHS matrix \(\mathbf{X}\) is the orbital gradient of the CC energy functional. For second-order response properties, the derivatives (subject to each perturbation requested) of this equation also need to be solved.

The discussion of relaxed properties vs unrelaxed ones are quite common. Generally, we believe that for electric-field dependent properties like polarizabilities, inclusion of orbital relaxation does not necessarily lead to more accurate values compared with experiments. However it is generally regarded essential to include orbital relaxation for magnetic-field dependent properties like NMR and EPR parameters.

Finally, it is also worth noting that this orbital relaxtion scheme does not recover the MO response of using HF orbitals in a CC wavefunction setup, which relies on the orbital-optimized coupled-cluster (OOCC) methods, as described in the previous section.

3.12.4.4. Input Options

The input parameters for the wavefunctions are controlled by the %autoci block, while the property-related parameters are controlled by respective property blocks, like elprop and eprnmr. Some useful options are gathered below.

%autoci
    citype
        CCSD(T)
        CCSD
    STol      1e-06 # residue convergence tolerance
    MaxIter   50    # maximum number of iterations
    MaxDIIS   5     # depth of the DIIS memory
    density         # need at least unrelaxed density (see details below)
        unrelaxed
        relaxed
end

%elprop
    polar         true    # dipole-dipole polarizability
    polardipquad  true    # dipole-quadrupole polarizability
    polarquadquad true    # quadrupole-quadrupole polarizability
end

%eprnmr
    gtensor       true    # calculate the g-tensor
    NMRShielding  1       # for chosen nuclei - specified with the Nuclei keyword
                  2       # for all nuclei, equivalent to the 'NMR' simple input
end

Output files: Calculation outputs generated from AutoCI-Response modules are re-directed into separate outputfiles, with the names basename.autociresp_efield, basename.autociresp_bfield, etc.

3.12.4.5. Input Examples

Highly accurate NMR shielding for HF:

! AUTOCI-CCSD(T) NMR pcsseg-4 cc-pwcv5z/c ExtremeSCF

%autoci
  STol 1e-10
  density relaxed
end

%eprnmr
  Tol 1e-12
end

%coords
  CTyp xyz
  Charge 0
  Mult 1
  Units bohrs
  coords
    H  0.00000000  0.00000000  1.64411926
    F  0.00000000  0.00000000 -0.08721704
  end
end

G-Tensor for NH radical, with the SOC operator being the effective nuclear charge (see SOCoptions):

! UHF cc-pVTZ G-Tensor ExtremeSCF AUTOCI-CCSD(T)
! NoFrozenCore Bohrs

%rel 
  SOCType 1
  Zeff[1] = 1.0   # H
  Zeff[7] = 4.55  # N
end

%eprnmr
  Tol 1e-10
end

%autoci
  STol 1e-9
  Density relaxed
end

# NH radical
* xyz 0 3
N 0.0  0.0  0.0
H 0.0  0.0  1.957794
*  

By changing the %rel block to the following we can do the same calculation, using a more accurate SOC description (SOMF):

%rel 
  # these two are defaults
  SOCType 3
  SOCFlags 1,4,3,0
end

3.12.5. Fully Internally Contracted MRCI (FIC-MRCI)

Starting point for any multireference approach is a reference wavefunction that consists of multiple determinants or configurations state functions (CSFs). In many instances this is the complete active space SCF (CASSCF) wavefunction. In the uncontracted MRCI approach, as implemented in the orca_mrci module, the wavefunction is expanded in terms of excited CSFs that are generated by considering excitations with respect to all reference CSFs. The methodology scales with the number of reference CSFs and hence is restricted to small reference spaces. Moreover, the configuration driven algorithm used in orca_mrci keeps all integrals in memory, which further limits the overall size of the molecule.

Internal contraction as proposed by Meyer and Siegbahn avoids these bottlenecks [416, 417]. Here, excited CSFs are generated by applying the excitation operator to the reference wavefunction as whole. The fully internally contracted MRCI (FIC-MRCI) uses the same internal contraction scheme as the FIC-NEVPT2 (aka PC-NEVPT2). The entire methodology as well as a comparison with the conventional uncontracted MRCI is reported in our article [418]. The CEPA0, ACPF and AQCC variants are straight forward adoptions [419]. The residue of the FIC-MRCI ansatz

(3.130)\[ R_K=\braket{\Phi^{pr}_{qs}| \hat H - E^\text{CAS}-\lambda E_{c}|\Psi^\text{FIC} }, \]

is modified by the factor

\[\begin{split}\lambda = \left\{\begin{array}{cc} 1 & \text{MRCI} \\ 0 & \text{CEPA0} \\ \frac{2}{N_{e} } & \text{ACPF}\\ \frac{1-(N_e-3)(N_e-2) }{N_e\cdot(N_e-1) } & \text{AQCC} \end{array} \right.\end{split}\]

Here, \(E_{c}\) is the correlation energy and \(\Phi^{pr}_{qs}\) denote the internally contracted CSF that arise from the action of the spin-traced excitation operators on the CAS-CI reference wave function

\[\Phi^{pr}_{qs}=E^p_qE^r_s \ket{\Psi^\text{CAS} }.\]

In case of ACPF and AQCC, the \(\lambda\) factor explicitly depends on the number of correlated electrons (\(N_{e}\)).

The general input structure is like that of the CASSCF module, e.g., the following example input reads an arbitrary set of orbitals and starts the FIC-MRCI calculation. The internal contracted formalism requires CAS-CI reduced densities up to fourth-order, which can be expensive to construct. By default, the density construction is sped up using the prescreening (PS) approximation reported in Section N-Electron Valence State Perturbation Theory (NEVPT2).[420]

!def2-tzvp moread allowrhf noiter nofrozencore
%moinp "start.gbw"  # could be from CASSCF

%autoci
citype   FIC-MRCI      # Fully internally contracted MRCI (singles, doubles)
         FIC-MRCEPA(0) # CEPA0 version of FIC-MRCI
         FIC-MRACPF    # ACPF  version of FIC-MRCI
         FIC-MRAQCC    # AQCC  version of FIC-MRCI
         FIC-DDCI3     # FIC-MRCI without the IJAB excitation

# CAS-CI reference wavefunction
nel 2
norb 2
mult 1,3
nroots 3,1

nthresh 1e-6  # removal of linear dependencies in the IC-CSFs
D3TPre        1e-14 # default density truncation of the 3-RDM
D4TPre        1e-14 # default density truncation of the 4-RDM

# Davidson correction for the FIC-MRCI
DavidsonOpt 0  # none (default)
            1  # Davidson correction
end

Currently, the program is capable of computing total energies and vertical excitation energies. More features will be available with future releases.

3.12.5.1. Input Example

# potentially "!FIC-MRCEPA(0)" omitting the %autoci block.
!def2-SVP  

# get CASSCF orbitals and define states
%casscf
  nel 8
  norb 6
  mult 3
  nroots 1
end

# run FIC-MRCEPA(0)
%autoci
  # CAS settings are automatically copied from the %casscf block, if 
  # not defined in %autoci
 
  citype FIC-MRCEPA(0)
end

# O2 molecule
*xyz 0 3
O 0 0 0 
O 0 0 1.28
*

3.12.6. Fully Internally Contracted MRCC (FIC-MRCC)

Several approaches have been taken towards extending CC theory to work with genuinely multiconfigurational reference wave functions [421], yet none of these approaches has found widespread adoption. As of 2011, the internally contracted MRCC theory has had a revival, with a rigorous theoretical investigation of several approximations that also proved its orbital invariance [422] and a first report of a polynomial-scaling code obtained through automatic equation generation [423].

Our implementation in ORCA is akin to the previously published formulations in Refs. [422, 423], although everything is formulated rigorously in terms of the spin-free excitation operators \(\hat E_q^p = \hat a^{p\alpha }\hat a_{q\alpha } + \hat a^{p\beta }\hat a_{q\beta }\), using an improved version of the ORCA-AGE code generator.[360] To begin with, the ansatz for the wave function is

(3.131)\[\begin{aligned} |\Psi \rangle = \text{e}^{\hat T}|0\rangle \;, \end{aligned}\]

where \(|0\rangle\) denotes a zeroth order CASSCF reference wave function and the cluster operator can be written as (Einstein’s summation convention implied)

\[\begin{split} \hat T = \frac{1}{2}t_{ab}^{ij}\hat E_i^a\hat E_j^b + t_{ab}^{it}\hat E_i^a\hat E_t^b + \frac{1}{2}t_{ab}^{tu}\hat E_t^a\hat E_u^b + t_{at}^{ij}\hat E_i^a\hat E_j^t + t_{au}^{it}\hat E_i^a\hat E_t^u \\ + t_{ua}^{it}\hat E_i^u\hat E_t^a + t_{av}^{tu}\hat E_t^a\hat E_u^v + \frac{1}{2} t_{tu}^{ij}\hat E_i^t\hat E_j^u + t_{uv}^{it}\hat E_i^u\hat E_t^v\;. \end{split}\]

Note that we do not use normal order in the cluster operator.

Inserting the ansatz from Eq. (3.131) into the Schrödinger equation and pre-multiplying with the inverse exponential, we obtain the similarity transformed Hamiltonian and the energy expression,

\[\begin{aligned} \bar H|0\rangle = \text{e}^{ - \hat T}\hat H\text{e}^{\hat T}|0\rangle = E|0\rangle \;.\end{aligned}\]

In our code, the similarity-transformed Hamiltonian is truncated after the quadratic terms since that approximation has been found to only have minor impact on the accuracy of the method [422],

\[\begin{aligned} \bar H \approx \hat H + [\hat H,\hat T] + \frac{1}{2}[[\hat H,\hat T],\hat T]\;.\end{aligned}\]

The residual conditions are subsequently obtained by projecting with contravariant excited functions \(\langle \tilde \Phi _P|\) onto the Schrödinger equation,

\[\begin{aligned} { r_P} = \langle \tilde \Phi _P|\bar H|0\rangle \;.\end{aligned}\]

For a definition of the contravariant projection functions, we refer to Ref. [418] since this FIC-MRCC implementation uses the same contravariant functions as the published FIC-MRCI implementation. Despite using contravariant projection functions, this is not sufficient to remove all linear dependencies from the set of projection functions \(\{ \tilde\Phi_P \}\), i.e. the metric matrix

\[\begin{aligned} S_{PQ} &= \braket{\tilde\Phi_P | \tilde\Phi_Q} \ne \delta_{PQ}\end{aligned}\]

has off-diagonal elements within excitation classes and between classes with the same number of inactive and virtual indices (ITAU and ITUA). Hence, the set of projection functions needs to be orthonormalized, which is achieved with Löwdin’s canonic orthogonalization in the AUTOCI module.[1]

3.12.6.1. Input Example

The FIC-MRCC module can be started by specifying the CIType keyword in the %autoci block or by adding FIC-MRCC to the simple input line of an ORCA input file. The following example computes the singlet ground state energy of four hydrogen atoms arranged as a square with a side length of \(2 a_0\), which is commonly known as the H4 model [424].

! cc-pVTZ Bohrs  # it is possible to add the `FIC-MRCC' keyword here
                 # and omit the %autoci block below

%maxcore 10000

# get CASSCF orbitals and define states.
%casscf
  nel 2
  norb 2
  mult 1
  nroots 1
end

# run FIC-MRCC
%autoci  
  # CAS settings are automatically copied from the %casscf block, if 
  # not defined in %autoci
  citype fic-mrcc
  
  # optional parameters
  nthresh 1e-6        # default removal of linear dependencies in the IC-CSFs
  D3TPre        1e-14 # default density truncation of the 3-RDM
  D4TPre        1e-14 # default density truncation of the 4-RDM
  D5TPre        1e-14 # default density truncation of the 5-RDM
  Density  linearized  # compute the linearized density (not default).
end

* int 0 1
  H 0 0 0 0.0  0.0 0.0
  H 1 0 0 2.0  0.0 0.0
  H 2 1 0 2.0 90.0 0.0
  H 1 2 3 2.0 90.0 0.0
*

In this example, ORCA will first run a state-specific CASSCF calculation, and then immediately continue with the fic-MRCC calculation on top of the CASSCF solution from the first step. It is, however, not required to always run a CASSCF calculation before the autoci module. Any ORCA gbw/mp2nat/… file is accepted through %moinp, although that route requires the user to specify the active space in the autoci block. autoci will then compute a CASCI solution with the provided input orbitals and use that information to drive the correlated calculations.

Please be aware that FIC-MRCC is a very extensive theory, which leads to long run times. The computational effort depends mainly on the number of orbitals, the number of total electrons and the size of the active space. On modestly modern hardware, calculations of \(\sim 300\) orbitals with a CAS(2,2) should be readily achievable. For larger active spaces, such as a CAS(6,6), calculations with a total of \(\sim 200\) orbitals will also complete within a day.

3.12.7. Fully Internally Contracted NEVPT (FIC-NEVPT4(SD)/FIC-NEVPT3/FIC-NEVPT2)

Perturbation theory is a cost-effective alternative to the previously mentioned FIC-MRCI and FIC-MRCC methods. Most often, a reasonable description of electron correlation is already achieved at the second order of the expansion. The CASSCF module features an efficient implementation of the second-order methods, that are able to tackle large molecules and active spaces. For more details and references on the second-order variants, we refer to the dedicated section N-Electron Valence State Perturbation Theory (NEVPT2).
In the following section, the focus is on higher orders of the n-electron valence state perturbation theory (NEVPT), that are available within the AUTOCI framework. Here, all integrals are computed and kept on disk, which limits the application to approximately 700 atomic orbitals (AOs). The appearance of higher-order reduced density matrices is a common feature of internally contracted methods. Like FIC-NEVPT2, the FIC-NEVPT3 and FIC-NEVPT4(SD) require at most the fourth-order reduced density matrix. Its computation can be accelerated with the prescreening method discussed in section N-Electron Valence State Perturbation Theory (NEVPT2).[420, 425] In terms of computation time, the FIC-NEVPT4(SD) is as expensive as FIC-NEVPT3 or a single iteration of the FIC-MRCI method.

3.12.7.1. FIC-NEVPT3

The FIC-NEVPT3 methodlogy is desribed in detail by Angeli et al.[426] With the first order wave function \(\Psi^{(1)}\), the third order energy is given by the expression

\[E^{(3)}=\langle \Psi^{(1)}|V|\Psi^{(1)} \rangle.\]

3.12.7.2. FIC-NEVPT4(SD)

The fourth-order perturbation theory requires the second wave function, which, in principle, involves triple and quadruple excitation defining the second-order interacting space. The inclusion of the latter is computationally demanding. In the following we report an approximate solution, coined FIC-NEVPT4(SD), which is restricted to the first-order interacting space - the single and double excitations with respect to the reference wave function. This method is reported and analyzed in more detail by Kempfer et al.[427] Denoting the second order wave function with \(\Psi^{(2)}\), the fourth order energy is readily expressed as

\[E^{(4)} = \langle\Psi^{(1)}|V|\Psi^{(2)}\rangle - E^{(2)}\langle\Psi^{(1)}|\Psi^{(1)}\rangle.\]

The second contribution in the expression above, \(E^{(2)}\langle\Psi^{(1)}|\Psi^{(1)}\rangle\), is a renormalization term. As discussed in the paper,[427] the renormalization term is dropped (default) as it destroys the size consistancy in the absence of the exact quadruple excitations.

3.12.7.3. Input Example

# potentially "!FIC-NEVPT4(SD)"/"!FIC-NEVPT3" instead of the %autoci block.
!def2-tzvp 

# get CASSCF orbitals and define the states
%casscf
    nel 8
    norb 6
    mult 3
    nroots 1
end


# run NEVPT4(SD)
%autoci
# CAS settings are automatically copied from the %casscf block, if 
# not defined in %autoci

 citype fic-nevpt4(sd)

# optional parameters
nthresh 1e-6        # default removal of linear dependencies in the IC-CSFs
D3TPre        1e-14 # default density truncation of the 3-RDM
D4TPre        1e-14 # default density truncation of the 4-RDM

Density  Unrelaxed  # compute the unrelaxed density (not default).
end

# O2 molecule
* xyz 0 3
O 0 0 0
O 0 0 1.28
*

3.12.8. FIC Ansatz: Unrelaxed Densities and Natural Orbitals

The FIC-NEVPT4(SD), the FIC-MRCC and the FIC-MRCI variants have a common set excitation parameters in their wavefunction. For these methods, (linearized) unrelaxed densities are available (density unrelaxed). With !KeepDens the densities are stored in the density container (.densities file on disk), which can be read with the orca_plot program. The name of the densities reflect its multiplicity, root count and the actual method. Below is an example snippet of the “orca_plot jobname.densities” call:

---------------------
List of density names
---------------------

Index:                                                   Name of Density
------------------------------------------------------------------------
    0:                                                      jobname.scfp
    1:                                                      jobname.scfr
    2:                                         Tdens-CAS.mult.3.root.0.p
    3:                                         Tdens-CAS.mult.3.root.1.p
    4:                                         Tdens-CAS.mult.3.root.2.p
    5:                                         Tdens-CAS.mult.3.root.0.r
    6:                                         Tdens-CAS.mult.3.root.1.r
    7:                                         Tdens-CAS.mult.3.root.2.r
    8:                     jobname.mult.3.root.0.FIC-NEVPT4(SD).autoci.p
    9:                     jobname.mult.3.root.1.FIC-NEVPT4(SD).autoci.p
   10:                     jobname.mult.3.root.2.FIC-NEVPT4(SD).autoci.p

Here, the densities 8-10 are unrelaxed densities of the FIC-NEVPT4(SD) method for the first three triplet roots. With the unrelaxed densities, natural orbitals (natorbs true) can be generated. The resulting .gbw files ending with (.nat) can be used to carry out further analysis, e.g. using the MultiWFN program.

!FIC-NEVPT4(sd)
...
%autoci
  Density  Unrelaxed  # compute the unrelaxed density
  natorbs true        # generate natural orbitals
end

3.12.9. Keywords

Table 3.44 Simple input keywords for AUTOCI methods.

Keyword	Description
AUTOCI-<method>	Activate single-reference correlation module from AUTOCI. Disables MDCI.
FIC-MRCI
FIC-MRCEPA(0)
FIC-MRACPF
FIC-MRAQCC
FIC-DDCI3
FIC-MRCC
FIC-NEVPT2
FIC-NEVPT3
FIC-NEVPT4(SD)

Table 3.45 %autoci block input keywords AUTOCI methods.

Keyword	Options	Description
CIType	CID	Correlation method
	CISD
	CISDT
	CCD
	CCSD
	CCSDT
	CEPA(0)
	QCISD
	CC2
	CC3
	CCSD[T]
	CCSD(T)
	CCSDT-1
	CCSDT-2
	CCSDT-3
	CCSDT-4
	FIC-MRCI
	FIC-MRCEPA(0)
	FIC-MRACPF
	FIC-MRAQCC
	FIC-DDCI3
	FIC-MRCC
	FIC-NEVPT2
	FIC-NEVPT3
	FIC-NEVPT4(SD)
	MP2
	MP3
	MP4(SDQ)
	MP4
	MP5
STol	1e-6	Residue convergence tolerance
MaxIter	50	Maximum number of iterations
MaxDiis	5	Depth of the DIIS memory
DIISStartIter	2	Apply DIIS starting at iteration 1 (counting starts with 0).
ExcludeHigherExcDIIS	false	Exclude triples and higher excitations from DIIS procedure
NEl	6	Number of active electrons (CAS only)
NOrb	7	Number of active orbitals (CAS only)
Mult	1	Requested multiplicity block (CAS only)
NRoots	1	Number of roots for mult block (CAS only)
Irrep	0	Requested irrep for mult block (CAS only)
NThresh	1e-6	Threshold for lin. dependencies in the IC-CSFs basis (CAS only)
D3TPre	1e-14	Density truncation in D3 (CAS only)
D4TPre	1e-14	Density truncation in D4 (CAS only)
PrintLevel	3	Amount of printing (1-7)
TrafoType	NoRI	Type of integral transformation
	RITrafo
KeepInts	false	Keep the transformed integrals on disk
UseOldInts	false	Use the transformed integrals found on disk
RunROHFasUHF	false	Invokes AUTOCI UHF modules with orbitals from ROHF calculation
Density	none	Type of 1-body density requested
	linearized
	unrelaxed
	relaxed
Density2	none	Type of 2-body density requested
	linearized
	unrelaxed
	relaxed
NatOrbs	false	Generate natural orbitals (not available for all methods)

N.B. For a ROHF reference, only CISD calculations can be performed in the current version. However, it is possible to run UHF calculations with an ROHF reference by setting the RunROHFasUHF flag to true. Note that this only makes sense when the reference is indeed ROHF, e.g. (calculating the isotropic part of the HFC at CCSD level, running AutoCI UHF CCSD module, with orbitals obtained from the ROHF SCF calculation):

! ROHF def2-svp tightscf pmodel AUTOCI-CCSD

%autoci
    RunROHFasUHF true
end

* xyz 0 2
Cu  0.0  0.0  0.0
*

%eprnmr nuclei = all Cu {aiso} end

If one wishes to experiment with the module itself and the reference wavefunction stays constant, it is possible to store the transformed MO integrals on disk (keepints) and reuse them (useoldints). The program checks only whether the dimension of the integrals on disk match the problem actually solved, i.e. the user is responsible for valid data.

---

[1]

This is similar to scheme A from Ref. [423].