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

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. 682. This feature is available both for canonical and RI/CD implementations. Relevant keywords are CC_EOM_2PA (turns on the calculation, controls NTO calculation), CC_STATE_TO_OPT (used for EOM-EOM transitions); additional customization can be performed using the $2pa section. The quantity printed out is the microscopic cross-section ${\delta }^{TPA}$ (also known as rotationally averaged 2PA strength), as defined in Eq. (30) of Ref. 682.

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): $2h\nu =E_{ex}$ For resonant or near-resonant calculations (i.e., when one of the photons is in resonance with an excited state), damped response theory should be used; such calculations are deployed by adding DAMPED_EPSILON in the $2pa section (the value of $ϵ$ is specified in hartrees, 0.001 is usually sufficient).

$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.

To gain insight into computed cross sections for 2PA, one can perform NTO analysis of the response one-particle density matrices ⁶⁸⁴. To activate NTO analysis of the 2PA response one-particle transition density matrices, set STATE_ANALYSIS = TRUE, MOLDEN_FORMAT = TRUE (to export the orbitals as MolDen files), NTO_PAIRS (specifies the number of orbitals to print). The NTO analysis will be performed for the full 2PA response one-particle transition density matrices as well as the normalized $\omega$DMs (see Ref. 684 for more details).

Capabilities for computing other non-linear properties include RIXS (see section 7.10.6.1) and static and dynamic polarizabilities (section 7.10.18.4).

Calculations of spin–orbit couplings (SOCs) for EOM-CC wave functions is available in CCMAN2^{242, 765}. 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. There are two versions of the code that can be engaged by using CALC_SOC. The old code computes SOC as a matrix element between a pair of states that need to be explicitly computed; thus, it is able to produce only one matrix element per transition. Although this approach can be used to extract the whole matrix through rotations of the molecule⁷⁶⁹, the relative phases of these matrix elements are lost. The new version of the code solves this issue and also improves performance⁷⁶⁵. The new code actively uses Wigner–Eckart’s theorem to compute the entire spin–orbit matrix. We partition excitation operators by their spin properties. Since the spin–orbit operator is a triplet operator, we can consider the following triplet excitation operators

	$T_{p,q}^{1,1}=-a_{p\alpha }^{†}{a}_{q\beta }$		(7.52)
	$T_{p,q}^{1,0}={\displaystyle \frac{1}{\sqrt{2}}}\left(a_{p\alpha }^{†}{a}_{q\alpha }-a_{p\beta }^{†}{a}_{q\beta }\right)$		(7.53)
	$T_{p,q}^{1,-1}=a_{p\beta }^{†}{a}_{q\alpha }$		(7.54)

If we apply Wigner–Eckart’s theorem, the following matrix elements are decomposed into the Clebsch–Gordan coefficients and reduced matrix elements:

$$⟨{\mathrm{\Psi }}_{I}SM|T_{p,q}^{1,k}|{\mathrm{\Psi }}_J{S}^{\prime }{M}^{\prime }⟩=⟨S^{\prime }{M}^{\prime }|1k|SM⟩⟨{\mathrm{\Psi }}_{I}S||T_{p,q}^{1,\cdot }||{\mathrm{\Psi }}_J{S}^{\prime }⟩,$$

(7.55)

where $S$ and $S^{\prime }$ denote spins, and $M$ and $M^{\prime }$ denote spin projections of ${\mathrm{\Psi }}_I$ and ${\mathrm{\Psi }}_J$, respectively. The reduced matrix elements form a new, spinless triplet transition density matrix, which we will denote as $u^{1,\cdot }$. The spin–orbit reduced matrix elements are obtained as a trace with $u^{1,\cdot }$:

	$⟨{\mathrm{\Psi }}_{I}S||H_{L_{-}}||{\mathrm{\Psi }}_J{S}^{\prime }⟩=-{\displaystyle \frac{1}{2}}{\displaystyle \sum _{p,q}}{h}_{L_{-};p,q}^{SO}{u}_{p,q}^{1,\cdot }$		(7.56)
	$⟨{\mathrm{\Psi }}_{I}S||H_{L_0}||{\mathrm{\Psi }}_J{S}^{\prime }⟩={\displaystyle \frac{\sqrt{2}}{2}}{\displaystyle \sum _{p,q}}{h}_{L_z;p,q}^{SO}{u}_{p,q}^{1,\cdot }$		(7.57)
	$⟨{\mathrm{\Psi }}_{I}S||H_{L_{+}}||{\mathrm{\Psi }}_J{S}^{\prime }⟩=+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{p,q}}{h}_{L_{+};p,q}^{SO}{u}_{p,q}^{1,\cdot },$		(7.58)

Applying Wigner–Eckart’s theorem again, all possible matrix elements are generated from reduced matrix elements. This algorithm is available for the EOM-EE/SF/IP/EA wave functions for MP2 and CCSD references, as well as between the CCSD reference and EOM-EE/SF.

   A(Ket)->B(Bra) transition SO matrices
 Analyzing Sz ans S^2 of the pair of states...
 Ket state:  Computed S^2 = 2.019877 will be treated as 2.000000 Sz = 0.000000
 Bra state:  Computed S^2 = 0.019216 will be treated as 0.000000 Sz = 0.000000
 Clebsch-Gordan coefficient: <1.000,0.000;1.000,0.000|0.000,0.000> = 0.577
 One-electron SO (cm-1)
 Reduced matrix elements:
 <S|| Hso(L-) ||S’> = (-2.660965,-22.180075)
 <S|| Hso(L0) ||S’> = (0.000000,-20.019517)
 <S|| Hso(L+) ||S’> = (-2.660965,22.180075)
 SOCC = 21.593626
 Actual matrix elements:
       |Sz=-1.00>           |Sz=0.00>           |Sz=1.00>
 <Sz=0.00|(1.536309,12.805673)(0.000000,-11.558274)(1.536309,-12.805673)
 Mean-field SO (cm-1)
 Reduced matrix elements:
 <S|| Hso(L-) ||S’> = (-2.069801,-9.910572)
 <S|| Hso(L0) ||S’> = (0.000000,-10.692261)
 <S|| Hso(L+) ||S’> = (-2.005570,9.957712)
 Singlet part of <S|| Hso(L0) ||S’> = (-0.000000,0.037101) (excluded from all matrix elements)
 L-/L+ Averaged reduced matrix elements:
 <S|| Hso(L-) ||S’> = (-2.037685,-9.934142)
 <S|| Hso(L+) ||S’> = (-2.037685,9.934142)
 SOCC = 10.328006
 Actual matrix elements:
       |Sz=-1.00>           |Sz=0.00>           |Sz=1.00>
 <Sz=0.00|(1.176458,5.735480)(0.000000,-6.173180)(1.176458,-5.735480)

The algorithm is activated by setting CALC_SOC to 1 or 2. The execution relies on a non-zero Clebsch–Gordan coefficient (otherwise, it is not possible to extract the necessary information from the transition density matrix). If the Clebsch–Gordan coefficient is zero, a warning will be printed, and the transition will be skipped. Usually it is possible to select the states, leading to non-zero Clebsch–Gordan coefficients, by selecting SF or EE, alpha or beta electrons in case of EA or IP. Determination of state multiplicities is based on $S^2$. In case of a strong spin contamination, a warning is printed, and the closest multiplicity is assumed.

In the case of open-shell references, the resulting reduced matrix elements may violate symmetries implied by the spin–orbit operator. Setting CALC_SOC to 1 will produce these reduced matrix elements, but the actual spin–orbit matrix will be produced from the averaged reduced matrix elements, in which the desired symmetry is restored. Setting CALC_SOC to 2 will form spin–orbit matrix without averaging of reduced matrix elements. Although in practice the difference between these two options is small, the averaging may simplify the analysis of splittings.

Note:  The ECP version is not yet implemented. Attempts of using ECP will, most likely, result in unphysically low SOC values. A warning will be printed if an ECP is used.

To activate the NTO analysis of the spinless one-particle transition density matrices ⁷⁶⁷ and to compute spin–orbit integrals over NTOs, set STATE_ANALYSIS = TRUE , MOLDEN_FORMAT = TRUE, in addition to CALC_SOC = 1 or 2.

The legacy code has the full two-electron treatment and the mean-field approximation. To enable the SOC calculation, transition properties between EOM states must be enabled via CC_TRANS_PROP, and SOC requested setting CALC_SOC to 3. 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 absence of CC_STATE_TO_OPT, SOCs are computed between the reference state and all EOM-EE or EOM-SF states.

Note:  In a spin-restricted case, such as EOM-EE calculations using closed-shell reference state, SOCs between the singlet and triplet EOM manifolds cannot be computed (only SOCs between the reference state and EOM triplets can be calculated). To compute SOCs between EOM-EE singlets and EOM-EE triplets, run the same job with UNRESTRICTED = TRUE, such that triplets and singlets appear in the same manifold.

Calculations of non-adiabatic (derivative) couplings (NACs) for EOM-CC wave functions is available in CCMAN2. We employ Szalay’s approach in which couplings are computed by a modified analytic gradient code, via “summed states”:⁹⁴⁴

$$h_{IJ}^x\equiv ⟨{\mathrm{\Psi }}_I|H^x|{\mathrm{\Psi }}_J⟩=\frac{1}{2}\left(G_{(I+J)}-G_I-G_J\right),$$

(7.59)

where, $G_I$, $G_J$, and $G_{IJ}$ are analytic gradients for states $I$, $J$, and a fictitious summed state $|{\mathrm{\Psi }}_{I+J}⟩\equiv |{\mathrm{\Psi }}_I⟩+|{\mathrm{\Psi }}_J⟩$. Currently, NACs for EE/IP/EA are available.²⁵⁰

Note:  Note that the individual components of the NAC vector depend on the molecular orientation.

Calculation of the and dynamical dipole polarizabilities for CCSD and EOM-EE/SF wave functions is available in CCMAN2. Polarizabilities are calculated as second derivatives of the CCSD or EOM-CCSD energy⁶⁸³ as well as using the expectation-value approach, i.e., the sum-over-states expression for the exact wave function⁶⁷⁹. In the analytic derivative approach (or the response-theory formalism), only the response of the cluster amplitudes is taken into the account; orbital relaxation is not included. CC_FULLRESPONSE must be set to 2 to turn off the orbital relaxation but allow amplitude response in EOM calculations, if the polarizabilities are to be computed with this approach (both CC_POL and EOM_POL must equal 1 or 2). CC_POL and EOM_POL = 3 or 4 coupled with CC_FULLRESPONSE = 0 will enable expectation-value polarizabilities. Note that CC_REF_PROP and CC_EOM_PROP/ must be set to TRUE for these calculations.

Dynamical polarizabilities are enabled using the $dyn_pol section, which specifies the frequency range. Without this section, only the static polarizabilities will be computed.

$dyn_pol
   500000 10000 6  !Scans 500 cm$^{-1}$ to 550 cm$^{-1}$
                   !      in steps of 10 cm$^{-1}$
$end

Currently, this feature is available for the canonical implementation only.

Note:  EOM-CCSD polarizabilities are available for the EE and SF wave functions only.

Angular momentum matrix elements can be useful for analyzing wave functions, assigning state characters, or identifying term symmetry labels in linear molecules. As other one-electron properties, these matrix elements are evaluated using one-particle state and transition density matrices. Angular momentum is a gauge-dependent property, meaning that its value depends on the origin of the gauge. The gauge position can be controlled with $gauge_origin input section. The default position of the gauge corresponds to the $(0,0,0)$ of the coordinate system.

$gauge_origin            Position of the gauge origin for
1.00  0.200  0.500       angular momentum calculations in angstroms.
$end