(sec:autoci.detailed)= # 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. (sec:autoci.introduction.detailed)= ## Introduction All the theories are obtained by the means of automated programming within the ORCA-AGE (Automated Generator Environment for ORCA).{cite}`krupicka2017ORCA_AGE,lechner2024code` 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. ## 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. (sec:autoci.available.detailed)= ## 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 ```orca ! AUTOCI-CCSD ``` (sec:autoci.gradients.detailed)= ## 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).{cite}`lechner2024code` An example is shown below how to optimise a geometry using AUTOCI's gradients at the CCSD level of theory ```orca ! 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 ```orca ! 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 ```orca %autoci density relaxed end ``` For geometry optimisations, both the unrelaxed and relaxed densities are computed automatically and do not need to be requested explicitly. (sec:autoci.autociresp.detailed)= ## 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: ```orca ! 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: ```orca ! 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 * ``` ```orca ! 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). (sec:autoci.ficmrci.detailed)= ## 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 {cite}`siegbahn_direct_1980,meyer_configuration_1977`. 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 {cite}`sivalingam_comparison_2016`. The CEPA0, ACPF and AQCC variants are straight forward adoptions {cite}`saitowfic`. The residue of the FIC-MRCI ansatz $$ R_K=\braket{\Phi^{pr}_{qs}| \hat H - E^\text{CAS}-\lambda E_{c}|\Psi^\text{FIC} }, $$ (residualfic) is modified by the factor $$\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.$$ 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 {ref}`sec:nevpt2.detailed`. ```orca !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. (sec:autoci.ficmrcc.detailed)= ## Fully Internally Contracted MRCC Several approaches have been taken towards extending CC theory to work with genuinely multiconfigurational reference wave functions {cite}`Lyakh2012`, 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 {cite}`Evangelista2011` and a first report of a polynomial-scaling code obtained through automatic equation generation {cite}`Hanauer2011`. Our implementation in ORCA is akin to the previously published formulations in Refs. {cite}`Evangelista2011,Hanauer2011`, 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.{cite}`lechner2024code` To begin with, the *ansatz* for the wave function is $$\begin{aligned} |\Psi \rangle = \text{e}^{\hat T}|0\rangle \;, \end{aligned}$$ (eq:MRCC-ansatz) 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. {eq}`eq:MRCC-ansatz` 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 {cite}`Evangelista2011`, $$\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. {cite}`sivalingam_comparison_2016` 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] (sec:autoci.ficmrcc.example.detailed)= ### 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 {cite}`Jankowski1980`. ```orca ! 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. {cite}`Hanauer2011`.