(sec:casscf.detailed)= # The Complete Active Space Self-Consistent Field (CASSCF) Module (sec:casscf.general.detailed)= ## General Description The complete active space self-consistent field (CASSCF) method is a special form of a multiconfigurational SCF method and can be thought of as an extension of the Hartree-Fock method. It is a very powerful method to study static correlation effects and a solid basis for MR-CI and MR-PT treatments. It can be applied to the ground state and excited states or averages thereof. The implementation in ORCA is fairly general and reasonably efficient. However, CASSCF calculations are fairly complex and ultimately require a lot of insight from the user in order to be successful. In addition to detailed description here, the manual explores some typical examples in section {ref}`sec:energygradients.casscf.natorbitals.typical`. Furthermore, the manual is supplemented with a tutorial for CASSCF that covers many practical tips on the calculation design and usage of the program. **The wavefunction.** The wavefunction of a given CASSCF state is written as $$\left| \Psi_I^S \right\rangle= \sum_{k} { C_{kI} \left| \Phi_k^S \right\rangle}. $$ (eqn:148) Here, $\left| \Psi_I^S \right\rangle$ is the CASSCF $N$-electron wavefunction for state $I$ with total spin S. The set of $\left| \Phi_k^S \right\rangle$ is a set of configuration state functions (for example linear combination of Slater determinants) each adapted to a total spin $S$. The expansion coefficients $C_{kI}$ represent the first set of variational parameters. Each CSF is constructed from a common set of orthonormal molecular orbitals $\psi_{i} \left({ \mathrm{\mathbf{r} }} \right)$ which are in turn expanded in basis functions $\psi_{i} \left({ \mathrm{\mathbf{r} }} \right)=\sum\nolimits_\mu{c_{\mu i} \phi_{\mu } \left({ \mathrm{\mathbf{r} }} \right)}$. The MO coefficients $c_{\mu i}$ form the second set of variational parameters. **The energy.** The energy of the CASSCF wavefunction is given by the Rayleigh quotient $$E\left({ \mathrm{\mathbf{c} },\mathrm{\mathbf{C} }} \right)=\frac{\left\langle { \Psi _{I}^{S} \left|{ \hat{{H} }_{\text{BO} } } \right|\Psi_{I}^{S} } \right\rangle}{\left\langle { \Psi_{I}^{S} \left|{ \Psi_{I}^{S} } \right.} \right\rangle}, $$ (eqn:149) and represents an upper bound to the true total energy. However, CASSCF calculations are *not* designed to provide values for total energy which are close to the exact energy. The purpose of a CASSCF calculation is to provide a qualitatively correct wavefunction, which forms a good starting point for a treatment of dynamic electron correlation. The CASSCF method is fully variational in the sense that the energy is made stationary with respect to variations in both sets of MO and CI coefficients. At convergence, the gradient of the energy with respect to the MO and CI coefficients vanishes $$\frac{\partial E\left({ \mathrm{\mathbf{c} },\mathrm{\mathbf{C} }} \right)}{\partial c_{\mu i} }=0, $$ (eqn:150) $$\frac{\partial E\left({ \mathrm{\mathbf{c} },\mathrm{\mathbf{C} }} \right)}{\partial C_{kI} }=0. $$ (eqn:151) **Orbital spaces.** In CASSCF calculations, the MO space is divided into three user defined subspaces: - The "inactive orbitals" are the orbitals which are doubly occupied in all configuration state functions (labels $i, j, k, l$). - The "active orbitals" are the orbitals with variable occupation numbers in the various CSFs (labels $t, u, v, w$). - The "external orbitals" (labels $a, b, c, d$) Note that in older publications, the inactive and active orbitals are distinguished and referred to as "internal" orbitals. The wavefunction and energy is invariant with respect to unitary transformations within the three subspaces. The special feature of a CASSCF wavefunction is that a fixed number of electrons is assigned to each subspace. The internal subspace is of course completely filled but the CSFs in the active space constitute a full-CI of $n$-electrons in $m$-orbitals. The CSF list is constructed such, however, that a wavefunction of well defined total spin (and potential space) symmetry results. Such a wavefunction is referred to as a **CASSCF($n$,$m$)** wavefunction. The CSF list grows extremely quickly with the number of active orbitals and the number of active electrons (basically factorially). Depending on the system, the limit of feasibility is roughly around $\sim$14 active orbitals or about one million CSFs in the active space. Larger active spaces are tractable with approximate CI solver such as the Iterative-Configuration-Expansion CI (ICE-CI) described in {ref}`sec:iceci.detailed` or the Density Matrix Renormalization Group (DMRG) discussed in {ref}`sec:dmrg.detailed`.[^1] Since the orbitals within the subspaces are only defined up to a unitary transformation, the program needs to make some canonicalization choice. In ORCA, the final orbitals by default are: 1. natural orbitals in the active space, 2. orbitals which diagonalize the CASSCF Fock matrix in the internal space and 3. orbitals which diagonalize the CASSCF Fock matrix in the external space. **State averaging.** In many circumstances, it is desirable to optimize the orbitals not for a single state but for the average of several states. In order to see what is done, the energy for state $I$ is re-written as: $$E_{I} \left({ \mathrm{\mathbf{c} },\mathrm{\mathbf{C} }} \right)=\sum\limits_{pq} { \Gamma _{q}^{p\left( I \right)} h_{pq} } +\sum\limits_{pqrs} { \Gamma _{qs}^{pr\left( I \right)} \left({ pq\left|{ rs} \right.} \right)} $$ (eqn:152) Here, $\Gamma_{q}^{p\left( I \right)}$ and $\Gamma_{qs}^{pr\left( I \right)}$ are the one-and two-particle reduced electron density matrices for this state (labels $p, q, r, s$ span the internal and active subspaces): $$\Gamma_{q}^{p\left( I \right)} =\left\langle { \Psi_{I}^{S} \left|{E_{q}^{p} } \right|\Psi_{I}^{S} } \right\rangle$$ (eqn:153) $$\Gamma_{qs}^{pr\left( I \right)} =\frac{1}{2}\left\langle { \Psi_{I}^{S} \left|{ E_{q}^{p} E_{s}^{r} -\delta_{qr} E_{s}^{p} } \right|\Psi_{I}^{S} } \right\rangle$$ (eqn:154) The average energy is simply obtained from averaging the density matrices using arbitrary weights $w_{I}$ that are user defined but are constrained to sum to unity. $$\Gamma_{q}^{p\left({ av} \right)} =\sum\limits_I { w_{I} \Gamma _{q}^{p\left( I \right)} } $$ (eqn:155) $$\Gamma_{qs}^{pr\left({ av} \right)} =\sum\limits_I { w_{I} \Gamma_{qs}^{pr\left( I \right)} } $$ (eq:156) $$\sum\limits_I { w_{I} } =1 $$ (eqn:157) **Optimization of CASSCF wavefunctions.** In general, except for trivial cases, CASSCF wavefunctions are considerably more difficult to optimize than RHF (or UHF) wavefunctions. The underlying reason is that variations in **c** and **C** maybe strongly coupled and the energy functional may have many local minima in (**c**,**C**) space. Consequently, the choice of starting orbitals is of really high importance and the choice which orbitals and electrons are included in the active space has decisive influence on the success of a CASSCF study. In general, after transformation to natural orbitals, one can classify the active space orbitals by their occupation numbers which vary between 0.0 and 2.0. In general, convergence problems are almost guaranteed if orbitals with occupation numbers close to zero or close to 2.0 are included in the active space. Occupation numbers between 0.02 and 1.98 are typically very reasonable and should not lead to large convergence problems. The reason for the occurrence of convergence problems is that the energy is only very weakly dependent on rotations between internal and active orbitals if the active orbital is almost doubly occupied and similarly for the rotations between external and weakly occupied active orbitals. However, in some cases (for example in the study of potential energy surfaces) it may not be avoidable to include weakly or almost inactive orbitals in the active space and in these cases the use of the most powerful convergence aids is necessary (*vide infra*). As in the case of single-determinant wavefunctions (RHF, UHF, RKS, UKS) there are first and second order converging methods available. The first order CASSCF methods require the transformed integrals $(tu|vx)$ with $x$ belonging to any subspace. This is a very small subspace of the total transformed integral list and is readily held in central storage even for larger calculations. On the other hand, second order CASSCF methods require the integrals $(pq|xy)$ and $(px|qy)$ ($p,q=$ internal, active; $x,y=$ any orbital). This is a fairly large set of integrals and their generation is laborious in terms of CPU time and disk storage. Second order CASSCF calculations are therefore more limited in the size of the molecules which can be well treated. It would be possible to basically avoid the integral transformation also in the case of second-order CASSCF calculations and proceed to fully direct calculations. Such calculations may become quite time consuming since there may be a large number of Fock matrix builds necessary. The augmented Hessian method (Newton-Raphson) solves the eigenvalue problem: $$\begin{pmatrix} 0 & \mathbf{g} \\ \mathbf{g} & \mathbf{H} \end{pmatrix} \begin{pmatrix} 1 \\ \mathbf{t} \end{pmatrix} = \varepsilon \begin{pmatrix} 1 \\ \mathbf{t} \end{pmatrix} $$ (eqn:158) Here, $\mathrm{\mathbf{g} }$ is the orbital gradient (derivative of the total energy with respect to a non-redundant rotation between two orbitals) and $\mathrm{\mathbf{H} }$ is the orbital Hessian (second derivative of the energy with respect to two non-redundant orbital rotations). The vector $\mathrm{\mathbf{t} }$ (in intermediate normalization obtained from the CI like vector) summarizes the rotation angles. The angles are used to define the antisymmetric matrix ($X_{pq} =-X_{qp}$ is thus the rotation angle between orbitals $p$ and $q$): $$\mathrm{\mathbf{X} }= \begin{pmatrix} \mathbf{0} & \mathbf{t}\\ \mathbf{-t} & \mathbf{0} \end{pmatrix}, $$ (eqn:159) which is used to parametrize the unitary matrix $\mathrm{\mathbf{U} }=\exp \left({\mathrm{\mathbf{X} }} \right)$ which is used to update the orbitals according to: $$\mathrm{\mathbf{c} }^{\text{new} }=\mathrm{\mathbf{c} }^{\text{old} }\mathrm{\mathbf{U} } $$ (eqn:160) (where $\mathrm{\mathbf{c} }$ is an MO coefficient matrix). **Starting orbitals.** You cannot be careful enough with your starting orbitals. What type of initial guess works best depends on the system. Quite often it is not the magnitude of the initial gradient, but the similarity between initial and final active orbitals. The CASSCF tutorial discusses a number of guess options in more detail. Generally speaking, canonical orbitals HF orbitals from a RHF calculation are not good choice, as the identification and selection of the active space orbitals is often difficult. Usually DFT orbitals (quasi-restricted or RKS) perform better in this respect. Alternatively, if CASSCF orbitals from a previous run or a close-by geometry are available this is a good choice. For coordination chemistry complexes, the guess generated with `orca_mergefrag` (see the CASSCF tutorial), is probably the best choice - especially for heave metals. Natural orbitals from a simple correlation calculation like MP2 or a calculation with the MRCI module are usually a good choice and easily generated. For example: ```{literalinclude} ../../orca_working_input/C06S09_358.inp :language: orca ``` Now examine the occupation numbers of the natural orbitals (you will find that in the output of the MP2 part of the calculation): ``` Natural Orbital Occupation Numbers: N[ 0] = 2.00000000 N[ 1] = 2.00000000 N[ 2] = 1.98676733 N[ 3] = 1.97726840 N[ 4] = 1.97500109 N[ 5] = 1.96759239 N[ 6] = 1.96423113 N[ 7] = 1.93719340 N[ 8] = 0.05427454 N[ 9] = 0.02555886 N[ 10] = 0.02530580 N[ 11] = 0.01358500 N[ 12] = 0.01096092 N[ 13] = 0.01028129 N[ 14] = 0.00702048 N[ 15] = 0.00627820 ``` A rule of thumb is that orbitals with occupation numbers between 1.98 and 0.02 should be in the active space. Thus, in the present case we speculate that a 10 electrons in 8 orbitals active space would be appropriate for the CASSCF of the ground state. Let's try: ```{literalinclude} ../../orca_working_input/C06S09_359.inp :language: orca ``` If we run that calculation, it converges and produces the following: ``` MACRO-ITERATION 10: --- Inactive Energy E0 = -82.97337099 Eh E(CAS)= -113.889438276 Eh DE= -0.000000807 --- Energy gap subspaces: Ext-Act = -0.431 Act-Int = -0.240 N(occ)= 1.99763 1.99696 1.98360 1.97923 1.94253 0.05958 0.02153 0.01894 ||g|| = 0.000361782 Max(G)= 0.000189613 Rot=9,2 ---- THE CAS-SCF GRADIENT HAS CONVERGED ---- --- FINALIZING ORBITALS --- ---- DOING ONE FINAL ITERATION FOR PRINTING ---- --- Forming Natural Orbitals --- Canonicalize Internal Space --- Canonicalize External Space ``` From which we see that we had *two orbitals too many in the active space* with occupation numbers very close to two. The presence of barely correlated orbitals (occupation close to 0.0 or 2.0) can cause convergence problems. Their inclusion in the active space does not significantly change the energy and it might better to omit these orbitals from the start. In the present case, we re-run the CASSCF with 6 active electrons in six orbitals. The result is: ``` MACRO-ITERATION 2: --- Inactive Energy E0 = -101.16144179 Eh E(CAS)= -113.882700257 Eh DE= -0.012049926 --- Energy gap subspaces: Ext-Act = -0.411 Act-Int = -0.142 N(occ)= 1.98172 1.97921 1.94092 0.05983 0.02089 0.01743 ||g|| = 0.052811635 Max(G)= 0.025065586 Rot=19,7 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.160674186 Max(X)(5,4) = -0.128053569 --- SFit(Active Orbitals) MACRO-ITERATION 3: --- Inactive Energy E0 = -100.78371592 Eh E(CAS)= -113.885011169 Eh DE= -0.002310912 --- Energy gap subspaces: Ext-Act = -0.434 Act-Int = -0.199 N(occ)= 1.98150 1.97909 1.94143 0.05924 0.02108 0.01766 ||g|| = 0.017438409 Max(G)= 0.009231446 Rot=10,4 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.050337699 Max(X)(6,2) = -0.033671129 --- SFit(Active Orbitals) MACRO-ITERATION 4: --- Inactive Energy E0 = -100.72313195 Eh E(CAS)= -113.885258854 Eh DE= -0.000247685 --- Energy gap subspaces: Ext-Act = -0.438 Act-Int = -0.219 N(occ)= 1.98141 1.97918 1.94178 0.05886 0.02102 0.01776 ||g|| = 0.009726271 Max(G)= 0.004281706 Rot=9,2 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.031123960 Max(X)(22,9) = 0.015789781 --- SFit(Active Orbitals) MACRO-ITERATION 5: --- Inactive Energy E0 = -100.65264536 Eh E(CAS)= -113.885424851 Eh DE= -0.000165997 --- Energy gap subspaces: Ext-Act = -0.440 Act-Int = -0.238 N(occ)= 1.98140 1.97918 1.94202 0.05857 0.02105 0.01776 ||g|| = 0.006606671 Max(G)= 0.003548636 Rot=9,2 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.019988497 Max(X)(6,2) = -0.014410848 --- SFit(Active Orbitals) MACRO-ITERATION 6: --- Inactive Energy E0 = -100.56070274 Eh E(CAS)= -113.885549550 Eh DE= -0.000124699 --- Energy gap subspaces: Ext-Act = -0.440 Act-Int = -0.268 N(occ)= 1.98138 1.97925 1.94206 0.05849 0.02104 0.01778 ||g|| = 0.004483296 Max(G)= 0.002939015 Rot=9,2 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.011383690 Max(X)(5,4) = 0.005997355 --- SFit(Active Orbitals) MACRO-ITERATION 7: --- Inactive Energy E0 = -100.52522560 Eh E(CAS)= -113.885583124 Eh DE= -0.000033574 --- Energy gap subspaces: Ext-Act = -0.437 Act-Int = -0.283 N(occ)= 1.98132 1.97929 1.94192 0.05861 0.02103 0.01783 ||g|| = 0.002275031 Max(G)= -0.001215398 Rot=10,4 --- Orbital Update [SuperCI(PT)] --- Canonicalize Internal Space --- Canonicalize External Space --- SX_PT (Skipped TA=0 IT=0): ||X|| = 0.002106033 Max(X)(19,10) = -0.001056121 --- SFit(Active Orbitals) MACRO-ITERATION 8: --- Inactive Energy E0 = -100.52457962 Eh E(CAS)= -113.885584276 Eh DE= -0.000001152 --- Energy gap subspaces: Ext-Act = -0.438 Act-Int = -0.283 N(occ)= 1.98134 1.97931 1.94184 0.05868 0.02101 0.01781 ||g|| = 0.000752012 Max(G)= -0.000357510 Rot=13,4 ---- THE CAS-SCF GRADIENT HAS CONVERGED ---- --- FINALIZING ORBITALS --- ---- DOING ONE FINAL ITERATION FOR PRINTING ---- --- Forming Natural Orbitals --- Canonicalize Internal Space --- Canonicalize External Space MACRO-ITERATION 9: --- Inactive Energy E0 = -100.52457962 Eh --- All densities will be recomputed E(CAS)= -113.885584276 Eh DE= -0.000000000 --- Energy gap subspaces: Ext-Act = -0.858 Act-Int = -0.283 N(occ)= 1.98172 1.97932 1.94207 0.05845 0.02100 0.01743 ||g|| = 0.000752012 Max(G)= -0.000327367 Rot=12,4 -------------- CASSCF RESULTS -------------- Final CASSCF energy : -113.885584276 Eh -3098.9843 eV ``` The calculation converges very quickly and the occupation numbers show you that all of these orbitals are actually needed in the active space. The omission of the two orbitals from the active space came at an increase of the energy by $\sim$4 mEh which seems to be tolerable. Let's look what we have in the active space in figure {numref}`fig:CAS-H2CO`. (fig:CAS-H2CO)= :::{subfigure} ABC|DEF :subcaptions: below :class-grid: outline ![MO5](../../images/71.*) ![MO6](../../images/72.*) ![MO7](../../images/73.*) ![MO10](../../images/74.*) ![MO9](../../images/75.*) ![MO8](../../images/76.*) Orbitals of the active space for the CASSCF(6,6) calculation of H$_2$CO. ::: Thus, we can see that we got a fairly nice result: our calculation has correlated the in-plane oxygen lone pair, the C-O $\sigma$ and the C-O $\pi$ bond. For each strongly occupied bonding orbital, there is an accompanying weakly occupied antibonding orbital in the active space that is characterized by one more node. In particular, the correlating lone pair and the C-O $\sigma^{\ast }$ orbital would have been hard to find with any other procedure than the one chosen based on natural orbitals. We have now done it blindly and looked at the orbitals only after the CASSCF --- a better approach is normally to look at the starting orbitals *before* you enter a potentially expensive CASSCF calculation. If you have bonding/antibonding pairs in the active space plus perhaps the singly-occupied MOs of the system you probably have chosen a reasonable active space. We can play the game now somewhat more seriously and optimize the geometry of the molecule using a reasonable basis set: ```{literalinclude} ../../orca_working_input/C06S09_360.inp :language: orca ``` and get: ``` --------------------------------------------------------------------------- Redundant Internal Coordinates --- Optimized Parameters --- (Angstroem and degrees) Definition OldVal dE/dq Step FinalVal ---------------------------------------------------------------------------- 1. B(O 1,C 0) 1.2101 0.000259 -0.0002 1.2100 2. B(H 2,C 0) 1.0942 -0.000029 0.0001 1.0943 3. B(H 3,C 0) 1.0942 -0.000029 0.0001 1.0943 4. A(O 1,C 0,H 3) 122.07 0.000023 -0.00 122.07 5. A(H 2,C 0,H 3) 115.85 -0.000046 0.01 115.86 6. A(O 1,C 0,H 2) 122.07 0.000023 -0.00 122.07 7. I(O 1,H 3,H 2,C 0) -0.00 -0.000000 0.00 -0.00 ---------------------------------------------------------------------------- ``` Let us compare to MP2 geometries (this job was actually run first): ```{literalinclude} ../../orca_working_input/C06S10_361.inp :language: orca ``` ``` ------------------------------------------------------------------------- Redundant Internal Coordinates --- Optimized Parameters --- (Angstroem and degrees) Definition OldVal dE/dq Step FinalVal ---------------------------------------------------------------------------- 1. B(O 1,C 0) 1.2127 0.000374 -0.0002 1.2125 2. B(H 2,C 0) 1.0991 -0.000031 0.0001 1.0992 3. B(H 3,C 0) 1.0991 -0.000031 0.0001 1.0992 4. A(O 1,C 0,H 3) 121.77 0.000023 -0.00 121.77 5. A(H 2,C 0,H 3) 116.45 -0.000046 0.01 116.46 6. A(O 1,C 0,H 2) 121.77 0.000023 -0.00 121.77 7. I(O 1,H 3,H 2,C 0) -0.00 -0.000000 0.00 -0.00 ---------------------------------------------------------------------------- ``` The results are actually extremely similar (better than 1 pm agreement). Compare to RHF: ``` --------------------------------------------------------------------------- Redundant Internal Coordinates --- Optimized Parameters --- (Angstroem and degrees) Definition OldVal dE/dq Step FinalVal ---------------------------------------------------------------------------- 1. B(O 1,C 0) 1.1784 -0.000164 0.0001 1.1785 2. B(H 2,C 0) 1.0921 0.000010 -0.0000 1.0921 3. B(H 3,C 0) 1.0921 0.000010 -0.0000 1.0921 4. A(O 1,C 0,H 3) 121.93 -0.000003 -0.00 121.93 5. A(H 2,C 0,H 3) 116.13 0.000005 0.00 116.13 6. A(O 1,C 0,H 2) 121.93 -0.000003 -0.00 121.93 7. I(O 1,H 3,H 2,C 0) 0.00 0.000000 -0.00 -0.00 ---------------------------------------------------------------------------- ``` Thus, one can observe that the correlation brought in by CASSCF or MP2 has an important effect on the C$=$O distance ($\sim$4 pm), while the rest of the geometry is not much affected. **More on the technical use of the CASSCF program.** The most elementary input information which is always required for CASSCF calculations is the specification of the number of active electrons and orbitals. ```orca %casscf nel 4 # number of active space electrons norb 6 # number of active orbitals end ``` The CASSCF program in ORCA can average states of several multiplicities. The multiplicities are given as a list. For each multiplicity the number of roots should be specified: ```orca %casscf mult 1,3 # here: multiplicities singlet and triplet nroots 4,2 # four singlets, two triplets end ``` If the symmetry handling in ORCA is enabled (`! UseSym`) each multiplicity block must have an irreducible representation assigned. Numbers corresponding to the "irrep\" within a given symmetry are printed in the output of ORCA. ```orca %casscf mult 1,3 # here: multiplicities singlet and triplet irrep 0,1 # here: irrep for each mult. block (mandatory!) nroots 4,2 # four singlets, two triplets ``` Several roots and multiplicities usually imply a state average CASSCF (SA-CASSCF) calculation. Note that the program by default chooses equal weights for the multiplicity blocks. Roots within a given block have equal weight. Users can define a custom weighting scheme for the multiplicity blocks and roots: ```orca %casscf mult 1,3 # here: multiplicities singlet and triplet nroots 4,2 # four singlets, two triplets bweight 2,1 # singlets and triplets weighted 2:1 weights[0] = 0.5,0.2,0.2,0.2 # singlet weights weights[1] = 0.7,0.3 # triplet weights end ``` The program automatically normalizes these weights such that the sum over all weights is unity. If convergence on an excited state is desired then the `weights[0]` array may look like `0.0,0.0,1.0 `(this would optimize the orbitals for the third excited state. If several states cross during the orbital optimization this will ultimately cause convergence problems. We note passing that the converged orbitals of the state averaged procedure are a compromise for the set of states. ORCA by default only prints the SA-CASSCF gradient norm. State-specific gradients are summarized at the end of the calculation with the keyword `PrintGState`. ```orca %casscf ... printgstate true # optional printing of the state-specific orbital gradients end ``` **Orbital optimization methods.** In the following we discuss the available options for orbital optimization. A number of convergence problems can be resolved changing the guess orbitals. **The following keywords are optional and should only be used facing severe convergence difficulties.** Aside from the `SuperCI_PT` (default),{cite}`kollmarSXPT2019` several orbital optimization methods (list below) are implemented. ```orca # Keywords to be used as Orbstep/Switchstep SuperCI_PT # perturbative SuperCI (first order) SuperCI # SuperCI (first order) DIIS # DIIS (first order) KDIIS # KDIIS (first order) SOSCF # approx. Newton-Raphson (first order) NR # augmented Hessian Newton-Raphson # unfolded two-step procedure # - still not true second order ``` The different convergers have different strengths. First order method are cheap but typically require more iterations compared to second order methods. When the gradient is far off from convergence the program uses the converger defined as `orbstep` while close to convergence the `switchstep` is used. The actual criteria for switchstep are defined with the keywords `SwitchConv` and `SwitchIter`. ```orca %casscf OrbStep SuperCI # or any other from the list above SwitchStep DIIS # or any other from the list above SwitchConv 0.03 # gradient at which to switch SwitchIter 15 # iteration at which the switch takes place # irrespective of the gradient MaxIter 75 # Maximum number of macro-iterations end ``` Picking a convergence strategy, the program has to balance speed and robustness. The default strategy uses the `SuperCI_PT` as converger for `orbstep` and `switchstep`.{cite}`kollmarSXPT2019` This approach determines the elements $X_{pq}$ of the anti-Hermitean matrix used in the orbital update according to $$\mathbf{C}^{new}=\mathbf{C}^{old}e^{\mathbf{X} }$$ from first order perturbation theory using the Dyall-Hamiltonian {cite}`dyall1995journal` in zeroth order and a first-order perturbed wave function given as $\Psi^{(1) }=\sum_{pq}\Psi_p^q X_{qp}$ where the $\Psi_p^q$ represent singly excited functions obtained from the CASSCF wave function by excitation from orbital $\psi_p$ to orbital $\psi_q$. The `SuperCI_PT` is robust with respect to orbitals that are exactly doubly occupied or empty. Rotations with orbital close to this critical occupations can further be eliminated with the keyword `DThresh` (default=1e-6). However, the method is quiet aggressive in the orbital optimization. In some cases, such as basis set projection or `PATOM` guess (intrinsic basis set projection), the program might pick a step-size that is too big. Then restricting the step-size via the keyword `MaxRot` (default=0.2) might be useful. The keywords `DThresh` and `MaxRot` described below are specific to `SuperCI_PT`. For many users, `MaxRot` is less palpable than level shifting. Therefore, the present version allows level shifts as well. **In contrast to other convergers, level shifts are not needed and highly discouraged.** With the exception of `GradScaling` (*vide infra*), other damping techniques described further below do not apply to the `SuperCI_PT`. ```orca MaxRot 0.05 # cap stepsize for SuperCI_PT DThresh 1e-6 # thresh for critical occupation ``` In case of convergence problems with the default settings, it is recommended to try the combination of `orbstep SuperCI` and `switchstep DIIS`, which in conjuction with a large level shift (2 Eh), which may be immediately successful. The proposed scheme typically requires more iterations. Moreoever, in contrast to the SuperCI(PT), the `SuperCI`, `DIIS` and `KDIIS` should not be used when the active orbitals have an occupation of exactly 2.0 or 0.0! The `DIIS` may sometimes converge slowly or "trail" towards the end such that real convergence is never reached. The `KDIIS` {cite}`kollmar1996chem,kollmar1997int` --- based on perturbation theory --- is an approximation to the regular `DIIS` procedure avoiding redundant rotations. Both DIIS schemes avoid linear dependencies in the expansion space. ```orca MaxDIIS 15 # max. no of DIIS vectors to keep DIISThresh 1e-7 # overlap criteria for linear dependency ``` The combination of `SuperCI` and `DIIS` (switchstep) is particularly suited to protect the active space composition. Adjusting the level shift will do the job. **Here, level shift is the single most important lever to control convergence**. ```orca # default = dynamic level-shifting based on Ext-Act, Int-Act ShiftUp 2.0 # static up-shift the virtual orbitals ShiftDn 2.0 # static down-shift the internal orbitals MinShift 0.6 # minimum separation subspaces ``` Level-shift is particularly important if the active, inactive and virtual orbitals overlap in their orbital energies. The energy separation of the subspaces is printed in the output. Ideally, the entries `Ext-Act` and `Act-Int` should be positive and larger than 0.2 Eh. This will help the program to preserve your active space composition throughout the iterations. If no shift is specified in the input, ORCA will choose a level-shift to guarantee an energy separation between the subspaces (`MinShift`). ``` E(CAS)= -230.590325053 Eh DE= -0.000798832 --- Energy gap subspaces: Ext-Act = -0.244 Act-Int = -0.002 --- current l-shift: Up(Ext-Act) = 0.54 Dn(Act-Int) = 0.30 ``` In difficult cases the use of the Newton-Raphson method (`NR`) is recommended even if each individual iteration is considerably more expensive. It is strong towards the end but it would be a waste to start orbital optimization with the expensive `NR` method since its radius of quadratic convergence is quite small. The computationally cheaper alternative is the `SOSCF` procedure belonging to the family of quasi-Newton updates. Keep in mind that the Newton-Raphson is designed for optimization on a **convex surface** (Hessian is semidefinite). If the NR is switched on too early, there is a good chance that this condition is not fulfilled. In this case the program will complain about negative eigenvalues or diagonal elements of the Hessian as can be seen in the snippet below. The next optimization step is large and unpredictable. It is a wildcard that can get you closer to convergence or immediate divergence of the CASSCF procedure. ``` ||g|| = 0.771376945 Max(G)= 0.216712933 Rot=140,53 --- Orbital Update [ NR] Warning: NEGATIVE diagonal element D(81,53)= -4.733590 Warning: NEGATIVE diagonal element D(82,53)= -4.737955 ... ``` For larger system, the augmented Hessian equations are solved iteratively (NR iterations). The augmented Hessian is considered solved when the residual norm, $$, is small enough. Aside from the overall CASSCF convergence, negative eigenvalues affect these NR iterations. ``` --- Orbital Update [ NR] AugHess Tolerance (auto): Tol= 1.00e-07 AUGHESS-ITER 0: E= -0.174480747 = 0.558679452 AUGHESS-ITER 1: E= -0.308672359 = 0.468254671 AUGHESS-ITER 2: E= -0.434272813 = 0.286305469 AUGHESS-ITER 3: E= -0.439149451 = 0.286514628 AUGHESS-ITER 4: E= -0.605787445 = 0.191691955 AUGHESS-ITER 5: E= -0.607766529 = 0.310450670 AUGHESS-ITER 6: E= -0.611674930 = 0.141402593 AUGHESS-ITER 7: E= -0.623145299 = 0.394505306 AUGHESS-ITER 8: E= -0.658410333 = 0.166915094 AUGHESS-ITER 9: E= -0.790571374 = 4.722929453 AUGHESS-ITER 10: E= -0.790590554 = 4.716012014 AugHess: No convergence in the Davidson procedure ... ``` There are a number of refined NR settings that influence the convergence behavior on a non-convex energy surface. We mention the keywords for completeness and dis-encourage from changing the default settings. If overall convergence cannot be changed due to negative eigenvalues, it is recommended to delay the NR switchstep (`switchconv, switchiter`). This will require some trial and error, since the curvature of the surface is *a priori* not know. ```orca %casscf ... aughess Solver 0 # Davidson (default) 1 # Pople (pure NR steps) 2 # DIIS MaxIter 35 # max. no. of CI iters. MaxDim 35 # Davidson expansion space MaxDIIS 12 # max. number of DIIS vectors UseSubMatrixGuess true # diag a submatrix of the Hessian # as an initial guess NGuessMat 512 # size of initial guess matrix (part of # the Hessian exactly diagonalized) ExactDiagSwitch 512 # up to this dimension the Hessian # is exactly diagonalized (small problems) PrintLevel 1 # amount of output during AH iterations Tol 1e-6 # convergence tolerance Compress true # use compressed storage DiagShift 0.0 # shift of the diagonal elements of the # Hessian UseDiagPrec true # use the diagonal in updating SecShift 1e-4 # shift the higher roots in the Davidson # secular equations UpdateShift 0.5 # shift of the denominator in the # update of the AH coefficients end end ``` In general, convergence is strongly influenced by numerical noise, especially in the final iterations. One source of numerical noise is the incremental Fock build. Thus, it can help to enforce more frequent full Fock matrix formation. ```orca ResetFreq 1 # reset frequency for direct SCF ``` If the orbital change in the active space is small, the active Fock matrix in ORCA is approximated using the density matrix from the previous cycle saving a second Fock matrix build. However, this approximation might also be a source of numerical instability. The threshold "SwitchDens" can be set to zero to enable the exact build. The program default starts with a rather large value (1e-2) and will reduce this parameter automatically when necessary. ```orca switchdens 0.0001 # ~gtol * 0.1 ``` In all of the implemented orbital optimization schemes the step-size correlates with the gradient-norm. A constant damping factor can be set with the keyword `GradScaling`. Note, damping and level shifting techniques are not recommended for the default converger (`SuperCI_PT`). ```orca GradScaling 0.5 # constant damping in all steps ``` There are situations when the active space has been chosen carefully, but the initial gradient is still far off. To keep the "good\" active space, we can suppress all rotation but the inactive-external ones until the gradient-norm is small enough to continue safely. The threshold can be set with `FreezeIE` keyword. Once the components of the gradient in the inactive-external direction have a weight of less than `FreezeIE`, all constraints are lifted. ORCA by default freezes active rotations if the total gradient norm is larger than 1.0 and the active rotations have a weight of less than 5%. The feature can be turned off setting the threshold to zero. Similarly, rotations of the almost doubly occupied orbitals with the inactive orbitals can be damped using the threshold `FreezeActive`. Rotations of this type are damped as long as all their weight is smaller than `FreezeActive`. In contrast to the `ShiftDn`, it damps just the "troublesome\" parts of internal-active rotations. This option applies to all of the orbital optimization schemes but the `SuperCI_PT`. ```orca FreezeIE 0.4 # keep active space until int-ext rotation have # a contribution of less than 40% to the ||g|| FreezeActive 0.03 # keep almost doubly occupied orbitals as long as # their contribution is less than 3% to the ||g|| ``` If the calculation starts from a converged Hartree-Fock orbitals, the core orbitals should not change dramatically by the CASSCF optimization. Often trailing convergence is associated to rotations with low lying orbitals. Their contribution to the total energy is fairly small. With the keyword `FreezeGrad` these rotations can be omitted from the orbital optimization procedure. ```orca FreezeGrad 0.2 # omit hitting a gradient norm ||g|| <0.2 ``` The affected orbitals are printed at the startup of CASSCF. ``` FreezeGrad ... enabled if ||g|| is below 0.02 Note Convergence can be signaled if the reduced gradient reaches GTol Last frozen orbital ... 9 First deleted orbital ... 320 Once rotations with core and deleted orbitals are stabilized they will be damped. ``` By default rotations with frozencore (or deleted virtuals) are not omitted. If the option FreezeGrad is active, the ratio with respect to the total gradient is printed. ``` ||g|| = 0.001240414 Max(G)= -0.000431747 Rot=319,1 --- Option=FreezeGrad: ||g|| = 0.001081707 = 87.21% Omitting frozencore elements ``` **Using the RI Approximation.** Aside from the Fock matrices, integrals appearing in the orbital gradient and Hessian require substantial computation time. A good way to speed up the calculations at the expense of "only\" obtaining approximate results is to introduce the RI approximation. `TrafStep RI` approximated the aforementioned integrals. Here are sufficiently large auxiliary basis must be provided - ideally a /JK or /C. Further acceleration can be achieved approximating the Fock matrix construction with `!RIJCOSX` or `!RIJK` as described in section {ref}`sec:energygradients.casscf.approximation.typical`. More details can also be found in the CASSCF tutorial. Note that with ORCA 4.1, there are three destinct auxiliary basis slots, that need to be set if the auxiliary basis is defined via the `%basis` block. ```orca TrafoStep RI # RI used in transformation # Note: Needs an auxiliary basis for # AuxC slot. Exact # exact transformation (default) ``` **Monitoring the active space** (actconstraints)= During the iterations, the **active orbitals are chosen to maximize the overlap** with active orbitals from the previous iterations. Maximizing the overlap does not make any restrictions on the nature of the orbitals. Thus initially localized orbitals will stay localized and ordered, which is sometimes a desired feature e.g. in the density matrix renormalization group approach (DMRG). This feature is set with the keyword `ActConstraints` and is enabled by default (after the first 3 macroiterations). For some orbital optimization procedures, such as the SuperCI, natural orbitals are more advantageous. Therefore, the `ActConstraints` can be turned off in favor of natural orbital construction (see below). If the keyword is not set by the user, ORCA picks the best choice for the given orbital optimization step. ```orca ActConstraints 0 # no checks and no changes 1 # maximize overlap of active orbitals and check sanity. (default for DIIS) 2 # make natural orbitals in every iteration (default SuperCI) 3 # make canonical orbitals in every iteration 4 # localize orbitals ``` In addition to maximizing the overlap, **`"ActConstraints 1"` checks if the composition of the active space has changed** i.e. an orbital has been rotated out of the active space. In this case, ORCA aborts and stores the last valid set of orbitals. Below is an example error message. ``` --- Orbital Update [ DIIS] --- Failed to constrain active orbitals due to rotations: Rot( 37, 35) with OVL=0.960986 Rot( 38, 34) with OVL=0.842114 Rot( 43,104) with OVL=0.031938 ``` In the snippet above, the active space ranges from 37-43. The program reports that orbitals 37,38 and 43 have changed their character. The overlap of orbital 43 (active) with the previous set of active orbitals is just 3% and the program aborts. There are a number of reasons, why this happens in the calculation. If this error occurs constantly with the same orbitals, it is worthwhile to inspect the rotating partner orbitals (visualize). It might be sign that the active space is not balanced and should be extended. In many instances changing the level-shift or lowering switchconv is sufficient to protect the active space. In some cases, turning off the sanity check (`"ActConstraints 0"`) and re-rotating orbitals will bring CASSCF closer to convergence. Some problems can be avoided by a better design of the calculation. The CASSCF tutorial elaborates on the subject in more detail. There are situations such as parameter scans, where \"actconstraints\" is counter-productive and should be disabled. In other words, we want to allow changes in the active space composition. As an example, consider the rotation of ethylene around its double-bond represented by a CAS(2,2). Although the active space consists of the bonding and anti-bonding orbitals $\pi$-orbitals, their composition in terms of atomic orbitals changes from the eclipsed to the staggered conformation. Depending on the actual input settings (orbstep and number of scan points), this might trigger an abort. **Final orbitals options.** (actorbs)= Once the calculation has converged, ORCA will do a final macro-iteration, where the orbital are "finalized". For complete active spaces (CAS), these transformations do not alter the final energy and wavefunction. Note, that solutions from approximate CAS-CI solvers such as the ICE approach or the DMRG ansatz are affected by the final orbital transformation. The magnitude depends on the truncation level (e.g. `TGen`, `TVar` and `MaxM`) of the approximated wavefunction. The default final orbitals are canonical in the internal and external space with respect to state-averaged Fock operator. The active orbitals are chosen as natural orbitals. Other orbital choices are equally valid and can be selected for the individual subspaces. ```orca #internal space IntOrbs CanonOrbs # canonical LocOrbs # localized unchanged # no changes # partner orbitals for the active space based # on various concepts PMOS # based on orthogonalization tails. OSZ # based on oszillator orbital DOI # based on differential overlap #external space ExtOrbs CanonOrbs # canonical LocOrbs # localized unchanged # no changes # partner orbitals for the active space based # on various concepts PMOS # based on orthogonalization tails. OSZ # based on oszillator orbital DOI # based on differential overlap DoubleShell # based on the shell and angular momentum # of the highest active orbitals, the first few # virtual orbitals correspond to the doubled-shell. # All other virt. orbitals are canonicalized. # For 3d-metal complexes, these are the 4d orbitals! # For 4d-metal complexes, these are the 5d orbitals! # And so on... #active space ActOrbs NatOrbs # natural CanonOrbs # canonical LocOrbs # localized unchanged # no changes dOrbs # purify metal d-orbital and call the AILFT fOrbs # purify metal f-orbital and call the AILFT SDO # Single Determinant Orbitals: this is only possible if the # active space has a single hole or a single electron. # SDOs are then the unique choice of orbitals that simultaneously # turns each CASSCF root into a single determinant. ``` SDOs are specific for the active orbital space.{cite}`Lang2019` The set of options (`PMOS, OSZ, DOI, DoubelShell`) are specific for the inactive and external space. They aim to assist the extension of the current active space. All four options, re-design the first `NOrb` (number of active orbitals) next to the active space, while the remaining orbitals of the same subspace are canonical. The re-designed orbital are based on different concepts. - `PMOS` generates the bonding / anti-bonding partner orbitals for the chosen active space. It is based on the orthogonalization tail of the active orbitals. - `OSZ` generates a single orbital for each active orbital, that maximizes the dipole-dipole interaction. - `DOI` follows the same principle as OSZ, but the differential overlap is maximized instead. - `DoubleShell` is specific to the external space. The highest active MO or `DoubleShellMO` is analyzed. A set of orbitals with the same angular momentum but larger radial distribution is generated. Optionally, the four options above can be supplemented with a reference MO using the keyword `RefMO/DoubleShellMO`. The presence of `RefMO/DoubleShellMO` changes the default behavior. In case of `PMOS`, `OSZ` and `DOI`, all orbitals of the given subspace are chosen to maximize a single objective function with respect to the reference MO (must be active). This contrasts the default settings, where for each active MO an objective function is maximized and a single "best" orbital is picked. In other words, in the default setting, each active orbital has a corresponding "best" orbital in the selected subspace neighboring the active space. ```orca RefMO 17 # MO with number 17 (default =-1) DoubleShellMO 17 # MO with number 17 (default=-1) ``` The aforementioned options are aids and the resulting orbitals should be inspected prior extension of the active space. In particular the `PMOS` option is useful in the context of transition metal complexes to find suitable Ligand based orbitals. There are more options (`dorbs, forbs, DoubleShell`), that are specifically designed for coordination chemistry. A more detailed description is found in the CASSCF tutorial that supplements the manual. (PrintWFDets)= If the active space consists of a single set of metal d-orbitals, natural orbitals may be a mixture of the d-orbitals. The active orbitals are remixed to obtain \"pure\" d-orbitals (ligand field orbitals) if the `actorbs` is set to `dorbs`. The same holds for f-orbitals and the option `forbs`. Furthermore, the keyword `dorbs` automatically triggers the *ab initio* ligand field analysis (AILFT).{cite}`atanasov_modern_2011,lang_ailft_2020`The approach has been reported in a number of applications.{cite}`atanasov_detailed_2011,suturina_ailft_2015,aravena_periodic_2016,jung_ailft_2017,singh_covalency_2017,chakraborty_electronic_2017,chilkuri_revisiting_2017,chilkuri_aom_2020` Note that the actual representation depends on the chosen axis frame. It is thus recommended to align the molecule properly. For more details on the AILFT approach, we refer to the AILFT section ({ref}`sec:casscf.AbInitioLFT.detailed`), the original paper and the CASSCF tutorial, where examples are shown. For a few applications, a printing of the complete wavefunction is useful and can be requested. ```orca PrintWF 0 # (default) prints only the CFGs csf # Printing of the wavefunction in the basis of CSFs det # Printing of the wavefunction in the basis of Determinants ``` (cascistep)= The CI-step default setting is CSF based and is done in the present program by generating a partial "formula tape\" which is read in each CI iteration. The tape may become quite large beyond several hundred thousand CSFs which limits the applicability of the CASSCF module. The accelerated CI (ACCCI) has the same limitations, but uses a slightly different algorithm that handles multi-root calculations much more efficiently. For now, properties (spin-orbit coupling, g-Tensor\...) as well as NEVPT2 corrections are not available with ACCCI. Nevertheless, it is the recommended option to converge a CASSCF calculation with multiple roots. The resulting .gbw file may be used in a successive run to obtain properties or NEVPT2 corrections. Larger active spaces are tractable with the DMRG approach or the iterative configuration expansion (ICE) developed in our own group.{cite}`chilkuri_ice1,chilkuri_ice2` DMRG and ICE return approximate full CI results. The maximum size of the active space depends on the system and the required accuracy. Active spaces of 10--20 orbitals should be feasible with both approaches. The CASSCF tutorial covers examples with ACCCI and ICE as CI solvers. ```orca %casscf CIStep CSFCI # CSF based CI (default) ACCCI # CSF based CI solver with faster algorithm for multi-root calculations ICE # CSF based approximate CI -> ICE/CIPSI algorithm DMRGCI # density matrix renormalization group approach instead of the CI end ``` In the ICE approach, the computation of the coupling coefficients is time-consuming and by default repeated in every macro-iteration. To avoid the reconstruction, it is recommended to once generate a coupling coefficient library (`cclib`) and to use it for all of your ICE calculations. The details of the methodology and the `cclib` are described in the ICE section {ref}`sec:iceci.detailed`. Detailed settings for the conventional CI solvers (CSFCI, ACCCI, ICE) can be controlled in a sub-block. **Not all of the options and properties are available for CISteps apart from the default!** NEVPT2, transition densities and spin-dependent properties such as spin-orbit coupling are not yet available for `ACCCI` and `ICE`. ```orca %casscf ci MaxIter 64 # max. no. of CI iters. MaxDim 10 # Davidson expansion space = MaxDim * NRoots NGuessMat 512 # Initial guess matrix: 512x512 PrintLevel 3 # amount of output during CI iterations ETol 1e-10 # default 0.1*ETol in CASSCF RTol 1e-10 # default 0.1*ETol in CASSCF TGen 1e-4 # ICE generator thresh TVar 1e-11 # ICE selection thresh, default = TGen*1e-7 end ``` The `CI-step` DMRGCI interfaces to the **BLOCK** program developed in the group of G. K.-L. Chan {cite}`chan2002jcp,chan2004jcp,gosh2008jcp,sharma2012jcp`. A detailed description of the **BLOCK** program, its input parameters, general information and examples on the density matrix renormalization group (DMRG) approach, are available in the section {ref}`sec:dmrg.detailed` of the manual. The implementation of DMRG in **BLOCK** is fully spin-adapted. However, spin-densities and related properties are not available in the current version of the **BLOCK** code. To start a DMRG calculation add the keyword "`CIStep DMRGCI`" into a regular CASSCF input. ORCA will set default parameters and generate and input for the **BLOCK** program. In general, DMRG is not invariant to rotation in the active space. The program by **default** will run an automatic ordering procedure (`Fiedler`). More and refined options can be set in the `dmrg` sub-block of CASSCF --- see section {ref}`sec:dmrg.detailed` for a complete list of keywords. ```orca %casscf nel 8 norb 6 mult 3 CIStep DMRGCI # Detailed settings dmrg # more/refined options ... end end ``` It is highly recommended to start the calculation with split-localized orbitals. Any set of starting orbitals (gbw file) can be localized using the `orca_loc` program. Typing `orca_loc` in the shell will return a small help-file with details on how to setup an input for the localization. Examples for DMRG including the localization are in the corresponding section of the manual {ref}`sec:dmrg.detailed`. The utility program `orca_loc` is documented in section {ref}`sec:utilities.loc.detailed`. Split-localization refers to an independent localization of the internal and virtual part of the desired active orbitals. NOTE: - Let us stress again: it is strongly recommended to first LOOK at your orbitals and make sure that the ones that will enter the active space are really the ones that you want to be in the active space! Many problems can be solved by thinking about the desired physical contents of the reference space before starting a CASSCF. A poor choice of orbitals results in poor convergence or poor accuracy of the results! Choosing the active orbitals always requires chemical and physical insight into the molecules that you are studying! - Please try the program with default settings before playing with the more advanced options. If you encounter convergence problems, have a look into your output, read the warning and see how the gradient and energy evolves. Increasing `MaxIter` will not help in many cases. - Be careful with keywords such as `!tightscf`, `!verytightscf` and so on. These keywords set higher integral thresholds, which is a good idea, but also tighten the CASSCF convergence thresholds. If you do not need a tighter energy convergence, reset the criteria in the casscf block using `ETol`. For many applications an energy convergence beyond $10^{-7}$ is unnecessary. ## CASSCF Densities The one-particle electron and spin density can be stored on disk using the keyword `!KeepDens`. ORCA stores all densities in a container (**.densities** file on disk), which can be used in conjunction with `orca_plot` to plot the charge and spin densities. Please check Section {ref}`sec:utilities.plot.detailed` for more details on the procedure. The state-specific densities will have a name postfix that reflects the root, multiplicity and potentially irreducible representation of the state. Densities arising from a calculation with the spin-orbit coupling, will have an additional flag in the density container marking their origin (e.g. \"cas_qdsoc\" or \"nev_qdsoc\"). (sec:casscf.properties.detailed)= ## CASSCF Properties The CASSCF program is able to calculate UV transition, CD spectra, SOC, SSC, Zeeman splittings, EPR g-matrices and A-matrices (the latter implemented in the same way as in the DCD-CAS(2) method{cite}`Lang2019`), magnetization, magnetic susceptibility and MCD spectra. Note that the results for the Fermi contact contribution to A will not be reliable if the spin density is dominated by spin polarization, which is a dynamic correlation effect. The properties are exercised in more detail in the CASSCF tutorial. The techniques used to calculate SOC, and Zeeman splittings are identical to those implemented into the MRCI program. Input and keywords mimic the ones in the MRCI module described in section {ref}`sec:mrci.soc.detailed`. As an example, the input file to calculate g-values and HFC constants A of CO$^{+}$ is listed below: ```{literalinclude} ../../orca_working_input/C06S10_381.inp :language: orca ``` In addition to pseudo-spin 1/2 A-tensors for individual Kramers doublets, the CASSCF module also features the calculation of "intrinsic" HFC A-tensors for the whole lowest nonrelativistic spin multiplet in the effective Hamiltonian approach.{cite}`LangThesis` In contrast to the MRCI module, the CASSCF module also supports the calculation of susceptibility tensors at non-zero magnetic fields. The corresponding keywords are ```orca ... %casscf ... rel dosoc true dosusceptibility true susctensor_nfields 2 # number of user-defined magnetic fields susctensor_magfields[0] = 35000,0,0 # 1st user-defined magnetic field susctensor_magfields[1] = 70000,0,0 # 2nd user-defined magnetic field end end ``` This example input calculates the susceptibility tensor at the two (vector-valued!) magnetic fields (35000,0,0) and (70000,0,0) (in Gauss). Note that for practical reasons it is necessary to specify the number of user-defined magnetic fields using the keyword `susctensor_nfields`. Until ORCA 4.0 it was possible to access spin-spin couplings only via running CAS-CI type calculations in MRCI. Converged CASSCF orbitals can be read setting the following flags ```orca !MOREAD NOITER ALLOWRHF TZVPP TightSCF Bohrs %moinp "convergedCASSCF.gbw" %mrci ... TPre 0.0 citype mrci newblock 2 * excitations none refs CAS(9,8) end end soc DoSSC true # spin-spin coupling DoSOC true # spin-orbit coupling ... end end * xyz 1 2 C 0 0 0.0 O 0 0 2.3504 * ``` Starting with ORCA 4.1, spin-spin couplings are also directly accessible in the CASSCF module via the keyword `DoSSC true` in the `rel` subblock. Note that the calculation of SSC requires the definition of an auxiliary basis set (`AuxC` auxiliary basis set slot), since it is only implemented in conjunction with RI integrals. A common way to introduce dynamical correlation for the property computation, is to replace the energies entering the quasi-degenerate perturbation theory. If the NEVPT2 energy correction is computed in CASSCF, there will be additional printings where CASSCF energies are replaced by the more accurate NEVPT2 values. Alternatively, these diagonal energies can be taken from the input file similarly how it is described for the MRCI module. A more detailed documentation is presented in the MRCI property section. :::{note} - The program does NOT print the SOC matrix by default! To obtain SOCMEs at the CASSCF/NEVPT2/... levels, please set the `PrintLevel` in the `rel` block to at least 2. ::: (sec:casscf.AbInitioLFT.detailed)= ## 1- and 2-shell Abinitio Ligand Field Theory Starting from ORCA 5.0, ORCA features a 1- and 2-shell AILFT module. AILFT was originally developed for 1-shell d- and f- LFT problems {cite}`atanasov_modern_2011,lang_ailft_2020`. In ORCA 5.0 an extenion to 2-shell AILFT provides access to all common 1- and 2-shell AILFT problems namely: 1. Valence LFT problems, involving the d-, f-, sp-, ds- and df-shells 2. Core LFT problems involving the sd-, pd-, sf- and pf-shells become readily accesible. Requesting an CAS-AILFT calcultion withing the CASSCF module is provided in two ways: 1. Through the ActOrbs xOrbs keywords (e.g. xOrbs: dOrbs, fOrbs spOrbs, psOrbs, sdOrbs, dsOrbs, sfOrbs, fsOrbs, pdOrbs, dpOrbs, pfOrbs, fpOrbs, dfOrbs, fdOrbs) 2. Through the LFTCase keyword where particular LFT problems can be requested according to the above 1- and 2-shell combinations (e.g. LFTCase 3d, LFTCase 4f, LFTCase 1s3d, LFTCase 2p3d \...) Note: that the LFTCase keyword overwrites the ActOrbs keyword and as it will be discussed below provides a particular utility that simplifies the 2-shell AILFT input. A simple input for the Ni$^{2+}$ $d^8$ ion is provided below: ```orca !NEVPT2 def2-SVP def2-SVP/C %casscf nel 8 norb 5 ActOrbs dOrbs mult 3,1 nroots 10,15 rel dosoc true end end *xyz 2 3 Ni 0.0000000000 0.0000000000 0.0000000000 * ``` The programm after the CASSCF convergence will undergo few important steps and sanity checks which involve 1. an Orbital purification step 2. a Phase correction of the 1 and 2-electron integrals It is then important from the user's perspective to monitor that these steps have been succesfully performed. The relevant parts of the output are provided below: ``` ---- THE CAS-SCF GRADIENT HAS CONVERGED ---- --- FINALIZING ORBITALS --- ---- DOING ONE FINAL ITERATION FOR PRINTING ---- --- d-orbitals (depends on the molecular axis frame) L-Center: 0 Ni [0.000, 0.000, 0.000] --- The active space contains 5 d Orbitals : OK Setting 9 active MO to AO dz2 (11) Setting 10 active MO to AO dxz (12) Setting 11 active MO to AO dyz (13) Setting 12 active MO to AO dx2y2 (14) Setting 13 active MO to AO dxy (15) --- Canonicalize Internal Space --- Canonicalize External Space ... ============================= AB INITIO LIGAND FIELD THEORY d8 configuration 2 CI blocks MOs 9 to 13 ============================= Metal/Atom center is atom 0 orbital phases = 1.0 1.0 1.0 1.0 1.0 Metal/Atom d-orbital parts of active orbitals Shell 7 9 10 11 12 13 dz2 : 0.848522 -0.000000 0.000000 -0.000000 -0.000000 dxz : -0.000000 0.848522 0.000000 -0.000000 0.000000 dyz : 0.000000 0.000000 0.848522 -0.000000 -0.000000 dx2y2 : -0.000000 -0.000000 -0.000000 0.848522 0.000000 dxy : -0.000000 0.000000 -0.000000 0.000000 0.848522 Shell 8 9 10 11 12 13 dz2 : 0.300072 0.000000 -0.000000 0.000000 0.000000 dxz : 0.000000 0.300072 -0.000000 0.000000 -0.000000 dyz : -0.000000 -0.000000 0.300072 0.000000 0.000000 dx2y2 : 0.000000 0.000000 0.000000 0.300072 -0.000000 dxy : 0.000000 -0.000000 0.000000 -0.000000 0.300072 Adjusting phases of one-electron integrals ... done Adjusting phases of two-electron integrals ... done ``` In a subsequent step the program will 1. compute the AI Hamiltonian 2. construct the parameterized LFT Hamiltonian 3. and perform the fit The relevant output can be seen below: ``` Calculating ab initio Hamiltonian matrices ... ------------------------------------------------------------ NRoots (NEVPT2) for this block = 10 NEVPT2 correction for this block is calculated Full NEVPT2 Hamiltonian constructed Full NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ NRoots (NEVPT2) for this block = 15 NEVPT2 correction for this block is calculated Full NEVPT2 Hamiltonian constructed Full NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ Calculating fit matrices ... done Fitting ... done ``` In following the fitted 1-electron energies and SCP parameters also Racah parameters for 1-shells will be printed at the CASSCF and NEVPT2 levels of theory ``` ------------------------------ AILFT MATRIX ELEMENTS (CASSCF) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) : Orbital dz2 dxz dyz dx2-y2 dxy dz2 -8.111733 -0.000000 -0.000000 0.000000 -0.000000 dxz -0.000000 -8.111733 -0.000000 -0.000000 0.000000 dyz -0.000000 -0.000000 -8.111733 -0.000000 0.000000 dx2-y2 0.000000 -0.000000 -0.000000 -8.111733 -0.000000 dxy -0.000000 0.000000 0.000000 -0.000000 -8.111733 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F0dd(from 2el Ints) = 0.980960738 a.u. = 26.693 eV = 215296.0 cm**-1 (fixed) F2dd = 0.451725025 a.u. = 12.292 eV = 99142.2 cm**-1 F4dd = 0.280604669 a.u. = 7.636 eV = 61585.6 cm**-1 ------------------- Racah Parameters : ------------------- A(F0dd from 2el Ints) = 0.949782441 a.u. = 25.845 eV = 208453.2 cm**-1 B = 0.006037419 a.u. = 0.164 eV = 1325.1 cm**-1 C = 0.022270212 a.u. = 0.606 eV = 4887.7 cm**-1 C/B = 3.689 ------------------------------------------------------------------------------- ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) dz2 dxz dyz dx2-y2 dxy 1 0.000 0.0 -0.999978 -0.000164 -0.001934 0.005783 -0.002568 2 0.000 0.0 -0.005768 -0.000269 -0.000262 -0.999967 -0.005788 3 0.000 0.0 0.002600 0.000424 0.001046 0.005773 -0.999979 4 0.000 0.0 0.000241 -0.999246 -0.038831 0.000280 -0.000462 5 0.000 0.0 0.001930 0.038832 -0.999243 0.000246 -0.001022 Ligand field orbitals were stored in ni.3d.casscf.lft.gbw ... ------------------------------ AILFT MATRIX ELEMENTS (NEVPT2) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) : Orbital dz2 dxz dyz dx2-y2 dxy dz2 -8.118685 0.000000 0.000000 0.000005 -0.000000 dxz 0.000000 -8.118666 -0.000000 -0.000000 0.000000 dyz 0.000000 -0.000000 -8.118674 -0.000000 0.000000 dx2-y2 0.000005 -0.000000 -0.000000 -8.118676 0.000000 dxy -0.000000 0.000000 0.000000 0.000000 -8.118667 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F2dd = 0.415943380 a.u. = 11.318 eV = 91289.0 cm**-1 F4dd = 0.259145554 a.u. = 7.052 eV = 56875.9 cm**-1 ------------------- Racah Parameters : ------------------- B = 0.005550482 a.u. = 0.151 eV = 1218.2 cm**-1 C = 0.020567107 a.u. = 0.560 eV = 4514.0 cm**-1 C/B = 3.705 ------------------------------------------------------------------------------- ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) dz2 dxz dyz dx2-y2 dxy 1 0.000 0.0 -0.927589 0.002893 0.009447 0.373452 -0.003963 2 0.000 3.0 0.337599 0.005897 0.449370 0.827017 -0.010128 3 0.000 3.0 0.160011 -0.005672 -0.893243 0.420094 0.001383 4 0.001 4.4 0.000644 0.105816 -0.006612 -0.009602 -0.994317 5 0.001 4.6 0.001541 0.994348 -0.007084 -0.002573 0.105893 Ligand field orbitals were stored in ni.3d.nevpt2.lft.gbw ``` Note that: - At the CASSCF level F0 (and subsequently racah A) is computed from CASSCF 2-electron coulomb integrals - On the other hand at the NEVPT2 level F0 is not defined hence F0 and racah A are not printed. Below an alternative using the effective Slater exponnets will be provided. - The LFT orbitals are saved in \*.lft.gbw files which can be processed by the orca_plot to generate orbital visualization files. AILFT provides a Fit quality analysis (see the original paper {cite}`atanasov_modern_2011,lang_ailft_2020`) Note: That at the CASSCF level the AI matrix of free atoms and ions is exactly parameterized in the chosen LFT parameterization scheme. As a result the RMS AI-LFT fitting errors is expected to be practically zero. This is not the case when a correlation treatment is chosen like NEVPT2 and the errors are expected to be somewhat larger. The above is shown below: ``` Calculating statistical parameters ... done Reference energy AI-LFT = -38.134150221 au Reference energy AI = -38.134150221 au ------------------------------------------------ COMPARISON OF AB INITIO AND LIGAND FIELD RESULTS ------------------------------------------------ Block 1 --------- AI-Root 0: E(AI)= 0.000 eV -> LF-Root 0: 0.000 eV S= 0.998 Delta= -0.000 eV AI-Root 1: E(AI)= 0.000 eV -> LF-Root 1: 0.000 eV S= 0.981 Delta= -0.000 eV AI-Root 2: E(AI)= 0.000 eV -> LF-Root 2: 0.000 eV S= 0.980 Delta= -0.000 eV AI-Root 3: E(AI)= 0.000 eV -> LF-Root 3: 0.000 eV S= 0.773 Delta= 0.000 eV AI-Root 4: E(AI)= 0.000 eV -> LF-Root 4: 0.000 eV S= 0.774 Delta= -0.000 eV AI-Root 5: E(AI)= 0.000 eV -> LF-Root 5: 0.000 eV S= 0.985 Delta= -0.000 eV AI-Root 6: E(AI)= 0.000 eV -> LF-Root 6: 0.000 eV S= 0.986 Delta= -0.000 eV AI-Root 7: E(AI)= 2.464 eV -> LF-Root 7: 2.464 eV S= 0.998 Delta= -0.000 eV AI-Root 8: E(AI)= 2.464 eV -> LF-Root 8: 2.464 eV S= 0.998 Delta= -0.000 eV AI-Root 9: E(AI)= 2.464 eV -> LF-Root 9: 2.464 eV S= 1.000 Delta= -0.000 eV RMS error for this block = 0.000 eV = 0.0 cm**-1 Block 2 --------- AI-Root 0: E(AI)= 2.033 eV -> LF-Root 0: 2.033 eV S= 1.000 Delta= -0.000 eV AI-Root 1: E(AI)= 2.033 eV -> LF-Root 1: 2.033 eV S= 1.000 Delta= -0.000 eV AI-Root 2: E(AI)= 2.033 eV -> LF-Root 2: 2.033 eV S= 0.903 Delta= -0.000 eV AI-Root 3: E(AI)= 2.033 eV -> LF-Root 3: 2.033 eV S= 0.967 Delta= -0.000 eV AI-Root 4: E(AI)= 2.033 eV -> LF-Root 4: 2.033 eV S= 0.935 Delta= -0.000 eV AI-Root 5: E(AI)= 3.183 eV -> LF-Root 5: 3.183 eV S= 0.996 Delta= -0.000 eV AI-Root 6: E(AI)= 3.183 eV -> LF-Root 6: 3.183 eV S= 0.999 Delta= -0.000 eV AI-Root 7: E(AI)= 3.183 eV -> LF-Root 7: 3.183 eV S= 0.996 Delta= -0.000 eV AI-Root 8: E(AI)= 3.183 eV -> LF-Root 8: 3.183 eV S= 0.999 Delta= -0.000 eV AI-Root 9: E(AI)= 3.183 eV -> LF-Root 9: 3.183 eV S= 0.999 Delta= -0.000 eV AI-Root 10: E(AI)= 3.183 eV -> LF-Root 10: 3.183 eV S= 0.995 Delta= -0.000 eV AI-Root 11: E(AI)= 3.183 eV -> LF-Root 11: 3.183 eV S= 0.992 Delta= -0.000 eV AI-Root 12: E(AI)= 3.183 eV -> LF-Root 12: 3.183 eV S= 0.999 Delta= -0.000 eV AI-Root 13: E(AI)= 3.183 eV -> LF-Root 13: 3.183 eV S= 0.996 Delta= -0.000 eV AI-Root 14: E(AI)= 7.856 eV -> LF-Root 14: 7.856 eV S= 1.000 Delta= -0.000 eV RMS error for this block = 0.000 eV = 0.0 cm**-1 Total RMS error g= 0.000 eV = 0.0 cm**-1 Note: Dropped RMS error for the reference energy. Confidence intervals (95) computed from the square root of the diagonal elements of the covariance matrix: H(dz2 ,dz2 )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxz ,dz2 )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxz ,dxz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dyz ,dz2 )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dyz ,dxz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dyz ,dyz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dx2-y2,dz2 )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dx2-y2,dxz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dx2-y2,dyz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dx2-y2,dx2-y2)= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxy ,dz2 )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxy ,dxz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxy ,dyz )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxy ,dx2-y2)= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 H(dxy ,dxy )= 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 F2dd = 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 F4dd = 0.000000000 a.u. = 0.000 eV = 0.0 cm**-1 Pearson's correlation coefficient = 1.000 (should be close to 1) ``` ``` Calculating statistical parameters ... done Reference energy AI-LFT = -38.134345028 au Reference energy AI = -38.134150221 au ------------------------------------------------ COMPARISON OF AB INITIO AND LIGAND FIELD RESULTS ------------------------------------------------ Block 1 --------- AI-Root 0: E(AI)= -0.000 eV -> LF-Root 0: 0.000 eV S= 0.933 Delta= -0.000 eV AI-Root 1: E(AI)= 0.000 eV -> LF-Root 1: 0.000 eV S= 0.769 Delta= -0.000 eV AI-Root 2: E(AI)= 0.001 eV -> LF-Root 2: 0.000 eV S= 0.773 Delta= 0.000 eV AI-Root 3: E(AI)= 0.001 eV -> LF-Root 3: 0.000 eV S= 0.742 Delta= 0.001 eV AI-Root 4: E(AI)= 0.002 eV -> LF-Root 4: 0.000 eV S= 0.750 Delta= 0.001 eV AI-Root 5: E(AI)= 0.003 eV -> LF-Root 5: 0.001 eV S= 0.931 Delta= 0.002 eV AI-Root 6: E(AI)= 0.003 eV -> LF-Root 6: 0.001 eV S= 0.998 Delta= 0.003 eV AI-Root 7: E(AI)= 2.195 eV -> LF-Root 7: 2.266 eV S= 1.000 Delta= -0.070 eV AI-Root 8: E(AI)= 2.195 eV -> LF-Root 8: 2.266 eV S= 0.998 Delta= -0.070 eV AI-Root 9: E(AI)= 2.195 eV -> LF-Root 9: 2.266 eV S= 0.998 Delta= -0.070 eV RMS error for this block = 0.039 eV = 311.1 cm**-1 Block 2 --------- AI-Root 0: E(AI)= 1.812 eV -> LF-Root 0: 1.875 eV S= 0.825 Delta= -0.063 eV AI-Root 1: E(AI)= 1.812 eV -> LF-Root 1: 1.875 eV S= 0.938 Delta= -0.063 eV AI-Root 2: E(AI)= 1.812 eV -> LF-Root 2: 1.875 eV S= 1.000 Delta= -0.063 eV AI-Root 3: E(AI)= 1.812 eV -> LF-Root 3: 1.875 eV S= 1.000 Delta= -0.063 eV AI-Root 4: E(AI)= 1.812 eV -> LF-Root 4: 1.875 eV S= 0.773 Delta= -0.063 eV AI-Root 5: E(AI)= 2.987 eV -> LF-Root 5: 2.932 eV S= 0.955 Delta= 0.056 eV AI-Root 6: E(AI)= 2.987 eV -> LF-Root 6: 2.932 eV S= 0.910 Delta= 0.055 eV AI-Root 7: E(AI)= 2.988 eV -> LF-Root 7: 2.932 eV S= 0.874 Delta= 0.056 eV AI-Root 8: E(AI)= 2.988 eV -> LF-Root 8: 2.932 eV S= 0.792 Delta= 0.056 eV AI-Root 9: E(AI)= 2.988 eV -> LF-Root 9: 2.932 eV S= 0.796 Delta= 0.056 eV AI-Root 10: E(AI)= 2.988 eV -> LF-Root 10: 2.932 eV S= 0.808 Delta= 0.056 eV AI-Root 11: E(AI)= 2.988 eV -> LF-Root 11: 2.932 eV S= 0.971 Delta= 0.056 eV AI-Root 12: E(AI)= 2.988 eV -> LF-Root 12: 2.932 eV S= 0.994 Delta= 0.056 eV AI-Root 13: E(AI)= 2.989 eV -> LF-Root 13: 2.932 eV S= 0.994 Delta= 0.057 eV AI-Root 14: E(AI)= 7.122 eV -> LF-Root 14: 7.241 eV S= 1.000 Delta= -0.119 eV RMS error for this block = 0.064 eV = 519.6 cm**-1 Total RMS error g= 0.057 eV = 457.2 cm**-1 Note: Dropped RMS error for the reference energy. Confidence intervals (95) computed from the square root of the diagonal elements of the covariance matrix: H(dz2 ,dz2 )= 0.000523387 a.u. = 0.014 eV = 114.9 cm**-1 H(dxz ,dz2 )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dxz ,dxz )= 0.000523387 a.u. = 0.014 eV = 114.9 cm**-1 H(dyz ,dz2 )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dyz ,dxz )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dyz ,dyz )= 0.000523387 a.u. = 0.014 eV = 114.9 cm**-1 H(dx2-y2,dz2 )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dx2-y2,dxz )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dx2-y2,dyz )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dx2-y2,dx2-y2)= 0.000523387 a.u. = 0.014 eV = 114.9 cm**-1 H(dxy ,dz2 )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dxy ,dxz )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dxy ,dyz )= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dxy ,dx2-y2)= 0.000393965 a.u. = 0.011 eV = 86.5 cm**-1 H(dxy ,dxy )= 0.000523387 a.u. = 0.014 eV = 114.9 cm**-1 F2dd = 0.002351095 a.u. = 0.064 eV = 516.0 cm**-1 F4dd = 0.003038264 a.u. = 0.083 eV = 666.8 cm**-1 Pearson's correlation coefficient = 1.000 (should be close to 1) ``` Several utilities are offered for more specialized tasks that provide better control of the AILFT inputs and outputs: - Skipping orbital otimization or reading in previously computed orbitals can be requested in two ways: 1. By the !NoIter keyword in the command line 2. By the AILFT_SkipOrbOpt in the ailft block (see example below) - Estimating F0 SCPs or Racah A from single zeta Slater Exponents can be requested from the AILFT_EffectveSlaterExponents true keyword in the ailft block - For the above task the knowledge of the principle quantum numbers is required. This can be specidied in two ways: 1. By the AILFT_PQN x keyword in the ailft block (x=3 for 3d) 2. By the LFTCase x keyword (LFTCase 3d, ommiting in this way the ActOrbs dOrbs keyword) Let us see how all the above translates in the above example of the Ni$^{2+}$ $d^8$ ion ```orca !NoIter NEVPT2 def2-SVP def2-SVP/C %casscf nel 8 norb 5 LFTCase 3d mult 3,1 nroots 10,15 ailft #AILFT_SkipOrbOpt true (An alternative to NoIter) #AILFT_PQNs 3 (Works together with ActOrbs dOrbs as an alternative to LFTCase 3d) AILFT_SlaterExponents true end rel dosoc true end end *xyz 2 3 Ni 0.0000000000 0.0000000000 0.0000000000 * ``` By running the above input the fitted 1-electron energies and SCP parameters also Racah parameters for 1-shells will be printed at the CASSCF and NEVPT2 levels of theory, including F0s and Racah A as estimated from single zeta effective Slater exponents from the fitted F2dd SCPs. ``` ------------------------------ AILFT MATRIX ELEMENTS (CASSCF) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) : Orbital dz2 dxz dyz dx2-y2 dxy dz2 -7.974440 0.000000 -0.000000 0.000000 -0.000000 dxz 0.000000 -7.974440 -0.000000 -0.000000 0.000000 dyz -0.000000 -0.000000 -7.974440 -0.000000 0.000000 dx2-y2 0.000000 -0.000000 -0.000000 -7.974440 0.000000 dxy -0.000000 0.000000 0.000000 0.000000 -7.974440 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- ------------------------------------------------------------------------------- Computed Single Zeta Slater Effective Exponents ... ------------------------------------------------------------------------------- kd(SZ)(from F2dd) = 3.134434893 a.u. ------------------------------------------------------------------------------- Computed F0s from Single Zeta Slater Effective Exponents ... ------------------------------------------------------------------------------- F0dd(from F2dd kd(SZ)) = 0.809116820 a.u. = 22.017 eV = 177580.6 cm**-1 F2dd = 0.427107567 a.u. = 11.622 eV = 93739.3 cm**-1 F4dd = 0.264174069 a.u. = 7.189 eV = 57979.5 cm**-1 ------------------------------------------------------------------------------- ------------------- Racah Parameters : ------------------- A(from F2dd kd(SZ)) = 0.779764145 a.u. = 21.218 eV = 171138.4 cm**-1 B = 0.005721310 a.u. = 0.156 eV = 1255.7 cm**-1 C = 0.020966196 a.u. = 0.571 eV = 4601.5 cm**-1 C/B = 3.665 ------------------------------------------------------------------------------- ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) dz2 dxz dyz dx2-y2 dxy 1 0.000 0.0 -0.999981 0.000787 -0.001735 0.004390 -0.003810 2 0.000 0.0 0.003811 0.000493 0.000955 0.000497 -0.999992 3 0.000 0.0 0.004388 0.000169 0.000080 0.999990 0.000514 4 0.000 0.0 -0.000750 -0.999810 -0.019460 0.000174 -0.000514 5 0.000 0.0 -0.001754 -0.019459 0.999809 -0.000069 0.000939 Ligand field orbitals were stored in ni.3d.casscf.lft.gbw ... ------------------------------ AILFT MATRIX ELEMENTS (NEVPT2) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) : Orbital dz2 dxz dyz dx2-y2 dxy dz2 -7.974812 -0.000000 -0.000000 0.000002 -0.000000 dxz -0.000000 -7.974860 0.000000 0.000000 0.000000 dyz -0.000000 0.000000 -7.974856 -0.000000 0.000000 dx2-y2 0.000002 0.000000 -0.000000 -7.974837 0.000000 dxy -0.000000 0.000000 0.000000 0.000000 -7.974837 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- ------------------------------------------------------------------------------- Computed Single Zeta Slater Effective Exponents ... ------------------------------------------------------------------------------- kd(SZ)(from F2dd) = 3.126330433 a.u. ------------------------------------------------------------------------------- Computed F0s from Single Zeta Slater Effective Exponents ... ------------------------------------------------------------------------------- F0dd(from F2dd kd(SZ)) = 0.807024751 a.u. = 21.960 eV = 177121.5 cm**-1 F2dd = 0.426003229 a.u. = 11.592 eV = 93496.9 cm**-1 F4dd = 0.262001371 a.u. = 7.129 eV = 57502.7 cm**-1 ------------------------------------------------------------------------------- ------------------- Racah Parameters : ------------------- A(from F2dd kd(SZ)) = 0.777913487 a.u. = 21.168 eV = 170732.3 cm**-1 B = 0.005723406 a.u. = 0.156 eV = 1256.1 cm**-1 C = 0.020793760 a.u. = 0.566 eV = 4563.7 cm**-1 C/B = 3.633 ------------------------------------------------------------------------------- ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) dz2 dxz dyz dx2-y2 dxy 1 0.000 0.0 -0.000608 -0.999795 0.017518 0.010019 0.001244 2 0.000 0.9 0.000003 -0.017532 -0.999696 -0.003661 0.016920 3 0.001 4.9 0.065092 -0.009854 0.004898 -0.995782 0.063719 4 0.001 5.0 -0.004126 0.002173 0.016617 0.063640 0.997824 5 0.001 10.5 0.997871 0.000042 -0.000237 0.065225 -0.000030 Ligand field orbitals were stored in ni.3d.nevpt2.lft.gbw ``` It is also possible to treat only the High Spin states in the d- and f- 1-shell AILFT. Note that not all the cases can be treated as this renters the different SCP parameters undefined. In the beggining, AILFT will check whether such case is detected and will drop a warning message For example in the case of the Fe$^{2+}$ $d^6$ ion with an input like the following: ```orca !NoIter def2-SVP def2-SVP/C %casscf nel 6 norb 5 mult 5 nroots 5 LFTCase 3d end *xyz 2 3 Fe 0.0000000000 0.0000000000 0.0000000000 * ``` the following Warning will be printed in the beggining of the calculation ``` WARNING: In AILFT F2dd remains undefined when considering only the HS Multiplicity block Be Careful with the results! TIP: If possible use in addition a LS Multiplicity Block ``` Spin orbit coupling effects (SOC) can be introduced by parametrizing the effective SOC constant $\zeta$. As long as SOC is requested in the rel CASSCF block the respective requested shell effective SOC constant $\zeta$ will be computed at the end of every AILFT calculation Hence in the above examples one gets: ``` ---------------------------------------------- SPIN ORBIT COUPLING (based on CASSCF orbitals) ---------------------------------------------- AI-SOC-X integrals (cm-1) 0 1 2 3 4 0 0.000000 0.000000 1078.568273 -0.000000 -0.000000 1 -0.000000 0.000000 -0.000000 -0.000000 622.711683 2 -1078.568273 0.000000 0.000000 -622.711683 0.000000 3 0.000000 0.000000 622.711683 0.000000 -0.000000 4 0.000000 -622.711683 -0.000000 0.000000 0.000000 AI-SOC-Y integrals (cm-1) 0 1 2 3 4 0 -0.000000 -1078.568273 -0.000000 0.000000 0.000000 1 1078.568273 0.000000 -0.000000 -622.711683 -0.000000 2 0.000000 0.000000 0.000000 -0.000000 -622.711683 3 -0.000000 622.711683 0.000000 0.000000 -0.000000 4 -0.000000 0.000000 622.711683 0.000000 0.000000 AI-SOC-Z integrals (cm-1) 0 1 2 3 4 0 0.000000 -0.000000 0.000000 -0.000000 0.000000 1 0.000000 -0.000000 -622.711683 -0.000000 0.000000 2 -0.000000 622.711683 0.000000 -0.000000 -0.000000 3 0.000000 0.000000 0.000000 -0.000000 -1245.423365 4 -0.000000 -0.000000 0.000000 1245.423365 0.000000 Fit to the SOC matrix elements a = 15.000000 b = 1.158 eV = 9340.7 cm**-1 SOC constant zeta = 0.077 eV = 622.7 cm**-1 LF-SOC-X integrals (cm-1) 0 1 2 3 4 0 0.000000 -0.000000 1078.568273 -0.000000 -0.000000 1 0.000000 0.000000 -0.000000 -0.000000 622.711683 2 -1078.568273 0.000000 0.000000 -622.711683 -0.000000 3 0.000000 0.000000 622.711683 0.000000 -0.000000 4 0.000000 -622.711683 0.000000 0.000000 0.000000 LF-SOC-Y integrals (cm-1) 0 1 2 3 4 0 0.000000 -1078.568273 -0.000000 -0.000000 -0.000000 1 1078.568273 0.000000 -0.000000 -622.711683 -0.000000 2 0.000000 0.000000 0.000000 -0.000000 -622.711683 3 0.000000 622.711683 0.000000 0.000000 -0.000000 4 0.000000 0.000000 622.711683 0.000000 0.000000 LF-SOC-Z integrals (cm-1) 0 1 2 3 4 0 0.000000 -0.000000 -0.000000 -0.000000 -0.000000 1 0.000000 0.000000 -622.711683 -0.000000 -0.000000 2 0.000000 622.711683 0.000000 -0.000000 -0.000000 3 0.000000 0.000000 0.000000 0.000000 -1245.423365 4 0.000000 0.000000 0.000000 1245.423365 0.000000 RMS error of nonzero matrix elements = 0.0 cm**-1 -----SOC-CONSTANTS----- ---All Values in cm-1--- ZETA_D = 622.71 ------------------------ ``` Starting fron ORCA 5.0 it is also possible in addition to CASSCF and NEVPT2 to employ DCDCAS(2) and Hermitian QD-NEVPT2 Abinitio Hamiltonians in AILFT Example inputs are provided below for DCDCAS(2): ```orca !def2-SVP def2-SVP/C %casscf nel 8 norb 5 actorbs dorbs mult 3,1 nroots 10,15 dcdcas true corrorder 2 rel dosoc true end end *xyz 2 3 Ni 0.0000000000 0.0000000000 0.0000000000 * ``` and Hermitian QD-NEVPT2: ```orca !def2-SVP def2-SVP/C %casscf nel 8 norb 5 actorbs dorbs mult 3,1 nroots 10,15 PTMethod sc_nevpt2 PTSettings QDType QD_VanVleck end rel dosoc true end end *xyz 2 3 Ni 0.0000000000 0.0000000000 0.0000000000 * ``` Running the above inputs the respective DCDCAS(2) and Hermitian QD-NEVPT2 Hamiltonians will be processed: ``` Calculating ab initio Hamiltonian matrices ... ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 0 and order = 0 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 0 and order = 1 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 0 and order = 2 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 1 and order = 0 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 1 and order = 1 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ DCDCAS correction for this block/order is calculated DCDCAS Hamiltonian of block = 1 and order = 2 is passed DCDCAS Hamiltonian diagonalized ------------------------------------------------------------ ... Calculating ab initio Hamiltonian matrices ... ------------------------------------------------------------ Hermitian QD-NEVPT2 correction for this block is calculated Hermitian QD-NEVPT2 Hamiltonian of block = 0 is passed Hermitian QD-NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ Hermitian QD-NEVPT2 correction for this block is calculated Hermitian QD-NEVPT2 Hamiltonian of block = 1 is passed Hermitian QD-NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ ``` It should be noted that NEVPT2 and Hermitian QD-NEVPT2 AILFT require a complete saturation of the excitation space. This implies that if less roots than the required are requested the AILFT analyis will be skipped in these cases. This is on the contrary not the case in CASSCF or DCDCAS(2) in which AILFT can operate under incomplete saturation of the excitation space. ``` Calculating ab initio Hamiltonian matrices ... WARNING: Number of NEVPT2 roots for block 0 (5) is not equal to the number of CASCI CSFs (10)! Skipping AILFT analysis with NEVPT2 energies! WARNING: Number of NEVPT2 roots for block 1 (2) is not equal to the number of CASCI CSFs (15)! Skipping AILFT analysis with NEVPT2 energies! ``` In a similar fashion one can request a 2-shell AILFT calulation. For this purpose the recomended steps are the following: - In a first step the valence active space orbitals are optimized in the framework of SA-CASSCF calculation.e.g. the 3d MOs in a core 1s3d or 2p3d AILFT calculation, or the f MOs in an 4f5d AILFT calculation) - In a second step the relevant core or virtual orbitals are rotated into the active space and the chosen CASCI/AILFT problem is solved by saturating the excitation space with all the involved excitations/multiplicity. - In the most of the cases the excitation space of two multiplicities the High-Spin one and the subsequent Low-Spin one are enouph for a succesfull fitting of the parameters It should be noted that 2-shell AILFT ivolves a 2-step fitting process following a bottom up shell angular momentum approach : 1. At first when possible an intra-shell fitting is performed 2. In following the respective effective Slater exponents are derived 3. In a last step an inter-shell fitting is performed and all the computed/fitted parameters are printed This implies that: - the flag of computing effective Slater exponents is always on by default in 2-shell AILFT - the desired LFT problem is best requested by the LFTCase keywords (e.g. LFTCase 1s3d) Let look at the case of 1s3d LFT problem of the Ni$^{2+}$ $d^8$ ion. A relevant input is provided below: ```orca !NoIter NEVPT2 def2-SVP def2-SVP/C %method frozencore fc_none end %scf rotate {0,8,90} end end %casscf nel 10 norb 6 mult 3,1 nroots 100,100 LFTCase 1s3d rel dosoc true end end *xyz 2 3 Ni 0.0000000000 0.0000000000 0.0000000000 * ``` Like in 1-shell AILFT, 2-shells AILFT starts with a sanity check ``` ---- THE CAS-SCF GRADIENT HAS CONVERGED ---- --- FINALIZING ORBITALS --- ---- DOING ONE FINAL ITERATION FOR PRINTING ---- --- sd-orbitals (depends on the molecular axis frame) L-Center: 0 Ni [0.000, 0.000, 0.000] L-Center: 0 Ni Active Orbital 0 is the s orbital ; l = 0 ; active shell = 0 L-Center: 0 Ni Active Orbital 1 is one d orbital ; l = 2 ; active shell = 7 L-Center: 0 Ni Active Orbital 2 is one d orbital ; l = 2 ; active shell = 7 L-Center: 0 Ni Active Orbital 3 is one d orbital ; l = 2 ; active shell = 7 L-Center: 0 Ni Active Orbital 4 is one d orbital ; l = 2 ; active shell = 7 L-Center: 0 Ni Active Orbital 5 is one d orbital ; l = 2 ; active shell = 7 --- The active space contains 1 s orbitals and 5 d orbitals : OK Setting 8 active MO to AO s (0) Setting 9 active MO to AO dz2 (11) Setting 10 active MO to AO dxz (12) Setting 11 active MO to AO dyz (13) Setting 12 active MO to AO dx2y2 (14) Setting 13 active MO to AO dxy (15) --- Canonicalize Internal Space --- Canonicalize External Space ``` In following the AI-LFT Hamiltonians are constructed and the LFT parameters are fitted at the CASSCF and at the NEVPT2 levels of theory ``` ------------------------------ AILFT MATRIX ELEMENTS (CASSCF) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) with V(0,0) fixed : 0 1 2 3 4 5 0 -334.652557 -0.000000 0.000000 -0.000000 -0.000000 -0.000000 1 -0.000000 -10.085777 0.000000 0.000000 -0.000000 -0.000000 2 0.000000 0.000000 -10.085777 -0.000000 -0.000000 -0.000000 3 -0.000000 0.000000 -0.000000 -10.085777 0.000000 0.000000 4 -0.000000 -0.000000 -0.000000 0.000000 -10.085777 -0.000000 5 -0.000000 -0.000000 -0.000000 0.000000 -0.000000 -10.085777 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F0ss = 17.137989284 a.u. = 466.348 eV = 3761353.9 cm**-1 F0dd = 0.809116820 a.u. = 22.017 eV = 177580.6 cm**-1 F2dd = 0.427107567 a.u. = 11.622 eV = 93739.3 cm**-1 F4dd = 0.278548413 a.u. = 7.580 eV = 61134.3 cm**-1 F0sd = 1.285327742 a.u. = 34.976 eV = 282096.8 cm**-1 G2sd = 0.001928559 a.u. = 0.052 eV = 423.3 cm**-1 R2sddd = 0.003748289 a.u. = 0.102 eV = 822.7 cm**-1 ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) s dz2 dxz dyz dx2-y2 dxz 1 0.000 0.0 -1.000000 -0.000000 0.000000 -0.000000 -0.000000 0.000000 2 8831.911 71234174.2 0.000000 -0.006282 0.006583 -0.271975 0.962261 0.000000 3 8831.911 71234174.2 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 4 8831.911 71234174.2 0.000000 -0.075508 -0.016185 -0.959327 -0.271528 0.000000 5 8831.911 71234174.2 -0.000000 0.071391 -0.997365 0.008467 0.009682 0.000000 6 8831.911 71234174.2 -0.000000 0.994566 0.070404 -0.075159 -0.015232 -0.000000 Ligand field orbitals were stored ni.1s3d.casscf.lft.gbw ... ------------------------------ AILFT MATRIX ELEMENTS (NEVPT2) -------------------------------- Ligand field one-electron matrix VLFT (a.u.) with V(0,0) fixed : 0 1 2 3 4 5 0 -334.652557 -0.000000 0.000000 0.000000 0.000000 0.000000 1 -0.000000 -10.298799 0.000008 0.000006 -0.000002 -0.000002 2 0.000000 0.000008 -10.293700 -0.000042 -0.000008 -0.000006 3 0.000000 0.000006 -0.000042 -10.293839 -0.000001 0.000008 4 0.000000 -0.000002 -0.000008 -0.000001 -10.293927 -0.000033 5 0.000000 -0.000002 -0.000006 0.000008 -0.000033 -10.293511 ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F0ss = 16.677683435 a.u. = 453.823 eV = 3660328.4 cm**-1 F0dd = 0.806859685 a.u. = 21.956 eV = 177085.2 cm**-1 F2dd = 0.425916096 a.u. = 11.590 eV = 93477.8 cm**-1 F4dd = 0.277771367 a.u. = 7.559 eV = 60963.8 cm**-1 F0sd = 1.424742585 a.u. = 38.769 eV = 312694.9 cm**-1 G2sd = 0.025339297 a.u. = 0.690 eV = 5561.3 cm**-1 R2sddd = 0.003739506 a.u. = 0.102 eV = 820.7 cm**-1 ----------------------------------------- The ligand field one electron eigenfunctions: ----------------------------------------- Orbital Energy (eV) Energy(cm-1) s dz2 dxz dyz dx2-y2 dxz 1 0.000 0.0 1.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000 2 8826.114 71187421.2 -0.000000 0.999998 -0.001607 -0.001229 0.000405 0.000329 3 8826.247 71188489.9 -0.000000 0.000330 -0.040203 -0.027545 -0.995773 -0.077852 4 8826.249 71188507.4 -0.000000 0.001631 0.264542 0.963492 -0.035726 -0.020544 5 8826.254 71188542.9 0.000000 -0.001223 -0.962932 0.264866 0.034492 -0.037635 6 8826.258 71188582.5 0.000000 -0.000317 -0.034070 0.027728 -0.077265 0.996042 Ligand field orbitals were stored in ni.1s3d.nevpt2.lft.gbw ``` As discussed above saturation of the excitation space is a requirement also in the case of 2-shell AILFT. It is usually enough to specify a large number of roots for two multiplicites (e.g. 100 singlets and triplets in the above example) The exact number of roots will be automatically detected. ``` Multiplicity ... 3 #(Configurations) ... 15 #(CSFs) ... 15 #(Roots) ... 100 WARNING (ORCA_CASSCF): NRoots > NCSFs. Adjusting to maximum number of roots. Please check the output carefully! #(Roots) ... 15 ROOT=0 WEIGHT= 0.033333 ROOT=1 WEIGHT= 0.033333 ROOT=2 WEIGHT= 0.033333 ROOT=3 WEIGHT= 0.033333 ROOT=4 WEIGHT= 0.033333 ROOT=5 WEIGHT= 0.033333 ROOT=6 WEIGHT= 0.033333 ROOT=7 WEIGHT= 0.033333 ROOT=8 WEIGHT= 0.033333 ROOT=9 WEIGHT= 0.033333 ROOT=10 WEIGHT= 0.033333 ROOT=11 WEIGHT= 0.033333 ROOT=12 WEIGHT= 0.033333 ROOT=13 WEIGHT= 0.033333 ROOT=14 WEIGHT= 0.033333 BLOCK 2 WEIGHT= 0.5000 Multiplicity ... 1 #(Configurations) ... 21 #(CSFs) ... 21 #(Roots) ... 100 WARNING (ORCA_CASSCF): NRoots > NCSFs. Adjusting to maximum number of roots. Please check the output carefully! #(Roots) ... 21 ROOT=0 WEIGHT= 0.023810 ROOT=1 WEIGHT= 0.023810 ROOT=2 WEIGHT= 0.023810 ROOT=3 WEIGHT= 0.023810 ROOT=4 WEIGHT= 0.023810 ROOT=5 WEIGHT= 0.023810 ROOT=6 WEIGHT= 0.023810 ROOT=7 WEIGHT= 0.023810 ROOT=8 WEIGHT= 0.023810 ROOT=9 WEIGHT= 0.023810 ROOT=10 WEIGHT= 0.023810 ROOT=11 WEIGHT= 0.023810 ROOT=12 WEIGHT= 0.023810 ROOT=13 WEIGHT= 0.023810 ROOT=14 WEIGHT= 0.023810 ROOT=15 WEIGHT= 0.023810 ROOT=16 WEIGHT= 0.023810 ROOT=17 WEIGHT= 0.023810 ROOT=18 WEIGHT= 0.023810 ROOT=19 WEIGHT= 0.023810 ROOT=20 WEIGHT= 0.023810 ``` However very often the required number of states to be computed in the framework of NEVPT2 type of calculations are quite large. In these cases a Hamiltonian reduction process on the basis of the Restrictive Active Space (RAS) is required. In fact all LFT parameters can be determined by considering up to double excitations from the donor-shell. Let us consider the 2p3d case of the Fe$^{2+}$ $d^6$ ion. Saturation of the active space requires to consider 70 triplet and 378 singlet states. Restriction of the active space to only up to double excitations from the 2p-shell results in 65 quintet and 330 triplet states. The Hamiltonian reduction can be requested in the ailft block: ```orca ailft AILFT_Dimension 2 # Up to double excitations from the donor shell Other options: 1 Up to single excitations from the donor shell 0 No excitations from the donor shell end ``` Hence the relevant input can be now formulated as: ```orca !NoIter MOREAD DKH2 DKH-def2-TZVP def2-TZVP/C NEVPT2 %moinp "fe.3d.gbw" %pal nprocs 16 end %method frozencore fc_none end %scf rotate {2,6,90} {3,7,90} {4,8,90} end end %casscf nel 12 norb 8 mult 5,3 nroots 65,330 LFTCase 2p3d ailft AILFT_Dimension 2 end rel dosoc true GTensor false DoDTensor false end end *xyz 2 5 Fe 0.0000000000 0.0000000000 0.0000000000 * ``` Note that before running the above calculation: - An initial SA-CASSCF calculation has been performed on the valence states of Fe in the 3d active space. These orbitals (fe.3d.gbw) are read in - The computation of the g- and D- tensors is switched off. This is recomended if the magnetism analysis is not required - As core spectroscopy is targeted the frozen core is switched off At the NEVPT2 part the reduced AI and LFT Hamiltonians will be constructed ``` Calculating ab initio Hamiltonian matrices ... ------------------------------------------------------------ NRoots (NEVPT2) for this block = 65 NEVPT2 correction for this block is calculated Full NEVPT2 Hamiltonian constructed Full NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ ------------------------------------------------------------ NRoots (NEVPT2) for this block = 330 NEVPT2 correction for this block is calculated Full NEVPT2 Hamiltonian constructed Full NEVPT2 Hamiltonian diagonalized ------------------------------------------------------------ ``` As a result the CASSCF and NEVPT2 LFT parameters will be determined in the requested reduced basis ``` ------------------------------ AILFT MATRIX ELEMENTS (CASSCF) -------------------------------- ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F0pp = 3.889665545 a.u. = 105.843 eV = 853682.9 cm**-1 F2pp = 1.832901835 a.u. = 49.876 eV = 402275.5 cm**-1 F0dd = 0.767706319 a.u. = 20.890 eV = 168492.1 cm**-1 F2dd = 0.405248254 a.u. = 11.027 eV = 88941.7 cm**-1 F4dd = 0.264292339 a.u. = 7.192 eV = 58005.5 cm**-1 F0pd = 1.174187892 a.u. = 31.951 eV = 257704.5 cm**-1 F2pd = 0.220959621 a.u. = 6.013 eV = 48495.0 cm**-1 G1pd = 0.157243007 a.u. = 4.279 eV = 34510.9 cm**-1 G3pd = 0.089473762 a.u. = 2.435 eV = 19637.2 cm**-1 ... ------------------------------ AILFT MATRIX ELEMENTS (NEVPT2) -------------------------------- ------------------------------------------------- Slater-Condon Parameters (electronic repulsion) : ------------------------------------------------- F0pp = 3.527444471 a.u. = 95.987 eV = 774184.6 cm**-1 F2pp = 1.906123325 a.u. = 51.868 eV = 418345.7 cm**-1 F0dd = 0.808248242 a.u. = 21.994 eV = 177390.0 cm**-1 F2dd = 0.426649072 a.u. = 11.610 eV = 93638.6 cm**-1 F4dd = 0.278249395 a.u. = 7.572 eV = 61068.7 cm**-1 F0pd = 1.657417509 a.u. = 45.101 eV = 363761.1 cm**-1 F2pd = 0.198646995 a.u. = 5.405 eV = 43598.0 cm**-1 G1pd = 0.206111579 a.u. = 5.609 eV = 45236.3 cm**-1 G3pd = 0.128174221 a.u. = 3.488 eV = 28131.0 cm**-1 ... ``` In the above example inclusion of SOC will result in the computation of the effective SOC $\zeta$ constants of both p and d shells: ``` -----SOC-CONSTANTS----- ---All Values in cm-1--- ZETA_P = 65018.19 ZETA_D = 453.53 ------------------------ ``` One important feature of 1- and in particular of the 2-shell AILFT is that it is connected to the standalone orca_lft multiplet program. Hence every succesfull AILFT calculation will automatically construct relevant inputs for the orca_lft. For example in the avove 2p3d case of the Fe$^{2+}$ $d^6$ ion the following inputs will be constructed ``` fe.2p3d.casscf.lft.inp fe.2p3d.nevpt2.lft.inp ``` with the NEVPT2 one looking like this: ``` %lft #-----Parameters------ NEl= 12 Shell_PQN= 0, 2, 3, 0 Mult= 5, 3 NRoots= 65, 330 #-------------------- #---Slater-Condon Parameters--- #---All Values in eV--- PARAMETERS F0pp = 95.9866 F2pp = 51.8683 F0dd = 21.9936 F2dd = 11.6097 F4dd = 7.5716 F0pd = 45.1006 F2pd = 5.4055 G1pd = 5.6086 G3pd = 3.4878 end #-------------------- #---LFT-Matrix Elelemnts--- #---All Values in eV--- FUNCTIONS 0 0 " 0.0000" 1 0 " -0.0146" 1 1 " 0.0947" 2 0 " 0.0210" 2 1 " 0.0061" 2 2 " 0.1155" 3 0 " 0.0000" 3 1 " -0.0000" 3 2 " -0.0000" 3 3 "1086.2398" 4 0 " 0.0000" 4 1 " -0.0000" 4 2 " -0.0000" 4 3 " 0.0323" 4 4 "1086.2181" 5 0 " -0.0000" 5 1 " 0.0000" 5 2 " 0.0000" 5 3 " -0.0356" 5 4 " -0.0154" 5 5 "1086.1183" 6 0 " 0.0000" 6 1 " -0.0000" 6 2 " -0.0000" 6 3 " -0.0767" 6 4 " -0.0080" 6 5 " 0.0341" 6 6 "1086.1219" 7 0 " -0.0000" 7 1 " 0.0000" 7 2 " 0.0000" 7 3 " 0.0050" 7 4 " -0.0023" 7 5 " 0.0026" 7 6 " -0.0038" 7 7 "1086.0688" end #-------------------- #---SOC-CONSTANTS--- #---All Values in eV--- PARAMETERS ZETA_P = 8.06 ZETA_D = 0.06 end #-------------------- #---SPECTRA/PROPERTIES--- DoABS true #------------------------ end *xyz Charge Multiplicity Atom 0.00 0.00 0.00 * ``` Further details regardin orca_lft can be found in the orca_lft section ({ref}`sec:utilities.lft.detailed`) and the orca_lft tutorial. (sec:casscf.CoreExitedStates.detailed)= ## Core excited states with (C/R)ASCI/NEVPT2 Starting from ORCA 4.1, a CASCI/NEVPT2 protocol can be used to compute core excited spectra, namely X-ray absorption (XAS) and resonant inelastic scattering (RIXS) spectra. RASCI calculations can also be easily specified. The XAS/RIXS spectra calculations requires two steps: - In a first step one needs to optimize the valence active space orbitals in the framework of SA-CASSCF calculations, e.g. including valence excited states in the range between 6 to 15 eV. - In a second step the relevant core orbitals are rotated into the active space and the (C/R)ASCI/NEVPT2 problem is solved by saturating the excitation space with singly core-excited electronic configurations using the previously optimized sets of orbitals Further information can be found in reference{cite}`2018xascascinevpt2` A relevant input for Fe L-edge XAS calculation of a Fe(III) complex like $\textrm{Fe(acac)}_3$ is given below for CASCI/NEVPT2: ```orca %moinp "Fe_acac3_casscf.gbw" %scf rotate {4,89,90,0,0} {3,88,90,0,0} {2,87,90,0,0} end end %rel picturechange true FiniteNuc true end %method FrozenCore FC_NONE end # CASSCF/NEVPT2 on the valence and L-edge excited states %casscf nel 11 norb 8 mult 6,4 nroots 16,173 maxiter 1 # account for spin-orbit coupling rel DoSOC true end # adding the dynamical correlation with NEVPT2 PTMethod SC\_NEVPT2 end * xyz 0 6 ... * ``` For RASCI/NEVPT2 calculations the valence d AS is set to RAS2. The RAS3 space is usually set empty. The RAS1 space contains the previously rotated core orbitals. To generate a single core hole, the number of maximum holes in the RAS1 space must be set to 1. Accordingly, the maximum number of particles in the RAS3 space must be 0. The RASCI input should thus look the following ```orca %casscf nel 11 norb 8 ... refs ras(11:3 1/5/0 0) # (Nel: NRAS1 MaxHoles / NRAS2 / NRAS3 MaxParticles) end ... end ``` As it is explicitly described in the respective ROCIS section RIXS spectra can be requested by the following keywords: ```orca RIXS true # Request RIXS calculation (NoSOC) RIXSSOC true # Request RIXS calculation (with SOC) Elastic true # Request RIXS calculation (Elastic) ``` Please consult section {ref}`sec:rocis.rixs.detailed` for processing and analyzing the generated spectra. (sec:casscf.XES.detailed)= ## CASCI-XES Starting from ORCA 5.0 likewise to RASCI-XES (see section {ref}`sec:mrci.soc.x-ray.spectroscopy.detailed`) orca features a CASCI-XES protocol. Likewise to the RASCI-XES the CASCI-XES calculations requires two steps: - In a first step one needs to optimize the valence active space orbitals in the framework of SA-CASSCF calculations, e.g. including valence excited states in the range between 6 to 15 eV. - In a second step the relevant core orbitals e.g metal 1s and 3p are rotated into the active space and the CASCI problem is solved for the ionized system by saturating the excitation space with singly core-excited electronic configurations using the previously optimized sets of orbitals. In CASCI-XES this can be acheived by defining reference configurations. or via the RAS functionality. The XESSOC calculation is called by speciyfing the following keywords in the rel block: ```orca rel XESSOC true XASMOs Number of the rotated 1s MO end ``` Following a SA-CASSCF calculation: ```orca ! def2-SVP def2-SVP/C ZORA CPCM PAL8 ! NormalPrint ! NoLoewdin NoMulliken %casscf nel 5 norb 5 nroots 1, 24 mult 6, 4 #rel #dosoc true #end end * xyz -3 6 Fe 0 0 0 Cl 2.40 0 0 Cl -2.40 0 0 Cl 0 2.40 0 Cl 0 -2.40 0 Cl 0 0 2.40 Cl 0 0 -2.40 * ``` A relevant input for Fe XES calculation of a Fe(III) complex like FeCl$_6$ is given below: ```orca ! def2-SVP def2-SVP/C ZORA CPCM PAL8 ! NormalPrint ! MOREAD ! NoLoewdin NoMulliken %moinp "FeCl6_casscf.gbw" %scf #Rotate the 1s and 3p orbitals below the SOMOs by using the rotate option rotate {0,59,90} {36,60,90} {37,61,90} {38,62,90} end end %casscf nel 12 norb 9 nroots 1000, 1000 mult 7, 5 maxiter 1 refs 1 2 2 2 0 0 1 2 2 1 2 2 2 0 0 2 1 2 1 2 2 2 0 0 2 2 1 1 2 2 2 0 1 0 2 2 1 2 2 2 0 1 1 1 2 1 2 2 2 0 1 1 2 1 1 2 2 2 0 1 2 0 2 ... 2 2 2 2 2 0 0 2 0 2 2 2 2 2 0 1 0 1 2 2 2 2 2 0 1 1 0 2 2 2 2 2 0 2 0 0 2 2 2 2 2 1 0 0 1 2 2 2 2 2 1 0 1 0 2 2 2 2 2 1 1 0 0 2 2 2 2 2 2 0 0 0 end DoDipoleVelocity true DoHigherMoments true rel dosoc true XESSOC true XASMOs 59 DoDTensor false DoVelocity true end end %method SpecialGridAtoms 26 #Increase the radial integration accuracy on Fe SpecialGridIntAcc 7 #Requested radial integration accuracy values end * xyz -2 5 Fe 0 0 0 Cl 2.40 0 0 Cl -2.40 0 0 Cl 0 2.40 0 Cl 0 -2.40 0 Cl 0 0 2.40 Cl 0 0 -2.40 * ``` In ORCA 6 all these inputs between like MRCI, CASSCF and LFT have been unified so like in {ref}`sec:mrci.soc.x-ray.spectroscopy.detailed` one can also specify the respective RAS excitation space as following: ```orca %casscf nel 12 norb 9 nroots 1000, 1000 mult 7, 5 maxiter 1 refs ras(12:4 1/5/ 0 0) end DoDipoleVelocity true DoHigherMoments true rel dosoc true XESSOC true XASMOs 59 DoDTensor false DoVelocity true end end ``` In the above inputs one notes that the exact knowledge of the states to saturate the excitations space is not required. One only need to specify a large number (e.g. 1000) and the program will automatically detect the required 19 septet and 270 quintet states: ``` BLOCK 1 WEIGHT= 0.5000 Multiplicity ... 7 #(Configurations) ... 4 #(CSFs) ... 4 #(Roots) ... 1000 WARNING (ORCA_CASSCF): NRoots > NCSFs. Adjusting to maximum number of roots. Please check the output carefully! #(Roots) ... 4 ... BLOCK 2 WEIGHT= 0.5000 Multiplicity ... 5 #(Configurations) ... 89 #(CSFs) ... 105 #(Roots) ... 1000 WARNING (ORCA_CASSCF): NRoots > NCSFs. Adjusting to maximum number of roots. Please check the output carefully! #(Roots) ... 105 ``` By now running the above input for the 4 septet and the 81 quintet states the following output is generated ``` ------------------------------ Printing the XES spectrum ... ------------------------------ ------------------------------------------------------------------------------------- SPIN-ORBIT X-RAY EMISSION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS ------------------------------------------------------------------------------------- Transition Energy INT TX TY TZ 1 421 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 2 422 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 3 423 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 4 424 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 5 425 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 6 426 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 7 427 -> 5 7231.106 0.000000000000 0.00000 0.00000 0.00000 ... 2641 421 -> 25 7179.464 0.003027813485 0.00003 0.00001 0.00415 2642 422 -> 25 7179.464 0.003027769731 0.00114 0.00399 0.00001 2643 423 -> 25 7179.464 0.003027770766 0.00399 0.00114 0.00003 2644 424 -> 25 7179.464 0.000000000035 0.00000 0.00000 0.00000 2645 425 -> 25 7179.464 0.000000000000 0.00000 0.00000 0.00000 2646 426 -> 25 7179.464 0.000000000063 0.00000 0.00000 0.00000 2647 427 -> 25 7179.464 0.000000000001 0.00000 0.00000 0.00000 2648 428 -> 25 7179.520 0.000000000000 0.00000 0.00000 0.00000 2649 429 -> 25 7179.520 0.000000000000 0.00000 0.00000 0.00000 2650 430 -> 25 7179.520 0.000000000000 0.00000 0.00000 0.00000 2651 431 -> 25 7179.520 0.000000000000 0.00000 0.00000 0.00000 2652 432 -> 25 7179.520 0.000000000000 0.00000 0.00000 0.00000 2653 433 -> 25 7181.156 0.000000887886 0.00000 0.00000 0.00007 2654 434 -> 25 7181.156 0.000000887873 0.00001 0.00007 0.00000 2655 435 -> 25 7181.156 0.000000887873 0.00007 0.00001 0.00000 ... 54898 538 -> 420 7164.167 0.000077338518 0.00000 0.00000 0.00066 54899 539 -> 420 7164.167 0.000077338538 0.00047 0.00047 0.00000 54900 540 -> 420 7164.167 0.000077338525 0.00047 0.00047 0.00000 54901 541 -> 420 7164.186 0.000000000000 0.00000 0.00000 0.00000 54902 542 -> 420 7164.186 0.000000000000 0.00000 0.00000 0.00000 54903 543 -> 420 7164.186 0.000000000000 0.00000 0.00000 0.00000 54904 544 -> 420 7164.186 0.000000000000 0.00000 0.00000 0.00000 54905 545 -> 420 7164.186 0.000000000000 0.00000 0.00000 0.00000 54906 546 -> 420 7164.214 0.000000000000 0.00000 0.00000 0.00000 54907 547 -> 420 7164.214 0.000000000000 0.00000 0.00000 0.00000 54908 548 -> 420 7164.214 0.000000000000 0.00000 0.00000 0.00000 54909 549 -> 420 7164.214 0.000000000000 0.00000 0.00000 0.00000 54910 550 -> 420 7164.214 0.000050793492 0.00000 0.00000 0.00054 54911 551 -> 420 7164.214 0.000050793481 0.00020 0.00050 0.00000 54912 552 -> 420 7164.214 0.000050793480 0.00050 0.00020 0.00000 All Done ------------------------------------------------------------------------------------- ``` Finally by processing the .out file with orca_mapspc: ```orca orca_mapspc fecl6_xes.out XESSOC -x07140 -x17190 -w4.0 -eV -n10000 ``` and by plotting the resulted XES spectrum one will get the respective RASCI-XES spectrum presented in {numref}`fig:RASCI-XESSOC` Since Orca 6.0 the computed transition moments in the presence of SOC can be taken beyond the dipole approximation by using the OPS tool, check section ({ref}`sec:ops.detailed`) for details. All the possible approximations can be requested with the following commands: ```orca #Non-Relativistic/Relativistic treatment %casscf DoDipoleLength true DoDipoleVelocity true DoHigherMoments true DecomposeFoscLength true DecomposeFoscVelocity true DoFullSemiclassical true end ``` This wil generate a list of tables which for the case of XESSOC will look like the following: ``` ---------------------------------------------------------------------------------------------- SOC COMBINED ELECTRIC DIPOLE + MAGNETIC DIPOLE + ELECTRIC QUADRUPOLE X-RAY EMISSION SPECTRUM ---------------------------------------------------------------------------------------------- INT (fosc) --------------------------------------------------------- Transition Energy D2 m2 Q2 D2+m2+Q2 D2/TOT m2/TOT Q2/TOT (eV) (*1e6) (*1e6) ---------------------------------------------------------------------------------------------- ... ---------------------------------------------------------------------------------------------- SOC COMBINED ELECTRIC DIPOLE + MAGNETIC DIPOLE + ELECTRIC QUADRUPOLE X-RAY EMISSION SPECTRUM (Origin Adjusted) ---------------------------------------------------------------------------------------------- INT (fosc) --------------------------------------------------------- Transition Energy D2 m2 Q2 D2+m2+Q2 D2/TOT m2/TOT Q2/TOT (eV) (*1e6) (*1e6) ----------------------------------------------------------------------------------------------- ... ----------------------------------------------------------------------------------------------- SOC COMBINED ELECTRIC DIPOLE + MAGNETIC DIPOLE + ELECTRIC QUADRUPOLE X-RAY EMISSION SPECTRUM (Origin Independent, Velocity) ----------------------------------------------------------------------------------------------- INT (fosc) -------------------------------------------------------- Transition Energy D2 m2 Q2 D2+m2+Q2+DM+DO D2/TOT m2/TOT Q2/TOT (eV) (*1e6) (*1e6) ----------------------------------------------------------------------------------------------- ... ----------------------------------------------------------------------------------------------- SOC COMBINED ELECTRIC DIPOLE + MAGNETIC DIPOLE + ELECTRIC QUADRUPOLE X-RAY EMISSION SPECTRUM (Exact Formulation, Velocity) ------------------------------------------------------------------------------------------------ INT (fosc) --------------------------------------------------------------------------- Transition Energy D2 m2 Q2 Exact Osc. Str. D2/TOT m2/TOT Q2/TOT (eV) (*1e6) (*1e6) ------------------------------------------------------------------------------------------------ ``` Orca_mapspc can process all these files. The list of the relevant keywords is: ``` ABSQ #This will process non relativistic Origin Adjusted Absortion like spectra ABSOI #This will process the non relativistic Exact, Velocity Absortion like spectra SOCABSQ #This will process the SOC corrected Origin Adjusted Absortion like spectra SOCABSOI #This will process the SOC corrected Exact Velocity Absortion like spectra XESSOCQ #This will process the SOC corrected Origin Adjusted X-ray Emission spectra XESSOCOI #This will process the SOC corrected Exact Velocity X-ray Emission spectra ``` A more complete list can be found in ({ref}`sec:utilities.mapspc.detailed`) [^1]: For approximate full CI approaches, CASSCF is neither invariant to active-active rotations nor exactly size-consistent.