7.25. CI methods using generated code

The AUTOCI module is a replacement of the orca_mdci for cases, where manual implementation of the method 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.

7.25.1. Introduction

All the theories are obtained by the means of automated programming within the ORCA-AGE (Automated Generator Environment for ORCA).[474, 500] 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.

7.25.2. Input

All parameters applicable to the AUTOCI module are shown below.

%autoci
    # Algorithm selection
    citype      # Type of the CI expansion to be applied (one of following)
        CID     # Configuration Interaction with double substitutions
        CISD    # Configuration Interaction with singles and doubles
        CISDT   # Configuration Interaction with singles, doubles and triples
        CCD     # Coupled Cluster with double substitutions
        CCSD    # Coupled Cluster with single and double substitutions
        CCSDT   # Coupled Cluster with single, double and triple substitutions
        CEPA(0) # Linearized CCSD
        QCISD   # Quadratic CISD
        CC2     # Approximate CCSD
        CC3     # Approximate CCSDT
        CCSD[T]    # CCSD with perturbative [T] correction
        CCSD(T)    # CCSD with perturbative (T) correction
        CCSDT-1    # Approximate CCSDT, CCSDT-1
        CCSDT-2    # Approximate CCSDT, CCSDT-2
        CCSDT-3    # Approximate CCSDT, CCSDT-3
        CCSDT-4    # Approximate CCSDT, CCSDT-4
        FIC-MRCI   # Fully internally contracted MRCI
        FIC-MRCEPA(0)   # Fully internally contracted CEPA(0)
        FIC-MRACPF      # Fully internally contracted ACPF
        FIC-MRAQCC      # Fully internally contracted AQCC
        FIC-DDCI3       # FIC-MRCI without the IJAB excitation
        FIC-MRCC        # Fully internally contracted MRCC
        MP2             # Second order Moeller-Plesset perturbation theory
        MP3             # Third order Moeller-Plesset perturbation theory
        MP4(SDQ)        # MP4 without triple substitutions	
        MP4             # Fourth order Moeller-Plesset perturbation theory	
        MP5             # Fifth order Moeller-Plesset perturbation theory

    # converger details
    stol       1e-06 # residue convergence tolerance
    maxiter    50    # maximum number of iterations
    maxdiis    5     # depth of the DIIS memory
    diisstartiter 2  # Apply DIIS starting at iteration 1
    ExcludeHigherExcDIIS false # exclude triples and higher excitations from DIIS procedure

    # CAS settings similar to the CASSCF input
    nel        6         # number of active electrons (for CAS)
    norb       7         # number of active orbitals (for CAS)
    mult           1 # requested multiplicity block
    nroots         1 # number of roots for mult block
    irrep          0 # requested irrep for mult block
    nthresh       1e-6 # Threshold for lin. dependencies in the IC-CSFs basis
    D3TPre        1e-14 # Density truncation in D3
    D4TPre        1e-14 # Density truncation in D4

    # Algorithm details
    printlevel   3    # Amount of printing
    trafotype  0      # Type of integral transformation
        0             # Full canonical
        1             # Full using RI (RI basis needed)

    keepints		# Keep the transformed integrals on disk
    useoldints		# Use the transformed integrals found on disk
    RunROHFasUHF    # Invokes AUTOCI UHF modules with orbitals from ROHF calculation

    # Property calculations
    density           # type of density requested
         none         # No density calculation
         linearized   # Linear part of the density
         unrelaxed    # Unrelaxed 1-body density matrix
         relaxed      # Relaxed 1-body density matrix

end

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.

7.25.3. Available Properties

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	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 method can be called from the simple keyword by prepending AUTOCI- to the correlation method, for instance

! AUTOCI-CCSD

7.25.4. 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).[500]

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.

7.25.5. AUTOCI Response Properties via Analytic Derivatives

For single-reference methods (currently limited to CCSD and MP2), some response properties could be calculated via taking the analytic derivative of the wavefunctions computed by AUTOCI.

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 shown below.

%autoci
    citype          # only CCSD and MP2 available
        CCSD
        MP2
    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    # polarizability via analytic derivative
end

%eprnmr
    NMRShielding 2  # NMR shielding for all nuclei, equivalent to the 'NMR' simple input
end

N.B. For the response property calculations, the electron density and electron response density needs to be calculated. Currently for the analytic polarizability at CCSD level, only the unrelaxed density version is implemented:

! UHF AUTOCI-CCSD Def2-SVP NoFrozenCore

%elprop
    polar true
end

* xyz 0 1
O	 0.0000000000  0.0000000000 -0.1190150726
H	 0.7685504811  0.0000000000  0.4760602904
H	-0.7685504811  0.0000000000  0.4760602904
*

For MP2, both the polarizability and NMR shielding needs the relaxed densities:

! UHF AUTOCI-MP2 Def2-SVP NoFrozenCore

%elprop
    polar true
end

%autoci
    density relaxed
end

* xyz 0 1
O	 0.0000000000  0.0000000000 -0.1190150726
H	 0.7685504811  0.0000000000  0.4760602904
H	-0.7685504811  0.0000000000  0.4760602904
*

! UHF AUTOCI-MP2 Def2-SVP NMR NoFrozenCore

%autoci
    density relaxed
end

* xyz 0 1
O	 0.0000000000  0.0000000000 -0.1190150726
H	 0.7685504811  0.0000000000  0.4760602904
H	-0.7685504811  0.0000000000  0.4760602904
*

Please also note that for AUTOCI NMR calculations, the GIAOs (gauge-including atomic orbitals) are necessary (turned on by default). Also, there is no fronzen core options implemented yet (NoFrozenCore keyword is needed).

7.25.6. Fully Internally Contracted 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 [584, 796]. Here, excited CSFs are generated by applying the excitation operator to the reference wavefunction as whole. The fully internally contracted MRCI presented here (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 [805]. The CEPA0, ACPF and AQCC variants are straight forward adoptions [741]. The residue of the FIC-MRCI ansatz

(7.188)\[ 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 speed up using the prescreening (PS) approximation reported in Section N-Electron Valence State Pertubation Theory.

!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-10 # 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.

7.25.7. Fully Internally Contracted MRCC

Several approaches have been taken towards extending CC theory to work with genuinely multiconfigurational reference wave functions [538], 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 [249] and a first report of a polynomial-scaling code obtained through automatic equation generation [357].

Our implementation in ORCA is akin to the previously published formulations in Refs. [249, 357], 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.[500] To begin with, the ansatz for the wave function is

(7.189)\[\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{gathered} \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{gathered}\)

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

Inserting the ansatz from Eq. (7.189) 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 [249],

\[\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. [805] 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]

7.25.7.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 [418].

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

%maxcore 10000

%casscf
  nel 2
  norb 2
  mult 1
  nroots 1
end

%autoci  # CAS settings are automatically copied from the CASSCF block!
  citype fic-mrcc
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.

Note that in the current implementation, autoci_fic_mrcc requires large stack sizes, which may lead to segmentation faults immediately after starting the program. We recommend a stack size of at least 16 MiB. On Unix-like systems, this can be achieved by issuing the command ulimit -s 16384 in the terminal or in the submission script prior to running the fic-MRCC module.

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.

---

[1]

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