6.7 Coupled-Cluster Excited-State and Open-Shell Methods

6.7.1 Excited States via EOM-EE-CCSD and EOM-EE-OD

One can describe electronically excited states at a level of theory similar to that associated with coupled-cluster theory for the ground state by applying either linear response theory [351] or equation-of-motion methods [352]. A number of groups have demonstrated that excitation energies based on a coupled-cluster singles and doubles ground state are generally very accurate for states that are primarily single electron promotions. The error observed in calculated excitation energies to such states is typically 0.1–0.2 eV, with 0.3 eV as a conservative estimate, including both valence and Rydberg excited states. This, of course, assumes that a basis set large and flexible enough to describe the valence and Rydberg states is employed. The accuracy of excited state coupled-cluster methods is much lower for excited states that involve a substantial double excitation character, where errors may be 1 eV or even more. Such errors arise because the description of electron correlation of an excited state with substantial double excitation character requires higher truncation of the excitation operator. The description of these states can be improved by including triple excitations, as in the EOM(2,3) or EOM-CCSD(dT)/(fT) methods.

Q-Chem includes coupled-cluster methods for excited states based on the optimized orbital coupled-cluster doubles (OD), and the coupled cluster singles and doubles (CCSD) methods, described earlier. OD excitation energies have been shown to be essentially identical in numerical performance to CCSD excited states [353].

These methods, while far more computationally expensive than TDDFT, are nevertheless useful as proven high accuracy methods for the study of excited states of small molecules. Moreover, they are capable of describing both valence and Rydberg excited states, as well as states of a charge-transfer character. Also, when studying a series of related molecules it can be very useful to compare the performance of TDDFT and coupled-cluster theory for at least a small example to understand its performance. Along similar lines, the CIS(D) method described earlier as an economical correlation energy correction to CIS excitation energies is in fact an approximation to EOM-CCSD. It is useful to assess the performance of CIS(D) for a class of problems by benchmarking against the full coupled-cluster treatment. Finally, Q-Chem includes extensions of EOM methods to treat ionized or electron attachment systems, as well as di- and tri-radicals.

EOM-EE $\Psi (M_ s=0)=R(M_ s=0) \Psi _0(M_ s=0)$

$\begin{picture}(17,2)\thicklines \put( 0,1){\line(1,0){1}} \put( 0,2){\line(1,0){1}} \put(0.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put( 3,1){\line(1,0){1}} \put( 3,2){\line(1,0){1}} \put( 5,1){\line(1,0){1}} \put( 5,2){\line(1,0){1}} \put(3.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(3.48,1.8){\Large {\ensuremath{\downharpoonright }}} \put(5.48,0.8){\Large {\ensuremath{\downharpoonright }}} \put(5.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put( 8,1){\line(1,0){1}} \put( 8,2){\line(1,0){1}} \put(8.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(3.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(4.5,-0.7){$\Phi _ i^{a}$} \put(8.1,-0.7){$\Phi _{ij}^{ab}$} \end{picture}$

EOM-IP $\Psi (N)=R(-1) \Psi _0(N+1)$

$\begin{picture}(17,3)\thicklines \put( 0,1){\line(1,0){1}} \put( 0,2){\line(1,0){1}} \put( 0,3){\line(1,0){1}} \put(0.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(0.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put( 3,1){\line(1,0){1}} \put( 3,2){\line(1,0){1}} \put( 3,3){\line(1,0){1}} \put( 5,1){\line(1,0){1}} \put( 5,2){\line(1,0){1}} \put( 5,3){\line(1,0){1}} \put(3.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(3.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(5.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(5.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put( 8,1){\line(1,0){1}} \put( 8,2){\line(1,0){1}} \put( 8,3){\line(1,0){1}} \put(10,1){\line(1,0){1}} \put(10,2){\line(1,0){1}} \put(10,3){\line(1,0){1}} \put(12,1){\line(1,0){1}} \put(12,2){\line(1,0){1}} \put(12,3){\line(1,0){1}} \put(14,1){\line(1,0){1}} \put(14,2){\line(1,0){1}} \put(14,3){\line(1,0){1}} \put(8.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(8.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(10.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(10.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(12.48,0.8){\Large {\ensuremath{\downharpoonright }}} \put(12.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(12.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(14.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(14.48,1.8){\Large {\ensuremath{\downharpoonright }}} \put(14.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(3.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(4.5,-0.7){$\Phi _ i$} \put(8.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(11.5,-0.7){$\Phi _{ij}^ a$} \end{picture}$

EOM-EA $\Psi (N)=R(+1) \Psi _0(N-1)$

$\begin{picture}(17,3)\thicklines \put( 0,1){\line(1,0){1}} \put( 0,2){\line(1,0){1}} \put( 0,3){\line(1,0){1}} \put(0.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put( 3,1){\line(1,0){1}} \put( 3,2){\line(1,0){1}} \put( 3,3){\line(1,0){1}} \put( 5,1){\line(1,0){1}} \put( 5,2){\line(1,0){1}} \put( 5,3){\line(1,0){1}} \put(3.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(3.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(5.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(5.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put( 8,1){\line(1,0){1}} \put( 8,2){\line(1,0){1}} \put( 8,3){\line(1,0){1}} \put(10,1){\line(1,0){1}} \put(10,2){\line(1,0){1}} \put(10,3){\line(1,0){1}} \put(12,1){\line(1,0){1}} \put(12,2){\line(1,0){1}} \put(12,3){\line(1,0){1}} \put(14,1){\line(1,0){1}} \put(14,2){\line(1,0){1}} \put(14,3){\line(1,0){1}} \put(16,1){\line(1,0){1}} \put(16,2){\line(1,0){1}} \put(16,3){\line(1,0){1}} \put(8.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(8.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(10.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(10.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(10.48,2.8){\Large {\ensuremath{\downharpoonright }}} \put(12.48,0.8){\Large {\ensuremath{\downharpoonright }}} \put(12.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(12.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(14.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(14.48,1.8){\Large {\ensuremath{\downharpoonright }}} \put(14.18,2.8){\Large {\ensuremath{\upharpoonleft }}} \put(16.48,0.8){\Large {\ensuremath{\downharpoonright }}} \put(16.18,2.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(3.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(4.5,-0.7){$\Phi ^{a}$} \put(8.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(12.5,-0.7){$\Phi _ i^{ab}$} \end{picture}$

EOM-SF $\Psi (M_ s=0)=R(M_ s=-1) \Psi _0(M_ s=1)$

$\begin{picture}(17,2)\thicklines \put( 0,1){\line(1,0){1}} \put( 0,2){\line(1,0){1}} \put(0.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(0.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put( 3,1){\line(1,0){1}} \put( 3,2){\line(1,0){1}} \put( 5,1){\line(1,0){1}} \put( 5,2){\line(1,0){1}} \put( 7,1){\line(1,0){1}} \put( 7,2){\line(1,0){1}} \put( 9,1){\line(1,0){1}} \put( 9,2){\line(1,0){1}} \put(3.18,0.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(5.18,1.8){\Large {\ensuremath{\upharpoonleft }}\Large {\ensuremath{\downharpoonright }}} \put(7.18,0.8){\Large {\ensuremath{\upharpoonleft }}} \put(7.48,1.8){\Large {\ensuremath{\downharpoonright }}} \put(9.48,0.8){\Large {\ensuremath{\downharpoonright }}} \put(9.18,1.8){\Large {\ensuremath{\upharpoonleft }}} \put(3.0,0.4){$\underbrace{~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ }$} \put(4.5,-0.7){$\Phi _ i^{a}$} \end{picture}$

Figure 6.1: In the EOM formalism, target states $\Psi$ are described as excitations from a reference state $\Psi _0$ : $\Psi =R \Psi _0$ , where $R$ is a general excitation operator. Different EOM models are defined by choosing the reference and the form of the operator $R$ . In the EOM models for electronically excited states (EOM-EE, upper panel), the reference is the closed-shell ground state Hartree-Fock determinant, and the operator $R$ conserves the number of $\alpha$ and $\beta$ electrons. Note that two-configurational open-shell singlets can be correctly described by EOM-EE since both leading determinants appear as single electron excitations. The second and third panels present the EOM-IP/EA models. The reference states for EOM-IP/EA are determinants for $N+1$ / $N-1$ electron states, and the excitation operator $R$ is ionizing or electron-attaching, respectively. Note that both the EOM-IP and EOM-EA sets of determinants are spin-complete and balanced with respect to the target multi-configurational ground and excited states of doublet radicals. Finally, the EOM-SF method (the lowest panel) employs the hight-spin triplet state as a reference, and the operator $R$ includes spin-flip, i.e., does not conserve the number of $\alpha$ and $\beta$ electrons. All the determinants present in the target low-spin states appear as single excitations, which ensures their balanced treatment both in the limit of large and small HOMO-LUMO gaps.

6.7.2 EOM-XX-CCSD and CI Suite of Methods

Q-Chem features the most complete set of EOM-CCSD models [354] that enables accurate, robust, and efficient calculations of electronically excited states (EOM-EE-CCSD or EOM-EE-OD) [355, 356, 352, 353, 357]; ground and excited states of diradicals and triradicals (EOM-SF-CCSD and EOM-SF-OD [358, 357]); ionization potentials and electron attachment energies as well as problematic doublet radicals, cation or anion radicals, (EOM-IP/EA-CCSD) [359, 360, 361], as well as EOM-DIP-CCSD and EOM-2SF-CCSD. Conceptually, EOM is very similar to configuration interaction (CI): target EOM states are found by diagonalizing the similarity transformed Hamiltonian $\bar H= e^{-T}H e^{T}$ ,

$\begin{equation} \bar H R = ER, \end{equation}$

(6.30)

where $T$ and $R$ are general excitation operators with respect to the reference determinant $|\Phi _0\ensuremath{\rangle }$ . In the EOM-CCSD models, $T$ and $R$ are truncated at single and double excitations, and the amplitudes $T$ satisfy the CC equations for the reference state $|\Phi _0\ensuremath{\rangle }$ :

	$\displaystyle \ensuremath{\langle }\Phi _ i^ a\|\bar H \| \Phi _0\ensuremath{\rangle }$	$\displaystyle =$	$\displaystyle 0$		(6.31)
	$\displaystyle \ensuremath{\langle }\Phi _{ij}^{ab}\|\bar H \| \Phi _0\ensuremath{\rangle }$	$\displaystyle =$	$\displaystyle 0$		(6.32)

The computational scaling of EOM-CCSD and CISD methods is identical, i.e., ${\cal {O}}({N^6})$ , however EOM-CCSD is numerically superior to CISD because correlation effects are “folded in” in the transformed Hamiltonian, and because EOM-CCSD is rigorously size-intensive.

By combining different types of excitation operators and references $|\Phi _0\ensuremath{\rangle }$ , different groups of target states can be accessed as explained in Fig. 6.1. For example, electronically excited states can be described when the reference $|\Phi _0\ensuremath{\rangle }$ corresponds to the ground state wave function, and operators $R$ conserve the number of electrons and a total spin [352]. In the ionized/electron attached EOM models [360, 361], operators $R$ are not electron conserving (i.e., include different number of creation and annihilation operators)—these models can accurately treat ground and excited states of doublet radicals and some other open-shell systems. For example, singly ionized EOM methods, i.e., EOM-IP-CCSD and EOM-EA-CCSD, have proven very useful for doublet radicals whose theoretical treatment is often plagued by symmetry breaking. Finally, the EOM-SF method [358, 357] in which the excitation operators include spin-flip allows one to access diradicals, triradicals, and bond-breaking.

Q-Chem features EOM-EE/SF/IP/EA-CCSD methods for both closed and open-shell references (RHF/UHF/ROHF), including frozen core/virtual options. All EOM models take full advantage of molecular point group symmetry. Analytic gradients are available for RHF and UHF references, for the full orbital space, and with frozen core/virtual orbitals [362]. Properties calculations (permanent and transition dipole moments, $\ensuremath{\langle }S^2\ensuremath{\rangle }$ , $\ensuremath{\langle }R^2 \ensuremath{\rangle }$ , etc.) are also available. The current implementation of the EOM-XX-CCSD methods enables calculations of medium-size molecules, e.g., up to 15–20 heavy atoms. Using RI approximation 5.7.5 or Cholesky decomposition 5.7.6 helps to reduce integral transformation time and disk usage enabling calculations on much larger systems.

The CCMAN module of Q-Chem includes two implementations of EOM-IP-CCSD. The proper implementation [363] is used by default is more efficient and robust. The EOM_FAKE_IPEA keyword invokes is a pilot implementation in which EOM-IP-CCSD calculation is set up by adding a very diffuse orbital to a requested basis set, and by solving EOM-EE-CCSD equations for the target states that include excitations of an electron to this diffuse orbital. The implementation of EOM-EA-CCSD in CCMAN also uses this trick. Fake IP/EA calculations are only recommended for Dyson orbital calculations and debug purposes.

CCMAN2 features proper implementations of EOM-IP and EOM-EA (including Dyson orbitals).

The computational cost of EOM-IP calculations can be considerably reduced (with negligible decline in accuracy) by truncating virtual orbital space using FNO scheme (see Section 6.7.7).

Finally, a more economical CI variant of EOM-IP-CCSD, IP-CISD is also available. This is an N $^5$ approximation of IP-CCSD, and is recommended for geometry optimizations of problematic doublet states [364].

EOM and CI methods are handled by the CCMAN/CCMAN2 modules.

6.7.3 Spin-Flip Methods for Di- and Triradicals

The spin-flip method [358, 316, 365] addresses the bond-breaking problem associated with a single-determinant description of the wavefunction. Both closed and open shell singlet states are described within a single reference as spin-flipping, (e.g., $\alpha \to \beta$ excitations from the triplet reference state, for which both dynamical and non-dynamical correlation effects are smaller than for the corresponding singlet state. This is because the exchange hole, which arises from the Pauli exclusion between same-spin electrons, partially compensates for the poor description of the coulomb hole by the mean-field Hartree-Fock model. Furthermore, because two $\alpha$ electrons cannot form a bond, no bond breaking occurs as the internuclear distance is stretched, and the triplet wavefunction remains essentially single-reference in character. The spin-flip approach has also proved useful in the description of di- and tri-radicals as well as some problematic doublet states.

The spin-flip method is available for the CIS, CIS(D), CISD, CISDT, OD, CCSD, and EOM-(2,3) levels of theory and the spin complete SF-XCIS (see Section 6.2.4). An N $^7$ non-iterative triples corrections are also available. For the OD and CCSD models, the following non-relaxed properties are also available: dipoles, transition dipoles, eigenvalues of the spin-squared operator ( $\langle S^2\rangle$ ), and densities. Analytic gradients are also for SF-CIS and EOM-SF-CCSD methods. To invoke a spin-flip calculation the EOM_SF_STATES $rem should be used, along with the associated $rem settings for the chosen level of correlation (CORRELATION, and, optionally, EOM_CORR). Note that the high multiplicity triplet or quartet reference states should be used.

Several double SF methods have also been implemented [366]. To invoke these methods, use
EOM_DSF_STATES.

6.7.4 EOM-DIP-CCSD

Double-ionization potential (DIP) is another non-electron-conserving variant of EOM-CCSD [367, 368, 369]. In DIP, target states are reached by detaching two electrons from the reference state:

$\begin{equation} \Psi _ k = R_{N-2} \Psi _0 (N+2), \protect \label{DIP:WF} \end{equation}$

(6.33)

and the excitation operator $R$ has the following form:

$\displaystyle \label{eq:R-IP} R$	$\displaystyle =$	$\displaystyle R_1 + R_2,$	(6.34)
$\displaystyle R_1$	$\displaystyle =$	$\displaystyle 1/2 \sum _{ij} r_{ij} ji,$	(6.35)
$\displaystyle R_2$	$\displaystyle =$	$\displaystyle 1/6 \sum _{ijka} r_{ijk}^{a} a^\dag kji.$	(6.36)

As a reference state in the EOM-DIP calculations one usually takes a well-behaved closed-shell state. EOM-DIP is a useful tool for describing molecules with electronic degeneracies of the type “ $2n-2$ electrons on $n$ degenerate orbitals”. The simplest examples of such systems are diradicals with two-electrons-on-two-orbitals pattern. Moreover, DIP is a preferred method for four-electrons-on-three-orbitals wavefunctions.

Accuracy of the EOM-DIP-CCSD method is similar to accuracy of other EOM-CCSD models, i.e., 0.1–0.3 eV. The scaling of EOM-DIP-CCSD is ${\cal {O}}({N^6})$ , analogous to that of other EOM-CCSD methods. However, its computational cost is less compared to, e.g., EOM-EE-CCSD, and it increases more slowly with the basis set size.

An EOM-DIP calculation is invoked by using EOM_DIP_STATES, or EOM_DIP_SINGLETS and
EOM_DIP_TRIPLETS.

6.7.5 EOM-CC Calculations of Metastable States: Super-Excited Electronic States, Temporary Anions, and Core-Ionized States

While conventional coupled-cluster and equation-of-motion methods allow one to tackle electronic structure ranging from well-behaved closed shell molecules to various open-shell and electronically excited species [354], metastable electronic states, so-called resonances, present a difficult case for theory. By using complex scaling and complex absorbing potential techniques, we extended these powerful methods to describe autoionizing states, such as transient anions, highly excited electronic states, and core-ionized species [370, 371]. In addition, users can employ stabilization techniques using charged sphere and scaled atomic charges options [369]. These methods are only available within CCMAN2. The complex CC/EOM code is engaged by COMPLEX_CCMAN; the specific parameters should be specified in the $complex_ccman section.

COMPLEX_CCMAN

Requests complex-scaled or CAP-augmented CC/EOM calculations.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE

Engage complex CC/EOM code.

RECOMMENDATION:

Not available in CCMAN. Need to specify CAP strength or complex-scaling parameter in $complex_ccman section.

The $complex_ccman section is used to specify the details of the complex-scaled/CAP calculations, as illustrated below. If user specifies CS_THETA, complex scaling calculation is performed.

$complex_ccman
CS_THETA  10   Complex-scaling parameter theta=0.01, r->r exp(-i*theta)
CS_ALPHA  10   Real part of the scaling parameter alpha=0.01, r->alpha r exp(-itheta)
$end

Alternatively, for CAP calculations, the CAP parameters need to be specified.

$complex_ccman
CAP_ETA 1000  CAP strength in 10-5 a.u. (0.01)
CAP_X 2760  CAP onset along X in 10^-3 bohr (2.76 bohr)
CAP_Y 2760  CAP onset along Y in 10^-3 bohr (2.76 bohr)
CAP_Z 4880  CAP onset along Z in 10^-3 bohr (4.88 bohr)
CAP_TYPE 1  Use cuboid cap (CAP_TYPE=0 will use spherical CAP) 
$end

CS_THETA is specified in radian $\times$ 10 $^{-3}$ . CS_ALPHA, CAP_X/Y/Z are specified in a.u. $\times$ 10 $^{-3}$ , i.e., CS_THETA=10 means $\theta$ =0.01; CAP_ETA is specified in a.u. $\times$ 10 $^{-5}$ . When CAP calculations are performed, CC_EOM_PROP=TRUE by default; this is necessary for calculating first-order deperturbative correction.

Advanced users may find the following options useful. Several ways of conducing complex calculations are possible, i.e., complex scaling/CAPs can be either engaged at all levels (HF, CCSD, EOM), or not. By default, if COMPLEX_CCMAN is specified, the EOM calculations are conducted using complex code. Other parameters are set up as follows:

$complex_ccman
CS_HF=true 
CS_CCSD=true 
$end

Alternatively, the user can disable complex HF. These options are experimental and should only be used by advanced users. For CAP-EOM-CC, only CS_HF=TRUE and CS_CCSD=TRUE is implemented. Finally, only RHF references are supported in complex-scaled and CAP-augmented calculations.

6.7.6 Charge Stabilization for EOM-DIP and Other Methods

Unfortunately, the performance of EOM-DIP deteriorates when the reference state is unstable with respect to electron-detachment [368, 369], which is usually the case for dianion references employed to describe neutral diradicals by EOM-DIP. Similar problems are encountered by all excited-state methods when dealing with excited states lying above ionization or electron-detachment thresholds.

To remedy this problem, one can employ charge stabilization methods, as described in Refs. [368, 369]. In this approach (which can also be used with any other electronic structure method implemented in Q-Chem), an additional Coulomb potential is introduced to stabilize unstable wave functions. The following keywords invoke stabilization potentials: SCALE_NUCLEAR_CHARGE and ADD_CHARGED_CAGE. In the former case, the potential is generated by increasing nuclear charges by a specified amount. In the latter, the potential is generated by a cage built out of point charges comprising the molecule. There are two cages available: dodecahedral and spherical. The shape, radius, number of points, and the total charge of the cage are set by the user.

Note: A perturbative correction estimating the effect of the external Coulomb potential on EOM energy will be computed when target state densities are calculated, e.g., when CC_EOM_PROP is set to TRUE.

Note: Charge stabilization techniques can be used with other methods such as EOM-EE, CIS, and TDDFT to improve the description of resonances. It can also be employed to describe metastable ground states.

6.7.7 Frozen Natural Orbitals in CC and IP-CC Calculations

Large computational savings are possible if the virtual space is truncated using the frozen natural orbital (FNO) approach (see Section 5.10). Extension of the FNO approach to ionized states within EOM-CC formalism was recently introduced and benchmarked [278]. In addition to ground-state coupled-cluster calculations, FNOs can also be used in EOM-IP-CCSD, EOM-IP-CCSD(dT/fT) and EOM-IP-CC(2,3). In IP-CC the FNOs are computed for the reference (neutral) state and then are used to describe several target (ionized) states of interest. Different truncation scheme are described in Section 5.10.

6.7.8 EOM-CC Jobs: Controlling Guess Formation and Iterative Diagonalizers

An EOM-CC eigen problem is solved by an iterative diagonalization procedure that avoids full diagonalization and only looks for several eigen-states, as specified by the XX_STATES keywords.

The default procedure is based on the modified Davidson diagonalization algorithm, as explained in Ref. [357]. In addition to several keywords that control the convergence of algorithm, memory usage, and fine details of its execution, there are several important keywords that allow user to specify how the target state selection will be performed.

By default, the diagonalization looks for several lowest eigenstates, as specified by XX_STATES. The guess vectors are generated as singly excited determinants selected by using the Koopmans theorem; the number of guess vectors is equal to the number of target states. If necessary, the user can increase the number of singly excited guess vectors (EOM_NGUESS_SINGLES), include doubly excited guess vectors (EOM_NGUESS_DOUBLES), and even preconverge them (EOM_PRECONV_SINGLES and EOM_PRECONV_DOUBLES).

If a state of a particular character is desired (e.g., HOMO $\rightarrow$ LUMO+10 excitation or HOMO-10 ionization), the user can specify this by using EOM_USER_GUESS keyword and $eom_user_guess section, as illustrated by an example below. The algorithm will attempt to find an eigenstate that has the maximum overlap with this guess vector. Note that the point group symmetry of user defined guess vectors should be consistent with the symmetry of the transition specified by XX_STATES. The multiplicity of the state is determined as in the regular calculations, by using the EOM_XX_SINGLETS and EOM_EE_TRIPLETS keywords. This option is useful for looking for high-lying states such as core-ionized states.

The examples below illustrate how to use user-specified guess in EOM calculations:

$eom_user_guess
4  Calculate excited state corresponding to 4(OCC)->5(VIRT) transition. 
5
$end

In IP/EA calculations, only one orbital is specified:

$eom_user_guess
4  
$end

If IP_STATES is specified, this will invoke calculation of the EOM-IP state corresponding to the ionization from 4th occupied MO. If EA_STATES is requested, then EOM-EA equations will be solved for a root corresponding to electron-attachment to the 4th virtual MO.

The symmetries of the MOs should be consistent with the EE_STATES value, which should only request one state in the correct irrep. The orbitals are ordered by energy as printed at the beginning of the CCMAN2 output.

Alternatively, the user can specify an energy shift by EOM_SHIFT. In this case, the solver looks for the XX_STATES eigenstates that are closest to this energy; the guess vectors are generated accordingly, using the Koopmans theorem. This option is useful when highly excited states (i.e., interior eigenstates) are desired.

Finally, a new diagonalization technique, so-called GPLMR (or GPLHR), is available[372]. This solver is engaged by EOM_GPLMR keyword and is only available in CCMAN2. The GPLMR method usually converges in fewer iterations relative to Davidson and can use less memory, but it performs more floating point operations. Similarly to the Davidson procedure, GPLMR can be applied to look for either the lowest eigenstates, or for an interior set, as specified by EOM_SHIFT. The convergence is controlled by the same keywords as in the Davidson algorithm. One additional keyword controlling the maximum subspace size in the GPLMR solver is EOM_GPLMR_MSUBSIZE.

6.7.9 Equation-of-Motion Coupled-Cluster Job Control

It is important to ensure there are sufficient resources available for the necessary integral calculations and transformations. For CCMAN/CCMAN2 algorithms, these resources are controlled using the $rem variables CC_MEMORY, MEM_STATIC and MEM_TOTAL (see Section 5.13).

Equation-of-motion methods require a coupled-cluster reference state, which is computed when METHOD is set to EOMCCSD or EOMOD. In addition, a specification of the number of target states is required through XX_STATES (XX designates the type of the target states, e.g., EE, SF, IP, EA, DIP, DSF, etc.). Users must be aware of the point group symmetry of the system being studied and also the symmetry of the initial and target states of interest, as well as symmetry of transition. It is possible to turn off the use of symmetry by CC_SYMMETRY. If set to FALSE the molecule will be treated as having $C_1$ symmetry and all states will be of $A$ symmetry.

Note: In finite-difference calculations, the symmetry is turned off automatically, and the user must ensure that XX_STATES is adjusted accordingly.

Note: Mixing different EOM models in a single calculation is only allowed in Dyson orbitals calculations.

By default, the level of correlation of the EOM part of the wavefunction (i.e., maximum excitation level in the EOM operators $R$ ) is set to match CORRELATION, however, one can mix different correlation levels for the reference and EOM states by using EOM_CORR. To request a CI calculation, set CORRELATION=CI and select type of CI expansion by EOM_CORR. The table below shows default and allowed CORRELATION and EOM_CORR combinations.

CORRELATION	Default	Allowed	Target states	CCMAN/CCMAN2
	EOM_CORR	EOM_CORR
CI	none	CIS, CIS(D)	EE,SF	y/n
		CISD	EE,SF,IP	y/n
		SDT, DT	EE,SF,DSF	y/n
CIS(D)	CIS(D)	N/A	EE,SF	y/n
CCSD, OD	CISD		EE,SF,IP,EA,DIP	y/y
		SD(dT),SD(fT)	EE,SF, fake IP/EA	y/n
		SD(dT),SD(fT), SD(sT)	IP	y/n
		SDT, DT	EE,SF,IP,EA,DIP,DSF	y/n

Table 6.1: Default and allowed CORRELATION and EOM_CORR combinations as well as valid target state types. The last column shows if a method is available in CCMAN or CCMAN2.

The table below shows the correct combinations of CORRELATION and EOM_CORR for standard EOM and CI models.

Method	CORRELATION	EOM_CORR	Target states selection
CIS	CI	CIS	EE_STATES
			EE_SNGLETS,EE_TRIPLETS
SF-CIS	CI	CIS	SF_STATES
CIS(D)	CI	CIS(D)	EE_STATES
			EE_SNGLETS,EE_TRIPLETS
SF-CIS(D)	CI	CIS(D)	SF_STATES
CISD	CI	CISD	EE_STATES
			EE_SNGLETS,EE_TRIPLETS
SF-CISD	CI	CISD	SF_STATES
IP-CISD	CI	CISD	IP_STATES
CISDT	CI	SDT	EE_STATES
			EE_SNGLETS,EE_TRIPLETS
SF-CISDT	CI	SDT or DT	SF_STATES
EOM-EE-CCSD	CCSD		EE_STATES
			EE_SNGLETS,EE_TRIPLETS
EOM-SF-CCSD	CCSD		SF_STATES
EOM-IP-CCSD	CCSD		IP_STATES
EOM-EA-CCSD	CCSD		EA_STATES
EOM-DIP-CCSD	CCSD		DIP_STATES
			DIP_SNGLETS,DIP_TRIPLETS
EOM-2SF-CCSD	CCSD	SDT or DT	DSF_STATES
EOM-EE-(2,3)	CCSD	SDT	EE_STATES
			EE_SNGLETS,EE_TRIPLETS
EOM-SF-(2,3)	CCSD	SDT	SF_STATES
EOM-IP-(2,3)	CCSD	SDT	IP_STATES
EOM-SF-CCSD(dT)	CCSD	SD(dT)	SF_STATES
EOM-SF-CCSD(fT)	CCSD	SD(fT)	SF_STATES
EOM-IP-CCSD(dT)	CCSD	SD(dT)	IP_STATES
EOM-IP-CCSD(fT)	CCSD	SD(fT)	IP_STATES
EOM-IP-CCSD(sT)	CCSD	SD(sT)	IP_STATES

Table 6.2: Commonly used EOM and CI models. ’SINGLETS’ and ’TRIPLETS’ are only available for closed-shell references.

The most relevant EOM-CC input options follow.

EE_STATES

Sets the number of excited state roots to find. For closed-shell reference, defaults into EE_SINGLETS. For open-shell references, specifies all low-lying states.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any excited states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ excited states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EE_SINGLETS

Sets the number of singlet excited state roots to find. Valid only for closed-shell references.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any excited states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ excited states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EE_TRIPLETS

Sets the number of triplet excited state roots to find. Valid only for closed-shell references.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any excited states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ excited states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

SF_STATES

Sets the number of spin-flip target states roots to find.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any excited states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ SF states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

DSF_STATES

Sets the number of doubly spin-flipped target states roots to find.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any DSF states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ doubly spin-flipped states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

IP_STATES

Sets the number of ionized target states roots to find. By default, $\beta$ electron will be removed (see EOM_IP_BETA).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any IP states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ ionized states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EOM_IP_ALPHA

Sets the number of ionized target states derived by removing $\alpha$ electron (M $_ s$ = $-{{1}\over {2}}$ ).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any IP/ $\alpha$ states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ ionized states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EOM_IP_BETA

Sets the number of ionized target states derived by removing $\beta$ electron (M $_ s$ = ${{1}\over {2}}$ , default for EOM-IP).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any IP/ $\beta$ states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ ionized states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EA_STATES

Sets the number of attached target states roots to find. By default, $\alpha$ electron will be attached (see EOM_EA_ALPHA).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any EA states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ EA states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EOM_EA_ALPHA

Sets the number of attached target states derived by attaching $\alpha$ electron (M $_ s$ = ${{1}\over {2}}$ , default in EOM-EA).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any EA states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ EA states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

EOM_EA_BETA

Sets the number of attached target states derived by attaching $\beta$ electron (M $_ s$ = $-{{1}\over {2}}$ , EA-SF).

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any EA states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ EA states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

DIP_STATES

Sets the number of DIP roots to find. For closed-shell reference, defaults into DIP_SINGLETS. For open-shell references, specifies all low-lying states.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any DIP states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ DIP states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

DIP_SINGLETS

Sets the number of singlet DIP roots to find. Valid only for closed-shell references.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any singlet DIP states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ DIP singlet states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

DIP_TRIPLETS

Sets the number of triplet DIP roots to find. Valid only for closed-shell references.

TYPE:

INTEGER/INTEGER ARRAY

DEFAULT:

0

Do not look for any DIP triplet states.

OPTIONS:

$[i,j,k\ldots ]$

Find $i$ DIP triplet states in the first irrep, $j$ states in the second irrep etc.

RECOMMENDATION:

None

Note: It is a symmetry of a transition rather than that of a target state which is specified in excited state calculations. The symmetry of the target state is a product of the symmetry of the reference state and the transition. For closed-shell molecules, the former is fully symmetric and the symmetry of the target state is the same as that of transition, however, for open-shell references this is not so.

Note: For the XX_STATES options, Q-Chem will increase the number of roots if it suspects degeneracy, or change it to a smaller value, if it cannot generate enough guess vectors to start the calculations.

EOM_FAKE_IPEA

If TRUE, calculates fake EOM-IP or EOM-EA energies and properties using the diffuse orbital trick. Default for EOM-EA and Dyson orbital calculations in CCMAN.

TYPE:

LOGICAL

DEFAULT:

FALSE (use proper EOM-IP code)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

None. This feature only works for CCMAN.

Note: When EOM_FAKE_IPEA is set to TRUE, it can change the convergence of Hartree-Fock iterations compared to the same job without EOM_FAKE_IPEA, because a very diffuse basis function is added to a center of symmetry before the Hartree-Fock iterations start. For the same reason, BASIS2 keyword is incompatible with EOM_FAKE_IPEA. In order to read Hartree-Fock guess from a previous job, you must specify EOM_FAKE_IPEA (even if you do not request for any correlation or excited states) in that previous job. Currently, the second moments of electron density and Mulliken charges and spin densities are incorrect for the EOM-IP/EA-CCSD target states.

EOM_USER_GUESS

Specifies if user-defined guess will be used in EOM calculations.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE

Solve for a state that has maximum overlap with a trans-n specified in $eom_user_guess.

RECOMMENDATION:

The orbitals are ordered by energy, as printed in the beginning of the CCMAN2 output. Not available in CCMAN.

EOM_SHIFT

Specifies energy shift in EOM calculations.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$

corresponds to $n\cdot 10^{-3}$ hartree shift (i.e., 11000 = 11 hartree); solve for eigenstates around this value.

RECOMMENDATION:

Not available in CCMAN.

EOM_GPLMR

Specifies whether to engage GPLMR solver in EOM calculations.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE

Use GPLMR.

RECOMMENDATION:

Not available in CCMAN.

EOM_NGUESS_DOUBLES

Specifies number of excited state guess vectors which are double excitations.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$

Include $n$ guess vectors that are double excitations

RECOMMENDATION:

This should be set to the expected number of doubly excited states (see also EOM_PRECONV_DOUBLES), otherwise they may not be found.

EOM_NGUESS_SINGLES

Specifies number of excited state guess vectors that are single excitations.

TYPE:

INTEGER

DEFAULT:

Equal to the number of excited states requested

OPTIONS:

$n$

Include $n$ guess vectors that are single excitations

RECOMMENDATION:

Should be greater or equal than the number of excited states requested.

EOM_PRECONV_SINGLES

When not zero, singly excited vectors are converged prior to a full excited states calculation. Sets the maximum number of iterations for pre-converging procedure

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0

do not pre-converge

N

perform N Davidson iterations pre-converging singles.

RECOMMENDATION:

Sometimes helps with problematic convergence.

EOM_PRECONV_DOUBLES

When not zero, doubly excited vectors are converged prior to a full excited states calculation. Sets the maximum number of iterations for pre-converging procedure

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0

do not pre-converge

N

perform N Davidson iterations pre-converging doubles.

RECOMMENDATION:

Occasionally necessary to ensure a doubly excited state is found. Also used in DSF calculations instead of EOM_PRECONV_SINGLES

EOM_PRECONV_SD

When not zero, EOM vectors are pre-converged prior to a full excited states calculation. Sets the maximum number of iterations for pre-converging procedure

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

0

do not pre-converge

N

perform N Davidson iterations pre-converging singles and doubles.

RECOMMENDATION:

Occasionally necessary to ensure that all low-lying states are found. Also, very useful in EOM(2,3) calculations.

None

EOM_DAVIDSON_CONVERGENCE

Convergence criterion for the RMS residuals of excited state vectors

TYPE:

INTEGER

DEFAULT:

5

Corresponding to $10^{-5}$

OPTIONS:

$n$

Corresponding to $10^{-n}$ convergence criterion

RECOMMENDATION:

Use default. Should normally be set to the same value as EOM_DAVIDSON_THRESHOLD.

EOM_DAVIDSON_THRESHOLD

Specifies threshold for including a new expansion vector in the iterative Davidson diagonalization. Their norm must be above this threshold.

TYPE:

INTEGER

DEFAULT:

00105

Corresponding to 0.00001

OPTIONS:

$abcde$

Integer code is mapped to $abc\times 10^{-de}$

RECOMMENDATION:

Use default unless converge problems are encountered. Should normally be set to the same values as EOM_DAVIDSON_CONVERGENCE, if convergence problems arise try setting to a value less than EOM_DAVIDSON_CONVERGENCE.

EOM_DAVIDSON_MAXVECTORS

Specifies maximum number of vectors in the subspace for the Davidson diagonalization.

TYPE:

INTEGER

DEFAULT:

60

OPTIONS:

$n$

Up to $n$ vectors per root before the subspace is reset

RECOMMENDATION:

Larger values increase disk storage but accelerate and stabilize convergence.

EOM_DAVIDSON_MAX_ITER

Maximum number of iteration allowed for Davidson diagonalization procedure.

TYPE:

INTEGER

DEFAULT:

30

OPTIONS:

$n$

User-defined number of iterations

RECOMMENDATION:

Default is usually sufficient

EOM_GPLMR_MSUBSIZE

Specifies the number of Krylov-space residuals in GPLMR.

TYPE:

INTEGER

DEFAULT:

3

OPTIONS:

$n$

Generate $n$ residuals at each iteration.

RECOMMENDATION:

Use default. The convergence is faster for larger $n$ , but the memory usage and the overall cost will increase.

EOM_IPEA_FILTER

If TRUE, filters the EOM-IP/EA amplitudes obtained using the diffuse orbital implementation (see EOM_FAKE_IPEA). Helps with convergence.

TYPE:

LOGICAL

DEFAULT:

FALSE (EOM-IP or EOM-EA amplitudes will not be filtered)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

None

CC_FNO_THRESH

Initialize the FNO truncation and sets the threshold to be used for both cutoffs (OCCT and POVO)

TYPE:

INTEGER

DEFAULT:

None

OPTIONS:

range

0000-10000

$abcd$

Corresponding to $ab.cd$ %

RECOMMENDATION:

None

CC_FNO_USEPOP

Selection of the truncation scheme

TYPE:

INTEGER

DEFAULT:

1

OCCT

OPTIONS:

0

POVO

RECOMMENDATION:

None

SCALE_NUCLEAR_CHARGE

Scales charge of each nuclei by a certain value. The nuclear repulsion energy is calculated for the unscaled nuclear charges.

TYPE:

INTEGER

DEFAULT:

0 no scaling.

OPTIONS:

n a total positive charge of (1+n/100)e is added to the molecule.

RECOMMENDATION:

NONE

ADD_CHARGED_CAGE

Add a point charge cage of a given radius and total charge.

TYPE:

INTEGER

DEFAULT:

0 no cage.

OPTIONS:

0 no cage.

1 dodecahedral cage.

2 spherical cage.

RECOMMENDATION:

Spherical cage is expected to yield more accurate results, especially for small radii.

CAGE_RADIUS

Defines radius of the charged cage.

TYPE:

INTEGER

DEFAULT:

225

OPTIONS:

n radius is n/100 .

RECOMMENDATION:

None

CAGE_POINTS

Defines number of point charges for the spherical cage.

TYPE:

INTEGER

DEFAULT:

100

OPTIONS:

n n point charges are used.

RECOMMENDATION:

None

CAGE_CHARGE

Defines the total charge of the cage.

TYPE:

INTEGER

DEFAULT:

400 Add a cage charged +4e.

OPTIONS:

n total charge of the cage is n/100 a.u.

RECOMMENDATION:

None

6.7.10 Examples

Example 6.130 EOM-EE-OD and EOM-EE-CCSD calculations of the singlet excited states of formaldehyde

$molecule
0 1
O
C,1,R1
H,2,R2,1,A
H,2,R2,1,A,3,180.

R1=1.4
R2=1.0
A=120.
$end

$rem
method      eom-od
basis       6-31+g
ee_states   [2,2,2,2]
$end

@@@
$molecule
read
$end

$rem
method        eom-ccsd
basis         6-31+g
ee_singlets   [2,2,2,2]
ee_triplets   [2,2,2,2]
$end

Example 6.131 EOM-EE-CCSD calculations of the singlet excited states of PYP using Cholesky decomposition

$molecule
0 1
...too long to enter...
$end

$rem
method          eom-ccsd
basis           aug-cc-pVDZ
purecart        1112
n_frozen_core   fc
cc_t_conv       4
cc_e_conv       6
cholesky_tol    2    using CD/1e-2 threshold
ee_singlets     [2,2]
$end

Example 6.132 EOM-SF-CCSD calculations for methylene from high-spin $^3$ B $_2$ reference

$molecule
0 3
C
H  1 rCH
H  1 rCH 2 aHCH

rCH    = 1.1167 
aHCH   = 102.07
$end

$rem
method          eom-ccsd
basis           6-31G*
scf_guess       core
sf_states       [2,0,0,2]   Two singlet A1 states and singlet and triplet B2 states
$end

Example 6.133 EOM-IP-CCSD calculations for NO $_3$ using closed-shell anion reference

$molecule
-1 1
N
O 1 r1
O 1 r2   2 A2
O 1 r2   2 A2    3 180.0

r1    = 1.237
r2    = 1.237
A2    = 120.00
$end

$rem
method          eom-ccsd
basis           6-31G*
ip_states       [1,1,2,1]    ground and excited  states of the radical
$end

Example 6.134 EOM-IP-CCSD calculation using FNO with OCCT=99%.

$molecule
  0 1
  O
  H 1 1.0
  H 1 1.0 2 100.
$end

$rem
method          eom-ccsd
basis           6-311+G(2df,2pd)
ip_states       [1,0,1,1]
cc_fno_thresh   9900        99% of the total natural population recovered
$end

Example 6.135 DSF-CIDT calculation of methylene starting with quintet reference

$molecule
0 5
C 
H 1 CH
H 1 CH 2 HCH

CH  = 1.07
HCH = 111.0
$end

$rem
   method             cisdt
   basis              6-31G
   dsf_states         [0,2,2,0]
   eom_nguess_singles  0   
   eom_nguess_doubles  2
$end

Example 6.136 EOM-EA-CCSD job for cyano radical. We first do Hartree-Fock calculation for the cation in the basis set with one extremely diffuse orbital (EOM_FAKE_IPEA) and use these orbitals in the second job. We need make sure that the diffuse orbital is occupied using the OCCUPIED keyword. No SCF iterations are performed as the diffuse electron and the molecular core are uncoupled. The attached states show up as “excited” states in which electron is promoted from the diffuse orbital to the molecular ones.

$molecule
+1 1
C
N 1 bond

bond   1.1718
$end

$rem
jobtype       sp
method        hf
basis         6-311+G*
purecart      111
scf_convergence 8
eom_fake_ipea true
$end

@@@
$molecule
0 2
C
N 1 bond

bond   1.1718
$end

$rem
jobtype        sp
basis          6-311+G*
purecart       111
scf_guess      read
max_scf_cycles 0
method         eom-ccsd
cc_dov_thresh  2501   use threshold for CC iterations with problematic convergence
ea_states      [2,0,0,0]
eom_fake_ipea  true
$end

$occupied
1 2 3 4 5 6 14
1 2 3 4 5 6
$end

Example 6.137 EOM-DIP-CCSD calculation of electronic states in methylene using charged cage stabilization method.

$molecule
-2 1
C   0.000000     0.000000     0.106788
H  -0.989216     0.000000    -0.320363
H   0.989216     0.000000    -0.320363
$end

$rem
jobtype = sp
basis = 6-311g(d,p)
scf_algorithm = diis_gdm
symmetry = false
method = eom-ccsd
cc_symmetry = false
dip_singlets = [1]         ! Compute one EOM-DIP singlet state
dip_triplets = [1]         ! Compute one EOM-DIP triplet state
eom_davidson_convergence = 5
cc_eom_prop = true         ! Compute excited state properties
add_charged_cage = 2           ! Install a charged sphere around the molecule
cage_radius = 225  ! Radius = 2.25 A
cage_charge = 500  ! Charge = +5 a.u.
cage_points = 100  ! Place 100 point charges
cc_memory = 256            ! Use 256Mb of memory, increase for larger jobs
$end

Example 6.138 EOM-EE-CCSD calculation of excited states in NO $^-$ using scaled nuclear charge stabilization method.

$molecule
-1 1
 N  -1.08735    0.0000    0.0000
 O   1.08735    0.0000    0.0000
$end


$rem
jobtype = sp
input_bohr = true
basis = 6-31g
symmetry = false
cc_symmetry = false
method = eom-ccsd
ee_singlets = [2]          ! Compute two EOM-EE singlet excited states
ee_triplets = [2]          ! Compute two EOM-EE triplet excited states
cc_ref_prop = true         ! Compute ground state properties
cc_eom_prop = true         ! Compute excited state properties
cc_memory = 256            ! Use 256Mb of memory, increase for larger jobs
scale_nuclear_charge = 180 ! Adds +1.80e charge to the molecule
$end

Example 6.139 EOM-EE-CCSD calculation for phenol with user-specified guess requesting the EE transition from the occupied orbital number 24 (3 A") to the virtual orbital number 2 (23 A’)

$molecule
0 1
C         0.935445   -0.023376    0.000000
C         0.262495    1.197399    0.000000
C        -1.130915    1.215736    0.000000
C        -1.854154    0.026814    0.000000
C        -1.168805   -1.188579    0.000000
C         0.220600   -1.220808    0.000000
O         2.298632   -0.108788    0.000000
H         2.681798    0.773704    0.000000
H         0.823779    2.130309    0.000000
H        -1.650336    2.170478    0.000000
H        -2.939976    0.044987    0.000000
H        -1.722580   -2.123864    0.000000
H         0.768011   -2.158602    0.000000
$end

$rem
jobtype            SP           single point
METHOD             EOM-CCSD
BASIS              6-31+G(d,p)
cc_memory          3000		ccman2 memory
mem_static         250
CC_T_CONV          4		T-amplitudes convergence threshold
CC_E_CONV          6		Energy convergence threshold
ee_states          [0,1]        Calculate 1 A" states
EOM_DAVIDSON_CONVERGENCE    5   Convergence threshold for the Davidson procedure
EOM_DAVIDSON_THRESHOLD      5   Threshold for inclusion of new vectors to the subspace
eom_user_guess     true         Use user guess from $eom_user_guess section
$end

$eom_user_guess
24   Transition from the occupied orbital number 24(3 A") 
2    to the virtual orbital number 2 (23 A') 
$end

Example 6.140 Complex-scaled EOM-EE calculation for He. All roots of Ag symmetry are computed (full diagonalization)

$molecule
0 1
HE  0    0  0.0
$END

$rem
complex_ccman      1             engage complex_ccman
jobtype            SP            single point
METHOD             EOM-CCSD
BASIS              gen           use general basis
PURECART           1111
EE_SINGLETS        [2000,0,0,0,0,0,0,0]  compute all excitation energies of Ag symmetry
EOM_DAVIDSON_CONV      5
EOM_DAVIDSON_THRESH    5
EOM_NGUESS_SINGLES    2000       Number of guess singles
EOM_NGUESS_DOUBLES    2000       Number of guess doubles
cc_memory 5000
mem_total 3000
$end

$complex_ccman
CS_HF 1                         Use complex HF
CS_ALPHA 1000                   Set alpha equal 1
CS_THETA 300                    Set theta (angle) equals 0.3 (radian)
$end

$basis
He   0
S    4    1.000000
   2.34000000E+02    2.58700000E-03 
   3.51600000E+01    1.95330000E-02 
   7.98900000E+00    9.09980000E-02 
   2.21200000E+00    2.72050000E-01 
S    1    1.000000
   6.66900000E-01    1.00000000E+00 
S    1    1.000000
   2.08900000E-01    1.00000000E+00 
P    1    1.000000
   3.04400000E+00    1.00000000E+00 
P    1    1.000000
   7.58000000E-01    1.00000000E+00 
D    1    1.000000
   1.96500000E+00    1.00000000E+00 
S    1    1.000000
   5.13800000E-02    1.00000000E+00 
P    1    1.000000
   1.99300000E-01    1.00000000E+00 
D    1    1.000000
   4.59200000E-01    1.00000000E+00 
S    1    1.000000
   2.44564000E-02    1.00000000E+00
S    1    1.000000
   1.2282000E-02    1.00000000E+00
S    1    1.000000
   6.1141000E-03    1.00000000E+00
P   1   1.0
    8.130000e-02    1.0
P   1   1.0
    4.065000e-02    1.0
P   1   1.0
    2.032500e-02    1.0
D   1   1.0
    2.34375e-01    1.0
D   1   1.0
    1.17187e-01    1.0
D   1   1.0
    5.85937e-02    1.0
****
$end

Example 6.141 CAP-augmented EOM-EA-CCSD calculation for N2-. aug-cc-pVTZ basis augmented by the 3s3p3d diffuse functions placed in the COM. 2 EA states are computed for CAP strength eta=0.002

$molecule
0 1
N   0.0  0.0 -0.54875676501
N   0.0  0.0  0.54875676501
Gh  0.0  0.0  0.0
$end

$rem
complex_ccman      1             engage complex_ccman          
jobtype            SP            single point
METHOD             EOM-CCSD
BASIS              gen           use general basis
EA_STATES          [0,0,2,0,0,0,0,0]  compute electron attachment energies
cc_memory          5000          ccman2 memory
mem_total          2000
CC_EOM_PROP        true          compute excited state properties
$end

$complex_ccman
CS_HF 1                          Use complex HF
CAP_ETA 200                      Set strength of CAP potential 0.002
CAP_X 2760                       Set length of the box along x dimension
CAP_Y 2760                       Set length of the box along y dimension
CAP_Z 4880                       Set length of the box along z dimension
CAP_TYPE 1                       Use cuboid CAP
$end

$basis
N    0
aug-cc-pvtz
****
Gh   0
S    1    1.000000
   2.88000000E-02    1.00000000E+00
S    1    1.000000
   1.44000000E-02    1.00000000E+00
S    1    1.000000
   0.72000000E-02    1.00000000E+00
P    1    1.000000
   2.45000000E-02    1.00000000E+00
P    1    1.000000
   1.22000000E-02    1.00000000E+00
P    1    1.000000
   0.61000000E-02    1.00000000E+00
D    1    1.000000
   0.755000000E-01    1.00000000E+00
D    1    1.000000
   0.377500000E-01    1.00000000E+00
D    1    1.000000
   0.188750000E-01    1.00000000E+00
****
$end

6.7.11 Non-Hartree-Fock Orbitals in EOM Calculations

In cases of problematic open-shell references, e.g., strongly spin-contaminated doublet, triplet or quartet states, one may choose to use DFT orbitals. This can be achieved by first doing DFT calculation and then reading the orbitals and turning Hartree-Fock off. A more convenient way is just to specify EXCHANGE, e.g., if EXCHANGE=B3LYP, B3LYP orbitals will be computed and used in the CCMAN/CCMAN2 module.

6.7.12 Analytic Gradients and Properties for the CCSD and EOM-XX-CCSD Methods

Analytic gradients are available for the CCSD and all EOM-CCSD methods for both closed- and open-shell references (UHF and RHF only), including frozen core/virtual functionality [362] (see also Section 5.12).

Application limit: same as for the single-point CCSD or EOM-CCSD calculations.

Limitations: Gradients for ROHF and non-HF (e.g., B3LYP) orbitals are not yet available.

For the CCSD and EOM-CCSD wavefunctions, Q-Chem currently can calculate permanent and transition dipole moments, oscillator strengths, $\ensuremath{\langle }R^2\ensuremath{\rangle }$ (as well as XX, YY and ZZ components separately, which is useful for assigning different Rydberg states, e.g., $3p_ x$ vs. $3s$ , etc.), and the $\ensuremath{\langle }S^2\ensuremath{\rangle }$ values. Interface of the CCSD and EOM-CCSD codes with the NBO 5.0 package is also available. Furthermore, excited state analyses can be requested for EOM-CCSD excited states. Similar functionality is available for some EOM-OD and CI models.

The coupled-cluster package in Q-Chem can calculate properties of target EOM states including transition dipoles, two-photon absorption transition moments (and corss sections), and geometry optimizations. The target state of interest is selected by CC_STATE_TO_OPT $rem, which specifies the symmetry and the number of the EOM state.

Analysis of the EOM-CC wavefunctions can also be performed; see Section 10.2.7.

Users must be aware of the point group symmetry of the system being studied and also the symmetry of the excited (target) state of interest. It is possible to turn off the use of symmetry using the CC_SYMMETRY. If set to FALSE the molecule will be treated as having $C_1$ symmetry and all states will be of $A$ symmetry.

6.7.12.1 Transition moments and cross sections for two-photon absorption within EOM-EE-CCSD

Calculation of transition moments and cross sections for two-photon absorption for EOM-EE-CCSD wave functions is available in Q-Chem (CCMAN2 only). Both CCSD-EOM and EOM-EOM transitions can be computed. The formalism is described in Ref. [373]. This feature is available both for canonical and RI/CD implementations. Relevant keywords are CC_EOM_2PA (turns on the calculation), CC_STATE_TO_OPT (used for EOM-EOM transitions); additional customization can be performed using the $2pa section.

The $2pa section is used to specify the range of frequency-pairs satisfying the resonance condition. If $2pa section is absent in the input, the transition moments are computed for 2 degenerate photons with total energy matching the excitation energy of each target EOM state (for CCSD-EOM) or each EOM-EOM energy difference (for EOM-EOM transitions): $2 h\nu =E_{ex}$

$2pa                   Non-degenerate resonant 2PA
N_2PA_POINTS 6         Number of frequency pairs
OMEGA_1 500000 10000   Scans 500 cm$^{-1}$ to 550 cm$^{-1}$ 
                       in steps of 10 cm$^{-1}$
$end

N_2PA_POINTS is the number of frequency pairs across the spectrum. The first value associated with OMEGA_1 is the frequency $\times 1000$ in cm $^{-1}$ at the start of the spectrum and the second value is the step size $\times 1000$ in cm $^{-1}$ . The frequency of the second photon at each step is determined within the code as the excitation energy minus OMEGA_1.

6.7.12.2 Calculations of spin-orbit couplings using EOM-CC wavefunctions

Calculations of spin-orbit couplings (SOCs) for EOM-CC wavefunctions is available[374] in CCMAN2. We employ a perturbative approach in which SOCs are computed as matrix elements of the respective part of the Breit-Pauli Hamiltonian using zero-order non-relativistic wave functions. Both the full two-electron treatment and the mean-field approximation (a partial account of the two-electron contributions) are available for the EOM-EE/SF/IP/EA wavefunctions, as well as between the CCSD reference and EOM-EE/SF. To enable SOC calculation, transition properties between EOM states must be enabled via CC_TRANS_PROP, and SOC requested using CC_CALC_SOC. By default, one-electron and mean-field two-electron couplings will be computed. Full two-electron coupling calculation is activated by setting CC_EOM_PROP_TE.

As with other EOM transition properties, the initial EOM state is set by CC_STATE_TO_OPT, and couplings are computed between that state and all other EOM states. In the absense of CC_STATE_TO_OPT, SOCs are computed between the reference state and all EOM-EE or EOM-SF states.

6.7.13 EOM-CC Optimization and Properties Job Control

CC_STATE_TO_OPT

Specifies which state to optimize (or from which state compute EOM-EOM inter-state properties).

TYPE:

INTEGER ARRAY

DEFAULT:

None

OPTIONS:

[ $i$ , $j$ ]

optimize the $j$ th state of the $i$ th irrep.

RECOMMENDATION:

None

Note: The state number should be smaller or equal to the number of excited states calculated in the corresponding irrep.

Note: If analytic gradients are not available, the finite difference calculations will be performed and the symmetry will be turned off. In this case, CC_STATE_TO_OPT should be specified assuming C $_1$ symmetry, i.e., as [1,N] where N is the number of state to optimize (the states are numbered from 1).

CC_EOM_PROP

Whether or not the non-relaxed (expectation value) one-particle EOM-CCSD target state properties will be calculated. The properties currently include permanent dipole moment, the second moments $\ensuremath{\langle }X^2\ensuremath{\rangle }$ , $\ensuremath{\langle }Y^2\ensuremath{\rangle }$ , and $\ensuremath{\langle }Z^2\ensuremath{\rangle }$ of electron density, and the total $\ensuremath{\langle }R^2\ensuremath{\rangle }= \ensuremath{\langle }X^2\ensuremath{\rangle }+\ensuremath{\langle }Y^2\ensuremath{\rangle }+\ensuremath{\langle }Z^2\ensuremath{\rangle }$ (in atomic units). Incompatible with JOBTYPE=FORCE, OPT, FREQ.

TYPE:

LOGICAL

DEFAULT:

FALSE (no one-particle properties will be calculated)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

Additional equations (EOM-CCSD equations for the left eigenvectors) need to be solved for properties, approximately doubling the cost of calculation for each irrep. Sometimes the equations for left and right eigenvectors converge to different sets of target states. In this case, the simultaneous iterations of left and right vectors will diverge, and the properties for several or all the target states may be incorrect! The problem can be solved by varying the number of requested states, specified with XX_STATES, or the number of guess vectors (EOM_NGUESS_SINGLES). The cost of the one-particle properties calculation itself is low. The one-particle density of an EOM-CCSD target state can be analyzed with NBO package by specifying the state with CC_STATE_TO_OPT and requesting NBO=TRUE and CC_EOM_PROP=TRUE.

CC_TRANS_PROP

Whether or not the transition dipole moment (in atomic units) and oscillator strength for the EOM-CCSD target states will be calculated. By default, the transition dipole moment is calculated between the CCSD reference and the EOM-CCSD target states. In order to calculate transition dipole moment between a set of EOM-CCSD states and another EOM-CCSD state, the CC_STATE_TO_OPT must be specified for this state.

TYPE:

LOGICAL

DEFAULT:

FALSE (no transition dipole and oscillator strength will be calculated)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

Additional equations (for the left EOM-CCSD eigenvectors plus lambda CCSD equations in case if transition properties between the CCSD reference and EOM-CCSD target states are requested) need to be solved for transition properties, approximately doubling the computational cost. The cost of the transition properties calculation itself is low.

CC_EOM_2PA

Whether or not the transition moments and cross sections for two-photon absorption will be calculated. By default, the transition moments are calculated between the CCSD reference and the EOM-CCSD target states. In order to calculate transition moments between a set of EOM-CCSD states and another EOM-CCSD state, the CC_STATE_TO_OPT must be specified for this state.

TYPE:

INTEGER

DEFAULT:

0 (do not compute 2PA transition moments)

OPTIONS:

1

Compute 2PA using the fastest algorithm (use $\tilde{\sigma }$ -intermediates for canonical

and $\sigma$ -intermediates for RI/CD response calculations).

2

Use $\sigma$ -intermediates for 2PA response equation calculations.

3

Use $\tilde{\sigma }$ -intermediates for 2PA response equation calculations.

RECOMMENDATION:

Additional response equations (6 for each target state) will be solved, which increases the cost of calculations. The cost of 2PA moments is about 10 times that of energy calculation. Use default algorithm. Setting CC_EOM_2PA $>$ 0 turns on CC_TRANS_PROP.

CC_CALC_SOC

Whether or not the spin-orbit couplings between CC/EOM electronic states will be calculated. By default, the couplings are calculated between the CCSD reference and the EOM-CCSD target states. In order to calculate couplings between EOM states, CC_STATE_TO_OPT must specify the initial EOM state.

TYPE:

LOGICAL

DEFAULT:

FALSE (no spin-orbit couplings will be calculated)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

One-electron and mean-field two-electron SOCs will be computed by default. To enable full two-electron SOCs, two-particle EOM properties must be turned on (see CC_EOM_PROP_TE).

EOM_REF_PROP_TE

Request for calculation of non-relaxed two-particle EOM-CC properties. The two-particle properties currently include $\ensuremath{\langle }S^2\ensuremath{\rangle }$ . The one-particle properties also will be calculated, since the additional cost of the one-particle properties calculation is inferior compared to the cost of $\ensuremath{\langle }S^2\ensuremath{\rangle }$ . The variable CC_EOM_PROP must be also set to TRUE. Alternatively, CC_CALC_SSQ can be used to request $\ensuremath{\langle }S^2\ensuremath{\rangle }$ calculation.

TYPE:

LOGICAL

DEFAULT:

FALSE

(no two-particle properties will be calculated)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

The two-particle properties are computationally expensive since they require calculation and use of the two-particle density matrix (the cost is approximately the same as the cost of an analytic gradient calculation). Do not request the two-particle properties unless you really need them.

CC_FULLRESPONSE

Fully relaxed properties (including orbital relaxation terms) will be computed. The variable CC_EOM_PROP must be also set to TRUE.

TYPE:

LOGICAL

DEFAULT:

FALSE

(no orbital response will be calculated)

OPTIONS:

FALSE, TRUE

RECOMMENDATION:

Not available for non-UHF/RHF references. Only available for EOM/CI methods for which analytic gradients are available.

CC_SYMMETRY

Controls the use of symmetry in coupled-cluster calculations

TYPE:

LOGICAL

DEFAULT:

TRUE

OPTIONS:

TRUE

Use the point group symmetry of the molecule

FALSE

Do not use point group symmetry (all states will be of $A$ symmetry).

RECOMMENDATION:

It is automatically turned off for any finite difference calculations, e.g. second derivatives.

STATE_ANALYSIS

Activates excited state analyses.

TYPE:

LOGICAL

DEFAULT:

FALSE (no excited state analyses)

OPTIONS:

TRUE, FALSE

RECOMMENDATION:

Set to TRUE if excited state analysis is required, but also if plots of densities or orbitals are needed. For details see section 10.2.7.

6.7.14 Examples

Example 6.142 Geometry optimization for the excited open-shell singlet state, $^1B_2$ , of methylene followed by the calculations of the fully relaxed one-electron properties using EOM-EE-CCSD

$molecule
0 1
C
H  1 rCH
H  1 rCH 2 aHCH

rCH    = 1.083
aHCH   = 145.
$end

$rem
jobtype            OPT
METHOD             EOM-CCSD
BASIS              cc-pVTZ
SCF_GUESS          CORE
SCF_CONVERGENCE    9
EE_SINGLETS        [0,0,0,1]
EOM_NGUESS_SINGLES  2
cc_state_to_opt    [4,1]
EOM_DAVIDSON_CONVERGENCE  9    use tighter convergence for EOM amplitudes
$end

@@@
$molecule
READ
$end

$rem
jobtype            SP
METHOD             EOM-CCSD
BASIS              cc-pVTZ
SCF_GUESS          READ
EE_SINGLETS        [0,0,0,1]
EOM_NGUESS_SINGLES  2
CC_EOM_PROP         1  calculate properties for EOM states
CC_FULLRESPONSE     1  use fully relaxed properties
$end

Example 6.143 Property and transition property calculation on the lowest singlet state of CH $_2$ using EOM-SF-CCSD

$molecule
   0 3
   C
   H  1 rch
   H  1 rch 2 ahch

  rch  = 1.1167 
  ahch = 102.07
$end

$rem
   METHOD             eom-ccsd
   BASIS              cc-pvtz
   SCF_GUESS          core
   SCF_CONVERGENCE    9
   SF_STATES          [2,0,0,3]   Get three 1^B2 and two 1^A1 SF states
   CC_EOM_PROP        1
   CC_TRANS_PROP      1
   CC_STATE_TO_OPT    [4,1] First EOM state in the 4th irrep 
$end

Example 6.144 Geometry optimization with tight convergence for the $^2$ A $_1$ excited state of CH $_2$ Cl, followed by calculation of non-relaxed and fully relaxed permanent dipole moment and $\langle S^2\rangle$ .

$molecule
0 2
H
C 1 CH
CL 2 CCL 1 CCLH
H 2 CH 3 CCLH 1 DIH

CH=1.096247
CCL=2.158212
CCLH=122.0
DIH=180.0
$end

$rem
JOBTYPE            OPT
METHOD             EOM-CCSD
BASIS              6-31G*  Basis Set
SCF_GUESS          SAD
EOM_DAVIDSON_CONVERGENCE    9    EOM amplitude convergence
CC_T_CONV          9    CCSD amplitudes convergence
EE_STATES          [0,0,0,1]
cc_state_to_opt    [4,1]
EOM_NGUESS_SINGLES  2
GEOM_OPT_TOL_GRADIENT 2
GEOM_OPT_TOL_DISPLACEMENT 2
GEOM_OPT_TOL_ENERGY 2
$end

@@@
$molecule
READ
$end

$rem
JOBTYPE            SP
METHOD             EOM-CCSD
BASIS              6-31G*  Basis Set
SCF_GUESS          READ
EE_STATES          [0,0,0,1]
CC_NGUESS_SINGLES  2
CC_EOM_PROP        1   calculate one-electron properties
CC_EOM_PROP_TE     1   and two-electron properties (S^2)
$end

@@@
$molecule
READ
$end

$rem
JOBTYPE            SP
METHOD             EOM-CCSD
BASIS              6-31G*  Basis Set
SCF_GUESS          READ
EE_STATES          [0,0,0,1]
EOM_NGUESS_SINGLES  2
CC_EOM_PROP         1  calculate one-electron properties
CC_EOM_PROP_TE      1  and two-electron properties (S^2)CC_EXSTATES_PROP 1
CC_FULLRESPONSE     1  same as above, but do fully relaxed properties
$end

Example 6.145 CCSD calculation on three $A_2$ and one $B_2$ state of formaldehyde. Transition properties will be calculated between the third $A_2$ state and all other EOM states

$molecule
   0  1
   O
   C  1  1.4
   H  2  1.0  1  120
   H  3  1.0  1  120
$end

$rem
   BASIS             6-31+G
   METHOD            EOM-CCSD
   EE_STATES         [0,3,0,1]
   CC_STATE_TO_OPT   [2,3]
   CC_TRANS_PROP     true
$end

Example 6.146 EOM-IP-CCSD geometry optimization of X $^2B_2$ state of $\rm H_2O^+$ .

$molecule
0 1
          H       0.774767     0.000000     0.458565
          O       0.000000     0.000000    -0.114641
          H      -0.774767     0.000000     0.458565
$end

$rem
jobtype        opt
method         eom-ccsd
basis          6-311G
ip_states      [0,0,0,1]
cc_state_to_opt  [4,1]
$end

Example 6.147 Calculating resonant 2PA with degenerate photons.

$molecule
0 1
O
H  1 0.959
H  1 0.959 2 104.654
$end

$rem
method eom-ccsd
basis aug-cc-pvtz
ee_singlets [1,0,0,0]    1A_1 state
CC_TRANS_PROP 1          Compute transition properties
CC_EOM_2PA    1          Calculate 2PA cross sections using the fastest algorithm
$end

Example 6.148 Non-degenerate, resonant 2PA scan over a range of frequency pairs.

$molecule
0 1
O
H  1 0.959
H  1 0.959 2 104.654
$end

$rem
method eom-ccsd
basis aug-cc-pvtz
ee_singlets [2,0,0,0]    Two A_1 states
CC_TRANS_PROP 1          Calculate transition properties
CC_EOM_2PA    1          Calculate 2PA cross sections using the fastest algorithm
$end

$2pa                Non-degenerate resonant 2PA
n_2pa_points 11     Number of frequency pairs
omega_1 500000 5000 Scans 500 cm$^{-1}$ to 550 cm$^{-1}$ in steps of 5 cm$^{-1}$
$end

Example 6.149 Resonant 2PA with degenerate photons between two excited states.

$molecule
0 1
O
H  1 0.959
H  1 0.959 2 104.654
$end

$rem
method eom-ccsd
basis aug-cc-pvtz
ee_singlets     [2,0,0,0] Two A_1 states
STATE_TO_OPT    [1,1]     "Reference" state for transition properties is 1A_1 state 
CC_TRANS_PROP   1         Compute transition properties
CC_EOM_2PA      1         Calculate 2PA cross sections using the fastest algorithm
$end

Example 6.150 Computation of spin-orbit couplings between closed-shell singlet and $M_ s=1$ triplet state in NH using EOM-SF-CCSD

$molecule
0 3
N
H N 1.0450
$end

$rem
jobtype = sp
method = eom-ccsd
basis = 6-31g

sf_states = [1,2,0,0]
cc_trans_prop = true
cc_calc_soc = true
cc_state_to_opt = [1,1]
$end

6.7.15 EOM(2,3) Methods for Higher-Accuracy and Problematic Situations

In the EOM-CC(2,3) approach [375], the transformed Hamiltonian $\bar H$ is diagonalized in the basis of the reference, singly, doubly, and triply excited determinants, i.e., the excitation operator $R$ is truncated at triple excitations. The excitation operator $T$ , however, is truncated at double excitation level, and its amplitudes are found from the CCSD equations, just like for EOM-CCSD [or EOM-CC(2,2)] method.

The accuracy of the EOM-CC(2,3) method closely follows that of full EOM-CCSDT [which can be also called EOM-CC(3,3)], whereas computational cost of the former model is less.

The inclusion of triple excitations is necessary for achieving chemical accuracy (1 kcal/mol) for ground state properties. It is even more so for excited states. In particular, triple excitations are crucial for doubly excited states [375], excited states of some radicals and SF calculations (diradicals, triradicals, bond-breaking) when a reference open-shell state is heavily spin-contaminated. Accuracy of EOM-CCSD and EOM-CC(2,3) is compared in Table 6.7.15.

System	EOM-CCSD	EOM-CC(2,3)
Singly-excited electronic states	0.1–0.2 eV	0.01 eV
Doubly-excited electronic states	$\geq$ 1 eV	0.1–0.2 eV
Severe spin-contamination of the reference	$\sim$ 0.5 eV	$\leq$ 0.1 eV
Breaking single bond (EOM-SF)	0.1–0.2 eV	0.01 eV
Breaking double bond (EOM-2SF)	$\sim$ 1 eV	0.1–0.2 eV

Table 6.3: Performance of the EOM-CCSD and EOM-CC(2,3) methods

The applicability of the EOM-EE/SF-CC(2,3) models to larger systems can be extended by using their active-space variants, in which triple excitations are restricted to semi-internal ones.

Since the computational scaling of EOM-CC(2,3) method is ${\cal {O}}({N^8})$ , these calculations can be performed only for relatively small systems. Moderate size molecules (10 heavy atoms) can be tackled by either using the active space implementation or tiny basis sets. To achieve high accuracy for these systems, energy additivity schemes can be used. For example, one can extrapolate EOM-CCSDT/large basis set values by combining large basis set EOM-CCSD calculations with small basis set EOM-CCSDT ones.

Running the full EOM-CC(2,3) calculations is straightforward, however, the calculations are expensive with the bottlenecks being storage of the data on a hard drive and the CPU time. Calculations with around 80 basis functions are possible for a molecule consisting of four first row atoms (NO dimer). The number of basis functions can be larger for smaller systems.

Note: In EE calculations, one needs to always solve for at least one low-spin root in the first symmetry irrep in order to obtain the correlated EOM energy of the reference. The triples correction to the total reference energy must be used to evaluate EOM-(2,3) excitation energies.

Note: EOM-CC(2,3) works for EOM-EE, EOM-SF, and EOM-IP/EA. In EOM-IP, “triples” correspond to $3h2p$ excitations, and the computational scaling of EOM-IP-CC(2,3) is less.

6.7.16 Active-Space EOM-CC(2,3): Tricks of the Trade

Active space calculations are less demanding with respect to the size of a hard drive. The main bottlenecks here are the memory usage and the CPU time. Both arise due to the increased number of orbital blocks in the active space calculations. In the current implementation, each block can contain from 0 up to 16 orbitals of the same symmetry irrep, occupancy, and spin-symmetry. For example, for a typical molecule of C $_{\ensuremath{\mathrm{2v}}}$ symmetry, in a small/moderate basis set (e.g., TMM in 6-31G*), the number of blocks for each index is:

occupied: $(\alpha + \beta )\times ( a_1 + a_2 + b_1 + b_2) = 2\times 4 = 8$
virtuals: $(\alpha + \beta )\times ( 2 a_1 + a_2 + b_1 + 2 b_2) = 2\times 6 = 12$
(usually there are more than 16 $a_1$ and $b_2$ virtual orbitals).

In EOM-CCSD, the total number of blocks is $O^2V^2 = 8^2 \times 12^2 = 9216$ . In EOM-CC(2,3) the number of blocks in the EOM part is $O^3V^3 = 8^3 \times 12^3 = 884736$ . In active space EOM-CC(2,3), additional fragmentation of blocks occurs to distinguish between the restricted and active orbitals. For example, if the active space includes occupied and virtual orbitals of all symmetry irreps (this will be a very large active space), the number of occupied and virtual blocks for each index is 16 and 20, respectively, and the total number of blocks increases to $3.3\times 10^7$ . Not all of the blocks contain real information, some blocks are zero because of the spatial or spin-symmetry requirements. For the C $_{\ensuremath{\mathrm{2v}}}$ symmetry group, the number of non-zero blocks is about 10–12 times less than the total number of blocks, i.e., $3\times 10^6$ . This is the number of non-zero blocks in one vector. Davidson diagonalization procedure requires (2*MAX_VECTORS + 2*NROOTS) vectors, where MAX_VECTORS is the maximum number of vectors in the subspace, and NROOTS is the number of the roots to solve for. Taking NROOTS=2 and MAX_VECTORS=20, we obtain 44 vectors with the total number of non-zero blocks being $1.3\times 10^8$ .

In CCMAN implementation, each block is a logical unit of information. Along with real data, which are kept on a hard drive at all the times except of their direct usage, each non-zero block contains an auxiliary information about its size, structure, relative position with respect to other blocks, location on a hard drive, and so on. The auxiliary information about blocks is always kept in memory. Currently, the approximate size of this auxiliary information is about 400 bytes per block. It means, that in order to keep information about one vector ( $3\times 10^6$ blocks), 1.2 GB of memory is required! The information about 44 vectors amounts 53 GB. Moreover, the huge number of blocks significantly slows down the code.

To make the calculations of active space EOM-CC(2,3) feasible, we need to reduce the total number of blocks. One way to do this is to reduce the symmetry of the molecule to lower or C $_1$ symmetry group (of course, this will result in more expensive calculation). For example, lowering the symmetry group from C $_{\ensuremath{\mathrm{2v}}}$ to C $_{\ensuremath{\mathrm{s}}}$ would results in reducing the total number of blocks in active space EOM-CC(2,3) calculations in about $2^6 = 64$ times, and the number of non-zero blocks in about 30 times (the relative portion of non-zero blocks in C $_{\ensuremath{\mathrm{s}}}$ symmetry group is smaller compared to that in C $_{\ensuremath{\mathrm{2v}}}$ ).

Alternatively, one may keep the MAX_VECTORS and NROOTS parameters of Davidson’s diagonalization procedure as small as possible (this mainly concerns the MAX_VECTORS parameter). For example, specifying MAX_VECTORS = 12 instead of 20 would require 30% less memory.

One more trick concerns specifying the active space. In a desperate situation of a severe lack of memory, should the two previous options fail, one can try to modify (increase) the active space in such a way that the fragmentation of active and restricted orbitals would be less. For example, if there is one restricted occupied $b_1$ orbital and one active occupied $B_1$ orbital, adding the restricted $b_1$ to the active space will reduce the number of blocks, by the price of increasing the number of FLOPS. In principle, adding extra orbital to the active space should increase the accuracy of calculations, however, a special care should be taken about the (near) degenerate pairs of orbitals, which should be handled in the same way, i.e., both active or both restricted.

6.7.17 Job Control for EOM-CC(2,3)

EOM-CC(2,3) is invoked by METHOD=EOM-CC(2,3). The following options are available:

EOM_PRECONV_SD

Solves the EOM-CCSD equations, prints energies, then uses EOM-CCSD vectors as initial vectors in EOM-CC(2,3). Very convenient for calculations using energy additivity schemes.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$

Do $n$ SD iterations

RECOMMENDATION:

Turning this option on is recommended

CC_REST_AMPL

Forces the integrals, $T$ , and $R$ amplitudes to be determined in the full space even though the CC_REST_OCC and CC_REST_VIR keywords are used.

TYPE:

INTEGER

DEFAULT:

1

OPTIONS:

0

Do apply restrictions

1

Do not apply restrictions

RECOMMENDATION:

None

CC_REST_TRIPLES

Restricts $R_3$ amplitudes to the active space, i.e., one electron should be removed from the active occupied orbital and one electron should be added to the active virtual orbital.

TYPE:

INTEGER

DEFAULT:

1

OPTIONS:

1

Applies the restrictions

RECOMMENDATION:

None

CC_REST_OCC

Sets the number of restricted occupied orbitals including frozen occupied orbitals.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$

Restrict $n$ occupied orbitals.

RECOMMENDATION:

None

CC_REST_VIR

Sets the number of restricted virtual orbitals including frozen virtual orbitals.

TYPE:

INTEGER

DEFAULT:

0

OPTIONS:

$n$

Restrict $n$ virtual orbitals.

RECOMMENDATION:

None

To select the active space, orbitals can be reordered by specifying the new order in the $reorder_mosection. The section consists of two rows of numbers ( $\alpha$ and $\beta$ sets), starting from $1$ , and ending with $n$ , where $n$ is the number of the last orbital specified.

Example 6.151 Example $reorder_mosection with orbitals 16 and 17 swapped for both $\alpha$ and $\beta$ electrons

$reorder_mo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 17 16
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 17 16
$end

6.7.18 Examples

Example 6.152 EOM-SF(2,3) calculations of methylene.

$molecule
0 3
C 
H 1 CH
H 1 CH 2 HCH

CH  = 1.07
HCH = 111.0
$end

$rem
   method             eom-cc(2,3)
   basis              6-31G
   sf_states          [2,0,0,2]
   n_frozen_core      1
   n_frozen_virtual   1
   eom_preconv_sd      20 Get EOM-CCSD energies first (max_iter=20).
$end

Example 6.153 This is active-space EOM-SF(2,3) calculations for methane with an elongated CC bond. HF MOs should be reordered as specified in the $reorder_mosection such that active space for triples consists of sigma and sigma* orbitals.

$molecule
0 3
C
H 1 CH
H 1 CHX 2 HCH
H 1 CH 2 HCH 3 A120
H 1 CH 2 HCH 4 A120

CH=1.086
HCH=109.4712206
A120=120.
CHX=1.8
$end

$rem
   jobtype            sp
   method             eom-cc(2,3)
   basis              6-31G*
   sf_states          [1,0]
   n_frozen_core      1   
   eom_preconv_sd      20   does eom-ccsd first, max_iter=20
   cc_rest_triples   1   triples are restricted to the active space only
   cc_rest_ampl      0   ccsd and eom singles and doubles are full-space 
   cc_rest_occ        4   specifies active space
   cc_rest_vir        17  specifies active space
   print_orbitals     10  (number of virtuals to print)
$end

$reorder_mo
1 2 5 4 3 
1 2 3 4 5 
$end

Example 6.154 EOM-IP-CC(2,3) calculation of three lowest electronic states of water cation.

$molecule
0 1
          H       0.774767     0.000000     0.458565
          O       0.000000     0.000000    -0.114641
          H      -0.774767     0.000000     0.458565
$end

$rem
jobtype      sp
method       eom-cc(2,3)
basis        6-311G
ip_states    [1,0,1,1]
$end

6.7.19 Non-Iterative Triples Corrections to EOM-CCSD and CCSD

The effect of triple excitations to EOM-CCSD energies can be included via perturbation theory in an economical $N^7$ computational scheme. Using EOM-CCSD wavefunctions as zero-order wavefunctions, the second order triples correction to the $\mu$ th EOM-EE or SF state is:

$\begin{equation} \Delta E^{(2)}_{\mu } = - \frac{1}{36} \sum _{i,j,k} \sum _{a,b,c} \frac{ {\tilde\sigma _{ijk}^{abc}}(\mu ) {\sigma _{ijk}^{abc}}(\mu ) }{D_{ijk}^{abc} - \omega _{\mu }} \end{equation}$

(6.37)

where $i,j$ and $k$ denote occupied orbitals, and $a,b$ and $c$ are virtual orbital indices. $\omega _{\mu }$ is the EOM-CCSD excitation energy of the $\mu$ th state. The quantities $\tilde{\sigma }$ and $\sigma$ are:

	$\displaystyle \tilde{\sigma }_{ijk}^{abc}(\mu )$	$\displaystyle =$	$\displaystyle \langle \Phi _{0}\| ({L_{1}}_{\mu } + {L_{2}}_{\mu }) (H e^{(T_{1}+T_{2})})_{c} \| \Phi _{ijk}^{abc}\rangle$		(6.38)
	$\displaystyle \nonumber \sigma _{ijk}^{abc}(\mu )$	$\displaystyle =$	$\displaystyle \langle \Phi _{ijk}^{abc}\|[H e^{(T_{1}+T_{2})} ({R_{0}}_{\mu } + {R_{1}}_{\mu } + {R_{2}}_{\mu })]_{c} \| \Phi _{0}\rangle$

where, the $L$ and $R$ are left and right eigen-vectors for $\mu$ th state. Two different choices of the denominator, $D_{ijk}^{abc}$ , define the (dT) and (fT) variants of the correction. In (fT), $D_{ijk}^{abc}$ is just Hartree-Fock orbital energy differences. A more accurate (but not fully orbital invariant) (dT) correction employs the complete three body diagonal of $\bar H$ , $\langle \Phi _{ijk}^{abc}| (He^{(T_{1}+T_{2})})_{C}|\Phi _{ijk}^{abc}\rangle$ , $D_{ijk}^{abc}$ as a denominator. For the reference (e.g., a ground-state CCSD wavefunction), the (fT) and (dT) corrections are identical to the CCSD(2) $_{T}$ and CR-CCSD(T) $_{L}$ corrections of Piecuch and co-workers [376].

The EOM-SF-CCSD(dT) and EOM-SF-CCSD(fT) methods yield a systematic improvement over EOM-SF-CCSD bringing the errors below 1 kcal/mol. For theoretical background and detailed benchmarks, see Ref. Manohar:2008.

Similar corrections are available for EOM-IP-CCSD [378], where triples correspond to $3h2p$ excitations.

6.7.20 Job Control for Non-Iterative Triples Corrections

Triples corrections are requested by using METHOD or EOM_CORR:

METHOD

Specifies the calculation method.

TYPE:

STRING

DEFAULT:

No default value

OPTIONS:

EOM-CCSD(DT)

EOM-CCSD(dT), available for EE, SF, and IP

EOM-CCSD(FT)

EOM-CCSD(fT), available for EE, SF, and IP

EOM-CCSD(ST)

EOM-CCSD(sT), available for IP

RECOMMENDATION:

None

EOM_CORR

Specifies the correlation level.

TYPE:

STRING

DEFAULT:

None

No correction will be computed

OPTIONS:

SD(DT)

EOM-CCSD(dT), available for EE, SF, and IP

SD(FT)

EOM-CCSD(fT), available for EE, SF, and IP

SD(ST)

EOM-CCSD(sT), available for IP

RECOMMENDATION:

None

6.7.21 Examples

Example 6.155 EOM-EE-CCSD(fT) calculation of CH $^{+}$

$molecule
1 1
C
H  C CH

CH  = 2.137130
$end

$rem
input_bohr         true 
jobtype            sp
method             eom-ccsd(ft)
basis              general
ee_states          [1,0,1,1]
eom_davidson_max_iter  60  increase number of Davidson iterations
$end

$basis
 H   0
 S   3  1.00
       19.24060000         0.3282800000E-01
       2.899200000         0.2312080000
      0.6534000000         0.8172380000
 S   1  1.00
      0.1776000000          1.000000000
 S   1  1.00
      0.0250000000          1.000000000
 P   1  1.00
       1.00000000          1.00000000
 ****
 C   0
 S   6  1.00
       4232.610000         0.2029000000E-02
       634.8820000         0.1553500000E-01
       146.0970000         0.7541100000E-01
       42.49740000         0.2571210000
       14.18920000         0.5965550000
       1.966600000         0.2425170000
 S   1  1.00
       5.147700000          1.000000000
 S   1  1.00
      0.4962000000          1.000000000
 S   1  1.00
      0.1533000000          1.000000000
 S   1  1.00
      0.0150000000          1.000000000
 P   4  1.00
       18.15570000         0.1853400000E-01
       3.986400000         0.1154420000
       1.142900000         0.3862060000
      0.3594000000         0.6400890000
 P   1  1.00
      0.1146000000          1.000000000
 P   1  1.00
      0.0110000000          1.000000000
 D   1  1.00
      0.750000000          1.00000000
 ****
$end

Example 6.156 EOM-SF-CCSD(dT) calculations of methylene

$molecule
0 3
C 
H 1 CH
H 1 CH 2 HCH

CH  = 1.07
HCH = 111.0
$end

$rem
 method           eom-ccsd(dt)
 basis            6-31G
 sf_states        [2,0,0,2]
 n_frozen_core    1
 n_frozen_virtual 1
$end

Example 6.157 EOM-IP-CCSD(dT) calculations of Mg

$molecule
0 1
          Mg       0.000000     0.000000    0.000000
$end

$rem
jobtype      sp
method       eom-ccsd(dt)
basis        6-31g
ip_states    [1,0,0,0,0,1,1,1]
$end

6.7.22 Potential Energy Surface Crossing Minimization

Potential energy surface crossing optimization procedure finds energy minima on crossing seams. On the seam, the potential surfaces are degenerated in the subspace perpendicular to the plane defined by two vectors: the gradient difference

$\begin{equation} {\rm \bf g} = \frac{\partial }{\partial {\rm \bf q}} (E_1 - E_2) \end{equation}$

(6.39)

and the derivative coupling

$\begin{equation} {\rm \bf h} = \left\langle \Psi _1 \left\vert \frac{\partial {\rm \bf H}}{\partial {\rm \bf q}} \right\vert \Psi _2 \right\rangle \end{equation}$

(6.40)

At this time Q-Chem is unable to locate crossing minima for states which have non-zero derivative coupling. Fortunately, often this is not the case. Minima on the seams of conical intersections of states of different multiplicity can be found as their derivative coupling is zero. Minima on the seams of intersections of states of different point group symmetry can be located as well.

To run a PES crossing minimization, CCSD and EOM-CCSD methods must be employed for the ground and excited state calculations respectively.

6.7.22.1 Job Control Options

XOPT_STATE_1, XOPT_STATE_2

Specify two electronic states the intersection of which will be searched.

TYPE:

[INTEGER, INTEGER, INTEGER]

DEFAULT:

No default value (the option must be specified to run this calculation)

OPTIONS:

[spin, irrep, state]

spin = 0

Addresses states with low spin,

see also EE_SINGLETS.

spin = 1

Addresses states with high spin,

see also EE_TRIPLETS.

irrep

Specifies the irreducible representation to which

the state belongs, for $C_{2v}$ point group symmetry

irrep = 1 for $A_1$ , irrep = 2 for $A_2$ ,

irrep = 3 for $B_1$ , irrep = 4 for $B_2$ .

state

Specifies the state number within the irreducible

representation, state = 1 means the lowest excited

state, state = 2 is the second excited state, etc..

0, 0, -1

Ground state.

RECOMMENDATION:

Only intersections of states with different spin or symmetry can be calculated at this time.

Note: The spin can only be specified when using closed-shell RHF references. In the case of open-shell references all states are treated together, see also EE_STATES.

XOPT_SEAM_ONLY

Orders an intersection seam search only, no minimization is to perform.

TYPE:

LOGICAL

DEFAULT:

FALSE

OPTIONS:

TRUE

Find a point on the intersection seam and stop.

FALSE

Perform a minimization of the intersection seam.

RECOMMENDATION:

In systems with a large number of degrees of freedom it might be useful to locate the seam first setting this option to TRUE and use that geometry as a starting point for the minimization.

6.7.22.2 Examples

Example 6.158 Minimize the intersection of Ã $^{2}$ A $_{1}$ and $\~{B}$ $^{2}$ B $_{1}$ states of the NO $_{2}$ molecule using EOM-IP-CCSD method

$molecule
   -1  1
   N1
   O2  N1  rno
   O3  N1  rno  O2  aono

   rno  = 1.3040
   aono = 106.7
$end

$rem
   JOBTYPE                opt         Optimize the intersection seam
   UNRESTRICTED           true
   METHOD                 eom-ccsd
   BASIS                  6-31g
   IP_STATES              [1,0,1,0]   C2v point group symmetry
   EOM_FAKE_IPEA          1
   XOPT_STATE_1           [0,1,1]     1A1 low spin state
   XOPT_STATE_2           [0,3,1]     1B1 low spin state
   GEOM_OPT_TOL_GRADIENT  30          Tighten gradient tolerance
$END

Example 6.159 Minimize the intersection of Ã $^{1}$ B $_{2}$ and $\~{B}$ $^{1}$ A $_{2}$ states of the N $_{3}^{+}$ ion using EOM-CCSD method

$molecule
   1 1
   N1
   N2 N1 rnn
   N3 N2 rnn N1 annn

   rnn=1.46
   annn=70.0
$end

$rem
   JOBTYPE                opt
   METHOD                 eom-ccsd
   BASIS                  6-31g
   EE_SINGLES             [0,2,0,2]   C2v point group symmetry
   XOPT_STATE_1           [0,4,1]     1B2 low spin state
   XOPT_STATE_2           [0,2,2]     2A2 low spin state
   XOPT_SEAM_ONLY         true        Find the seam only
   GEOM_OPT_TOL_GRADIENT  100
$end

$opt
CONSTRAINT                Set constraints on the N-N bond lengths
   stre  1  2  1.46
   stre  2  3  1.46
ENDCONSTRAINT
$end

@@@

$molecule
   READ
$end

$rem
   JOBTYPE                opt         Optimize the intersection seam
   METHOD                 eom-ccsd
   BASIS                  6-31g

   EE_SINGLETS            [0,2,0,2]

   XOPT_STATE_1           [0,4,1]
   XOPT_STATE_2           [0,2,2]
   GEOM_OPT_TOL_GRADIENT  30
$end

6.7.23 Dyson Orbitals for Ionization from Ground and Excited States within EOM-CCSD Formalism

Dyson orbitals can be used to compute total photodetachment/photoionization cross sections, as well as angular distribution of photoelectrons. A Dyson orbital is the overlap between the N-electron molecular wavefunction and the N-1/N+1 electron wavefunction of the corresponding cation/anion:

	$\displaystyle \phi ^ d(1)$	$\displaystyle =$	$\displaystyle {{1}\over {N-1}}\int \Psi ^ N(1, \ldots , n) \Psi ^{N-1}(2, \ldots , n) d2 \ldots dn$		(6.41)
	$\displaystyle \phi ^ d(1)$	$\displaystyle =$	$\displaystyle {{1}\over {N+1}}\int \Psi ^ N(2, \ldots , n+1) \Psi ^{N+1}(1, \ldots , n+1) d2 \ldots d(n+1)$		(6.42)

For the Hartree-Fock wavefunctions and within Koopmans’ approximation, these are just the canonical HF orbitals. For correlated wavefunctions, Dyson orbitals are linear combinations of the reference molecular orbitals:

$\displaystyle \phi ^ d$	$\displaystyle =$	$\displaystyle \sum _ p \gamma _ p \phi _ p$	(6.43)
$\displaystyle \gamma _ p$	$\displaystyle =$	$\displaystyle \langle \Psi ^ N \| p^+ \| \Psi ^{N-1} \rangle$	(6.44)
$\displaystyle \gamma _ p$	$\displaystyle =$	$\displaystyle \langle \Psi ^ N \| p \| \Psi ^{N+1} \rangle$	(6.45)

The calculation of Dyson orbitals is straightforward within the EOM-IP/EA-CCSD methods, where cation/anion and initial molecule states are defined with respect to the same MO basis. Since the left and right CC vectors are not the same, one can define correspondingly two Dyson orbitals (left-right and right-left):

	$\displaystyle \gamma _ p^ R$	$\displaystyle =$	$\displaystyle \langle \Phi _0 e^{T_1+T_2} L^{EE} \|p^+\| R^{IP} e^{T_1+T_2} \Phi _0\rangle$		(6.46)
	$\displaystyle \gamma _ p^ L$	$\displaystyle =$	$\displaystyle \langle \Phi _0 e^{T_1+T_2} L^{IP} \|p\| R^{EE} e^{T_1+T_2} \Phi _0\rangle$		(6.47)

The norm of these orbitals is proportional to the one-electron character of the transition.

Dyson orbitals also offer qualitative insight visualizing the difference between molecular and ionized/attached states. In ionization/photodetachment processes, these orbitals can be also interpreted as the wavefunction of the leaving electron. For additional details, see Refs. Oana:2007,Oana:2009.

6.7.23.1 Dyson Orbitals Job Control

The calculation of Dyson orbitals is implemented for the ground (reference) and excited states ionization/electron attachment. To obtain the ground state Dyson orbitals one needs to run an EOM-IP/EA-CCSD calculation, request transition properties calculation by setting CC_TRANS_PROP=TRUE and CC_DO_DYSON = TRUE. The Dyson orbitals decomposition in the MO basis is printed in the output, for all transitions between the reference and all IP/EA states. At the end of the file, also the coefficients of the Dyson orbitals in the AO basis are available.

Two implementations of Dyson orbitals are currently available: (i) the original implementation in CCMAN; and (ii) new implementation in CCMAN2. The CCMAN implementation is using a diffuse orbital trick (i.e., EOM_FAKE_IPEA will be automatically set to TRUE in these calculations). Note: this implementation has a bug affecting the values of norms of Dyson orbitals (the shapes are correct); thus, using this code is strongly discouraged. CCMAN2 currently includes only EOM-IP and EOM-EA Dyson orbitals (EOM-EE and SF is not yet available). Also, plotting Dyson orbitals in CCMAN2 is not yet available. The CCMAN2 implementation can utilize user-specified EOM guess (using EOM_USER_GUESS keyword and $eom_user_guess section), which is recommended for highly excited states (such as core-ionized states).

CC_DO_DYSON

Whether the reference-state Dyson orbitals will be calculated for EOM-IP/EA-CCSD calculations.

TYPE:

LOGICAL

DEFAULT:

FALSE (the option must be specified to run this calculation)

OPTIONS:

TRUE/FALSE

RECOMMENDATION:

none

For calculating Dyson orbitals between excited states from the reference configuration and IP/EA states, CC_TRANS_PROP=TRUE and CC_DO_DYSON_EE = TRUE have to be added to the usual EOM-IP/EA-CCSD calculation. The IP_STATES keyword is used to specify the target ionized states. The attached states are specified by EA_STATES. The EA-SF states are specified by EOM_EA_BETA. The excited (or spin-flipped) states are specified by EE_STATES and SF_STATES The Dyson orbital decomposition in MO and AO bases is printed for each EE-IP/EA pair of states in the order: EE1 - IP/EA1, EE1 - IP/EA2, $\ldots$ , EE2 - IP/EA1, EE2 - IP/EA2, $\ldots$ , and so on. This feature is only available in CCMAN.

CC_DO_DYSON_EE

Whether excited-state or spin-flip state Dyson orbitals will be calculated for EOM-IP/EA-CCSD calculations.

TYPE:

LOGICAL

DEFAULT:

FALSE (the option must be specified to run this calculation)

OPTIONS:

TRUE/FALSE

RECOMMENDATION:

none

Dyson orbitals can be also plotted using IANLTY = 200 and the $plots utility (CCMAN only). Only the sizes of the box need to be specified, followed by a line of zeros:

$plots
comment
   10   -2   2
   10   -2   2
   10   -2   2
   0     0   0    0
$plots

All Dyson orbitals on the xyz Cartesian grid will be written in the resulting plot.mo file. For RHF(UHF) reference, the columns order in plot.mo is: $\phi ^{lr}_1\alpha \ (\phi ^{lr}_1\beta ) \ \phi ^{rl}_1\alpha \ (\phi ^{rl}_1\beta ) \ \phi ^{lr}_2\alpha \ (\phi ^{lr}_2\beta ) \ \ldots$

In addition, setting the MAKE_CUBE_FILES keyword to TRUE will create cube files for Dyson orbitals which can be viewed with VMD or other programs (see Section 10.6.4 for details). Other means of visualization (e.g., with MOLDEN_FORMAT=TRUE or GUI=2) are currently not available.

6.7.23.2 Examples

Example 6.160 Plotting grd-ex and ex-grd state Dyson orbitals for ionization of the oxygen molecule. The target states of the cation are $^2$ A $_ g$ and $^2$ B $_{2u}$ . Works for CCMAN only.

$molecule
0 3
O 0.000 0.000 0.000
O 1.222 0.000 0.000
$end

$rem
jobtype        sp
basis          6-31G*
method         eom-ccsd
ip_states      [1,0,0,0,0,0,1,0] Target EOM-IP states
cc_trans_prop  true  request transition OPDMs to be calculated
cc_do_dyson    true  calculate Dyson orbitals
IANLTY         200
$end

$plots
plots excited states densities and trans densities
   10   -2   2
   10   -2   2
   10   -2   2
   0     0   0    0
$plots

Example 6.161 Plotting ex-ex state Dyson orbitals between the 1st $^2A_1$ excited state of the HO radical and the the 1st A $_1$ and A $_2$ excited states of HO $^-$ . Works for CCMAN only.

$molecule
-1 1
H 0.000 0.000 0.000
O 1.000 0.000 0.000
$end

$rem
jobtype            SP
method             eom-ccsd
BASIS              6-31G*
ip_states          [1,0,0,0]      states of HO radical
ee_states          [1,1,0,0]      excited states of HO- 
CC_TRANS_PROP true          calculate transition properties  
CC_DO_DYSON_EE true         calculate Dyson orbitals for ionization from ex. states 
IANLTY 200                
$end

$plots
plot excited states densities and trans densities
   10   -2   2
   10   -2   2
   10   -2   2
   0     0   0    0
$plots

Example 6.162 Dyson orbitals for ionization of CO molecule; A $_1$ and B $_1$ ionized states requested.

$molecule
0 1
O
C O 1.131
$end

$rem
jobtype                    SP            single point
correlation                CCSD
basis                      cc-pVDZ
purecart                   111           5d, will be required for ezDyson
ip_states                  [1,0,1,0]     (A1,A2,B1,B2)
ccman2                     true
cc_do_dyson                true
cc_trans_prop              true          necessary for Dyson orbitals job
print_general_basis        true          will be required for ezDyson
$end

Example 6.163 Dyson orbitals for ionization of H $_2$ O; core (A $_1$ ) state requested — ionization from O(1s).


$molecule
0 1
O
H1 O 0.955
H2 O 0.955 H1 104.5
$end

$rem
jobtype                    SP            single point
correlation                CCSD
basis                      cc-pVTZ
purecart                   111           5d, will be required for ezDyson
ip_states                  [1,0,0,0]     (A1,A2,B1,B2)
eom_user_guess             1             on, further defined in $eom_user_guess
ccman2                     true
cc_do_dyson                true
cc_trans_prop              true          necessary for Dyson orbitals job
print_general_basis        true          will be required for ezDyson
$end

$eom_user_guess
1
$end

Example 6.164 Dyson orbitals for ionization of NO molecule using EOM-EA and a closed-shell cation reference; A $_1$ and B $_2$ states requested.

$molecule
+1 1
N   0.00000  0.00000  0.00000
O   0.00000  0.00000  1.02286
$end

$rem
jobtype                    SP            single point
correlation                CCSD
basis                      aug-cc-pVTZ
purecart                   111           5d, will be required for ezDyson
ea_states                  [1,0,0,1]     (A1,A2,B1,B2)
ccman2                     true
cc_do_dyson                true
cc_trans_prop              true          necessary for Dyson orbitals job
print_general_basis        true          will be required for ezDyson
$end

Example 6.165 Dyson orbitals for ionization of triplet O $_2$ and O $_2^-$ at slightly stretched (relative to the equibrium O $_2$ geometry); B $_{3g}$ states are requested.

$comment
EOM-IP-CCSD/6-311+G* and EOM-EA-CCSD/6-311+G* levels of theory, UHF reference.
Start from O2:
1) detach electron - ionizion of neutral (alpha IP).
2) attach electron, use EOM-EA w.f. as initial state - ionization of anion (beta EA).
$end 

$molecule
0 3
O   0.00000  0.00000  0.00000
O   0.00000  0.00000  1.30000
$end

$rem
jobtype                    SP            single point
correlation                CCSD
basis                      6-311(3+)G*
purecart                   2222          6d, will be required for ezDyson
eom_ip_alpha               [0,0,0,1,0,0,0,0]  (Ag,B1g,B2g,B3g,Au,B1u,B2u,B3u)
eom_ea_beta                [0,0,0,1,0,0,0,0]  (Ag,B1g,B2g,B3g,Au,B1u,B2u,B3u)
ccman2                     true
cc_do_dyson                true
cc_trans_prop              true          necessary for Dyson orbitals job
print_general_basis        true          will be required for ezDyson
$end

6.7.24 Interpretation of EOM/CI Wavefunction and Orbital Numbering

Analysis of the leading wavefunction amplitudes is always necessary for determining the character of the state (e.g., HOMO-LUMO excitation, open-shell diradical, etc.). The CCMAN module print out leading EOM/CI amplitudes using its internal orbital numbering scheme, which is printed in the beginning. The typical CCMAN EOM-CCSD output looks like:

Root 1 Conv-d yes Tot Ene= -113.722767530 hartree (Ex Ene 7.9548 eV), 
U1^2=0.858795, U2^2=0.141205 ||Res||=4.4E-07
Right U1:
       Value                   i            ->    a
       0.5358                  7( B2  ) B   ->   17( B2  ) B
       0.5358                  7( B2  ) A   ->   17( B2  ) A
      -0.2278                  7( B2  ) B   ->   18( B2  ) B
      -0.2278                  7( B2  ) A   ->   18( B2  ) A

This means that this state is derived by excitation from occupied orbital #7 (which has $b_2$ symmetry) to virtual orbital #17 (which is also of $b_2$ symmetry). The two leading amplitudes correspond to $\beta \rightarrow \beta$ and $\alpha \rightarrow \alpha$ excitation (the spin part is denoted by $A$ or $B$ ). The orbital numbering for this job is defined by the following map:

The orbitals are ordered and numbered as follows:
Alpha orbitals:
Number  Energy    Type    Symmetry  ANLMAN number  Total number:
  0    -20.613     AOCC      A1      1A1     1
  1    -11.367     AOCC      A1      2A1     2
  2     -1.324     AOCC      A1      3A1     3
  3     -0.944     AOCC      A1      4A1     4
  4     -0.600     AOCC      A1      5A1     5
  5     -0.720     AOCC      B1      1B1     6
  6     -0.473     AOCC      B1      2B1     7
  7     -0.473     AOCC      B2      1B2     8

  0      0.071     AVIRT     A1      6A1     9
  1      0.100     AVIRT     A1      7A1     10
  2      0.290     AVIRT     A1      8A1     11
  3      0.327     AVIRT     A1      9A1     12
  4      0.367     AVIRT     A1     10A1     13
  5      0.454     AVIRT     A1     11A1     14
  6      0.808     AVIRT     A1     12A1     15
  7      1.196     AVIRT     A1     13A1     16
  8      1.295     AVIRT     A1     14A1     17
  9      1.562     AVIRT     A1     15A1     18
 10      2.003     AVIRT     A1     16A1     19
 11      0.100     AVIRT     B1      3B1     20
 12      0.319     AVIRT     B1      4B1     21
 13      0.395     AVIRT     B1      5B1     22
 14      0.881     AVIRT     B1      6B1     23
 15      1.291     AVIRT     B1      7B1     24
 16      1.550     AVIRT     B1      8B1     25
 17      0.040     AVIRT     B2      2B2     26
 18      0.137     AVIRT     B2      3B2     27
 19      0.330     AVIRT     B2      4B2     28
 20      0.853     AVIRT     B2      5B2     29
 21      1.491     AVIRT     B2      6B2     30

The first column is CCMAN’s internal numbering (e.g., 7 and 17 from the example above). This is followed by the orbital energy, orbital type (frozen, restricted, active, occupied, virtual), and orbital symmetry. Note that the orbitals are blocked by symmetries and then ordered by energy within each symmetry block, (i.e., first all occupied $a_1$ , then all $a_2$ , etc.), and numbered starting from 0. The occupied and virtual orbitals are numbered separately, and frozen orbitals are excluded from CCMAN numbering. The two last columns give numbering in terms of the final ANLMAN printout (starting from 1), e.g., our occupied orbital #7 will be numbered as 1 $B_2$ in the final printout. The last column gives the absolute orbital number (all occupied and all virtuals together, starting from 1), which is often used by external visualization routines.

CCMAN2 numbers orbitals by their energy within each irrep keeping the same numbering for occupied and virtual orbitals. This numbering is exactly the same as in the final printout of the SCF wavefunction analysis. Orbital energies are printed next to the respective amplitudes. For example, a typical CCMAN2 EOM-CCSD output will look like that:

 EOMEE-CCSD transition 2/A1
 Total energy = -75.87450159 a.u.  Excitation energy = 11.2971 eV.
 R1^2 = 0.9396  R2^2 = 0.0604  Res^2 = 9.51e-08

 Amplitude    Orbitals with energies
  0.6486       1 (B2) A                  ->    2 (B2) A                   
              -0.5101                          0.1729                     
  0.6486       1 (B2) B                  ->    2 (B2) B                   
              -0.5101                          0.1729                     
 -0.1268       3 (A1) A                  ->    4 (A1) A                   
              -0.5863                          0.0404                     
 -0.1268       3 (A1) B                  ->    4 (A1) B                   
              -0.5863                          0.0404

which means that for this state, the leading EOM amplitude corresponds to the transition from the first b $_2$ orbital (orbital energy $-0.5101$ ) to the second b $_2$ orbital (orbital energy 0.1729).

0	do not pre-converge
N	perform N Davidson iterations pre-converging singles.

1	Compute 2PA using the fastest algorithm (use $\tilde{\sigma }$ -intermediates for canonical
	and $\sigma$ -intermediates for RI/CD response calculations).
2	Use $\sigma$ -intermediates for 2PA response equation calculations.
3	Use $\tilde{\sigma }$ -intermediates for 2PA response equation calculations.

TRUE	Use the point group symmetry of the molecule
FALSE	Do not use point group symmetry (all states will be of $A$ symmetry).

EOM-CCSD(DT)	EOM-CCSD(dT), available for EE, SF, and IP
EOM-CCSD(FT)	EOM-CCSD(fT), available for EE, SF, and IP
EOM-CCSD(ST)	EOM-CCSD(sT), available for IP

[spin, irrep, state]
spin = 0	Addresses states with low spin,
	see also EE_SINGLETS.
spin = 1	Addresses states with high spin,
	see also EE_TRIPLETS.
irrep	Specifies the irreducible representation to which
	the state belongs, for $C_{2v}$ point group symmetry
	irrep = 1 for $A_1$ , irrep = 2 for $A_2$ ,
	irrep = 3 for $B_1$ , irrep = 4 for $B_2$ .
state	Specifies the state number within the irreducible
	representation, state = 1 means the lowest excited
	state, state = 2 is the second excited state, etc..
0, 0, -1	Ground state.

TRUE	Find a point on the intersection seam and stop.
FALSE	Perform a minimization of the intersection seam.