Theory and implementation

approximations is about 0.05 eV, but it starts growing rapidly after the truncation of 15%. It means that only a couple of occupied orbitals can be dropped, and unfortunately the gain in the computation time is rather limited. The approximation is probably more efficient for large systems, for which, however, CC2 calculations are currently not feasible even with the present reduced-cost scheme.

4.2 ADC(2) excitation energies and transition

Table 4.5: Algorithm and working equations for calculating the sigma vector for singlet and triplet excitations

Step Scaling

1 Perform MP2 calculation and store intermediate Y_ai^Q Y_ai^Q =P

jbJ_jb^Q(2t^ab_ij −t^ab_ji) n²_occn²_virtn_aux 2 Calculate intermediates X^Q, X_ij^Q, and X_ia^Q

X^Q = 2P

iaJ_ia^Qr_i^a n_occn_virtn_aux

X_ij^Q =P

aJ_ia^Qr^a_j X_ia^Q =P

br_i^bJ_ba^Q n²_occnvirtnaux+noccn²_virtnaux

3 Add purely r_i^a-dependent contributions to σ_ia and construct F_ia if singlet: σ_ia ←P

QJ_ia^QX^Q n_occn_virtn_aux

F_ia =−P

jQJ_ja^QX_ij^Q+σ_ia n²_occn_virtn_aux+n_occn_virt if triplet: F_ia=P

jQJ_ja^QX_ij^Q n²_occn_virtn_aux

σ_ia ← −P

jQJ_ij^QX_ja^Q + (ε_a−ε_i)r_i^a n²_occn_virtn_aux 4 Compute ˆt^ab_ij, add corresponding contributions to σia

(ia|jb) =P

QJ_ia^QJ_jb^Q n²_occn²_virtn_aux

if singlet: ˆt^ab_ij = [2(ia|jb)−(ja|ib)]/D^ab_ij] n²_occn²_virt if triplet: ˆt^ab_ij = (ja|ib)/D^ab_ij n²_occn²_virt σ_ia ← ¹₂P

jbˆt^ab_ijF_jb n²_occn²_virt

if singlet: σ_ia ← ¹₂P

jb

P

kc[2(ia|jb)−(ja|ib)]r^c_kˆt^bc_jk 2n²_occn²_virt if triplet: σia ← ¹₂P

jb

P

kc(ja|ib)r^c_kˆt^bc_jk n²_occn²_virt 5 Construct the J^Q_ia list, build intermediate ˆr_ij^ab

J^Q_ia=X_ia^Q−P

jJ_ij^Qr_j^a n²_occn_virt+n_occn_virt

(ia|jb) =P

QJ^Q_iaJ_jb^Q (ia|jb) =P

QJ_ia^QJ^Q_jb 2n²_occn²_virtn_aux if singlet:

ˆ

r_ij^ab= (2r^ab_ij −r^ab_ji)/(D_ij^ab+ω_ADC(2)) n²_occn²_virt if triplet:

ˆ

r_ij^ab= [2(ia|jb)−(ja|ib)−(ja|ib)]/(D_ij^ab+ω_ADC(2)) n²_occn²_virt Y^Q_ia =P

jbˆr^ab_ijJ_jb^Q n²_occn²_virtn_aux

6 Add the remaining contributions to σ_ia σ_ia ←P

QbJ_ab^QY^Q_ib n_occn²_virtn_aux

σ_ia ← −P

QjJ_ij^QY^Q_ja n²_occn_virtn_aux

σia ← ¹₂P

j

P

Qb(J_ib^QY_jb^Q+J_jb^QY_ib^Q)r_j^a 2n²_occnvirtnaux+n²_occnvirt

σ_ia ← −¹₂P

b

P

Qj(J_jb^QY_ja^Q+J_ja^QY_jb^Q)r_i^b 2n_occn²_virtn_aux+n_occn²_virt Since we also calculate oscillator strengths in this subsection, the explicit form of the transition density matrix requires some considerations. The expression 3.60 is often

simplified [257,258] by discarding disconnected contributions and by neglecting the higher than fifth-power-scaling second-order terms. It can be shown that, analogously to CC2, the resulting ADC(2) density matrix is consistent with the LR-CC theory and correct up to the first order. Our working equations for the transition density matrix in the spatial orbital basis are given for its various blocks by

ρ_ab=X

ijc

r^ac_ij(2t^bc_ij −t^bc_ji), (4.12)

ρ_ij =−X

abk

r^ab_ik(2t^ab_jk−t^ab_kj), (4.13)

ρ_ai =X

bj

r^b_j(2t^ab_ij −t^ab_ji), (4.14) and

ρ_ia =r_i^a−X

k

r^a_kX

cbj

t^cb_kj(2t^cb_ij −t^cb_ji)−X

c

r^c_i X

bjk

t^cb_kj(2t^ab_kj −t^ab_jk) (4.15)

=r_i^a−X

k

r^a_kD^MP2_ik −X

c

r_i^cD_ac^MP2,

where D_ik^MP2 and D^MP2_ac are the elements, respectively, of the occupied-occupied and the virtual-virtual block of the MP2 one-particle density matrix.

The number of the variables can also be reduced in the case of the ADC(2) method, similar to CC2. For the singlet states no changes needed, but the construction of triplet density matrices to calculate the corresponding VNOs requires different expressions. First we derive the suitable working equations, then we make some comments about the trun-cation of the virtual MO space for both singlet and triplet states.

The expressions for the one-particle MP2 density matrix [Eq. (4.2)] and the CIS density matrix [Eq. (4.3)] are applicable without any changes for triplet states. However, the spin adaptation results in different expressions for matrixD^(D), whose virtual-virtual block is given as

D_ab^(D)=X

ijc

(c^ca_ijc^cb_ij +c^ca_ijc^cb_ij −c^ca_ijc^bc_ij), (4.16) where the c^ab_ij coefficients are still defined by Eq. (4.7), but they are built using the triplet CIS coefficients and excitation energy. Furthermore,c^ab_ij denotes the triplet-coupled CIS(D) doubles coefficient, which is evaluated as

c^ab_ij = P

c[(ac|bj)c^c_i −(ai|bc)c^c_j] +P

k[(kj|ai)c^b_k−(ki|bj)c^a_k]

D_ij^ab+ω_CIS . (4.17)

The algorithm and working equations to compute the triplet density matrix are given in detail in Table 4.6. The rate-determining step of the density matrix calculation is the

Table 4.6: Working equations for the evaluation of the triplet CIS(D) density matrix

Step Scaling

1 Add CIS contribution to the density matrix D^CIS(D)_ab ←P

ic^a_ic^b_i noccn²_virt 2 Compute intermediate Y_ai^Q

Y_ai^Q ←P

cJ_ac^Qc^c_i n_occn²_virtn_aux Y_ai^Q ← −P

jJ_ij^Qc^a_j n²_occnvirtnaux

3 Compute intermediate V_ij^ab V_ij^ab =P

QY_ai^QJ_bj^Q n²_occn²_virtn_aux 4 Calculate c^ab_ij coefficients

c^ab_ij = (V_ij^ab+V_ji^ba)/(D^ab_ij +ω_CIS) n²_occn²_virt 5 Compute intermediate X_ij^ab

X_ij^ab =c^ab_ij −c^ba_ij n²_occn²_virt 6 Add the first contribution to the density matrix

D^CIS(D)_ab ←P

ijcc^ca_ijX_ij^cb n²_occn³_virt 7 Calculate c^ab_ij coefficients

c^ab_ij = (V_ij^ab−V_ji^ba)/(D_ij^ab+ω_CIS) n²_occn²_virt 8 Add the second contribution to the density matrix

D^CIS(D)_ab ←P

ijcc^ca_ijc^cb_ij n²_occn³_virt

assembly of the four-index integral list of intermediate V(step 3), which is a fifth-power-scaling operation with costs proportional to n²_occn²_virtn_aux. The algorithm contains two additional fifth-power-scaling operations (steps 6 and 8) with computational expenses proportional to n²_occn³_virt, while there is only one such term in the case of MP2 or singlet CIS(D). The exchange term in Eq. (4.16) can also be neglected for the triplet density similar to the singlet case, whose benefits and effects have already been discussed in detail in Subsections 4.1.1 and 4.1.2.

As expected, the selected VNOs are only ideal for ADC(2) calculations if the CIS(D) wave function is a good approximation to that for ADC(2). The prerequisite for this is that the CIS method provides a qualitatively good description of the excited state. If the coefficients of double excitations are not stored, a good measure for the quality of the CIS and CIS(D) approximations is the overlap of the CIS and the single excitation part of the final ADC(2) wave functions. Our numerical experience shows that, if this overlap is relatively small, the ADC(2) calculations performed in the reduced VNO space may yield inadequate excitation energies and transition moments. In these cases, we found

that the dominant excitations of the ADC(2) wave function obtained in the reduced space correspond to that of the CIS wave function, and particular excitations that have con-siderable weight in the canonical ADC(2) are missing. This problem cannot be resolved by tightening the ε_VNO threshold as the occupancy of the VNOs which do not noticeably contribute to the CIS wave function are rather low. To overcome this problem we found the following approach to be useful.

In practice, for a particular excited state, the orbital energies of the canonical virtual orbitals contributing to the wave function are relatively close to each other. Thus, we can suppose that all the important orbitals will be included in the reduced space if it is augmented with the canonical virtual orbitals which are close in energy to the virtuals significantly contributing to the CIS wave function. For the selection of such orbitals two further thresholds were introduced. First, those canonical virtual orbitals are chosen for which at least one CIS coefficient in the normalized CIS wave function of an excitation involving them is greater than thresholdε_CIS. Second, all those canonical virtual orbitals are selected whose orbital energy is in a given interval. The lower limit of the interval is the orbital energy of the lowest MO selected in the first step minus threshold εE, while its upper limit is the orbital energy of the highest chosen MO plusε_E. The orbitals of the reduced VNO space are projected out of the selected canonical virtual orbitals, and the remaining orbitals are orthonormalized employing canonical orthogonalization dropping the eigenvectors of the overlap matrix with eigenvalues lower than 10⁻⁷. The resulting orthonormal MOs are added to the reduced VNO space, and the orbitals of the extended space are canonicalized. In practice it means that the Fock-matrix is transformed to the extended MO space, and the matrix is diagonalized. The resulting pseudo-canonical orbitals and the corresponding orbital energies will be used in the subsequent calculations.

For the sake of simplicity, from now on, this pseudo-canonical basis will be referred to as the VNO basis.

Significant gains can also be expected utilizing the NAF approximation for the ADC(2) method. In the following we will employ two different sets of NAFs. First, as in our previous CC2 study, we construct NAFs to compress the representation of the integrals in the NO basis. Here we build W of Eq. (3.66) using ˜J, the integral list transformed to the VNO basis. These NAFs will be referred to as restricted NO space NAFs or RS-NAFs for short. The RS-NAFs are suitable to reduce, for instance, the cost of certain assembly and σ-vector construction steps of the ADC(2) calculation (see, e.g., Table 4.5). Another option is to utilize the NAF approximation right at the beginning, as soon as the three-center integrals are transformed to the canonical HF MO basis. These NAFs will be called the complete MO space NAFs or shortly, CS-NAFs, and the corresponding J integrals will be denoted hereafter by ˆJ. We note here that, naturally, the CS-NAF basis is usually significantly larger than the RS-NAF basis but

still represents a compressed, system-specific auxiliary basis with which the expenses of the CIS calculation and density matrix evaluation steps can be cut roughly by half.

We briefly collect all the steps required for our reduced-cost ADC(2) algorithm, which is basically the generalization of our VNO and (RS-)NAF approximation-based CC2 scheme apart from the step that the VNOs selected using the ε_VNO threshold are supplemented with the close-lying canonical orbitals and CS-NAFs are introduced as well.

Our general algorithm is as follows.

1. Solve HF equations

2. Construct the canonical integral list J

3. Calculate matrix W using three-center integral list J [Eq. (3.66)], and transform the auxiliary index of J to the CS-NAF basis to obtain ˆJ

4. Solve CIS equations for all the excited states simultaneously using ˆJ 5. Loop over excited states

(a) Calculate the state-averaged one-particle density matrix D using ˆJ, and construct the truncated virtual space. Transform the MO indices of ˆJ to the VNO basis to obtain ˜J.

(b) Calculate matrix W using three-index integrals ˜J [Eq. (3.66)], and trans-form the auxiliary index of ˜J to the RS-NAF basis (J). Write integral list J to disk.

End loop

6. Loop over excited states

(a) Retrieve integral list J from disk

(b) Perform MP2 calculation and solve ADC(2) equations (Table4.5) End loop

We would like to point out that the CS-NAF approximation affects the accuracy of all the quantities computed subsequently, such as density matrices, VNOs, RS-NAFs, etc. These negligible errors are extensively investigated in the original paper [239].

In document 2019 D´avidMester EfficientExcited-StateQuantumChemicalApproachesforLargeMolecules (Pldal 42-48)