7.14. The Single Reference Correlation Module¶
ORCA features a variety of single-reference correlation methods for
single point energies (restricted to a RHF or RKS determinant in the
closed-shell case and a UHF or UKS determinant in the open-shell case;
quasi-restricted orbitals (QROs)[611] are also supported in the
open-shell case). They are all fairly expensive but maybe be used in
order to obtain accurate results in the case that the reference
determinant is a good starting point for the expansion of the many-body
wavefunction. The module is called orca_mdci for “matrix driven
configuration interaction”. This is a rather technical term to emphasize
that if one wants to implement these methods (CCSD, QCISD etc.)
efficiently, one needs to write them in terms of matrix operations which
pretty much every computer can drive at peak performance. Let us first
briefly describe the theoretical background of the methods that we have
implemented in ORCA.
7.14.1. Theory¶
We start from the full CI hierarchy in which the wavefunction is expanded as:
where \(\left| 0 \right\rangle\) is a single-determinant reference and S, D, T, Q, …denote the single, double, triple quadruple and higher excitations relative to this determinant at the spin-orbital level. As usual, labels \(i,j,k,l\) refer to occupied orbitals in \(\left| 0 \right\rangle\), \(a,b,c,d\) to unoccupied MOs and \(p,q,r,s\) to general MOs. The action of the second quantized excitation operators \(a_{i}^{a} =a_{a}^{\dagger} a_{i}\) on \(\left| 0 \right\rangle\) lead to excited determinants \(\left| \Phi_i^a \right\rangle\) that enter \(\left| \Psi \right\rangle\) with coefficients \(C_{a}^{i}\). The variational equations are:
Further equations coupling triples with singles through pentuples etc.
The total energy is the sum of the reference energy \(E_{0} =\left\langle {0\left| H \right|0} \right\rangle\) and the correlation energy
which requires the exact singles- and doubles amplitudes to be known. In order to truncate the series to singles- and doubles one may either neglect the terms containing the higher excitations on the right hand side (leading to CISD) or approximate their effect thereby losing the variational character of the CI method (CCSD, QCISD and CEPA methods). Defining the one- and two-body excitation operators as \(\hat{{C} }_{1} =\sum\nolimits_{ia} {C_{a}^{i} a_{i}^{a} }\), \(\hat{{C} }_{2} =\frac{1}{4}\sum\nolimits_{ijab} {C_{ab}^{ij} a_{ij}^{ab} }\) one can proceed to approximate the triples and quadruples by the disconnected terms:
As is well known, the CCSD equations contain many more disconnected contributions arising from the various powers of the \(\hat{{C} }_{1}\) operator (if one would stick to CC logics one would usually label the cluster amplitudes with \(t_{a}^{i} ,t_{ab}^{ij}\),…and the \(n\)-body cluster operators with \(\hat{{T} }_{n}\); we take a CI point of view here). In order to obtain the CEPA type equations from ((7.114)-(7.118)), it is most transparent to relabel the singles and doubles excitations with a compound label \(P\) for the internal indices (\(i\)) or (ij) and \(x\) for (\(a\)) or (ab). Then, the approximations are as follows:
Here the second line contains the approximation that only the terms in which either Qy or Rz are equal to Px are kept (this destroys the unitary invariance) and the fourth line contains the approximation that only “exclusion principle violating” (EPV) terms of internal labels are considered. The notation \(Qy\cup Px\) means “Qy joint with Px” (containing common orbital indices) and \(\varepsilon_{Q}\) is the pair correlation energy. The EPV terms must be subtracted from the correlation energy since they arise from double excitations that are impossible due to the fact that an excitation out of an occupied or into an empty orbital of the reference determinant has already been performed. Inserting eq. (7.122) into eq. (7.115) \(C_{x}^{P} E_{C}\) cancels and effectively is replaced by the “partial correlation energy” \(\sum\nolimits_{Q\cup P} { \varepsilon_{Q} }\).
The resulting equations thus have the appearance of a diagonally shifted
(“dressed”) CISD equation
\(\left\langle { \Phi_{P}^{x} \left|{ H-E_{0}
+\Delta } \right|0+S+D} \right\rangle=0\). If the second approximation
mentioned above is avoided Malrieu’s (SC) \(^{2}\)-CISD arises.
[38, 63, 646, 719, 911]
Otherwise, one obtains CEPA/3 with the shift:
CEPA/2 is obtained by \(-\Delta_{ab}^{ij} =\varepsilon_{ij}\) and CEPA/1 is the average of the CEPA/2 and CEPA/3. As mentioned by Ahlrichs, [11] no consensus appears to exist in the literature for the appropriate shift on the single excitations. If one proceeds straightforwardly in the same way as above, one obtains:
as the appropriate effect of the disconnected triples on the singles. It has been assumed here that only the singles \(\left|\Phi_i^a \right\rangle\) in \(\hat{{C} }_{1}\) contribute to the shift. If \(\left| 0 \right\rangle\) is a HF determinant, the effect of the disconnected triples in the doubles projection vanishes under the same CEPA approximations owing to Brillouin’s theorem. Averaged CEPA models are derived by assuming that all pair correlation energies are equal (except \(\varepsilon_{ii} =0)\). As previously discussed by Gdanitz [291], the averaging of CEPA/1 yields \(\frac{2}{n}E_{\text{C} }\) and CEPA/3 \(E_{\text{C} } \frac{4n-6}{n\left({ n-1} \right)}\) where \(n\) is the number of correlated electrons. These happen to be the shifts used for the averaged coupled-pair functional (ACPF [292]) and averaged quadratic coupled-cluster (AQCC [840]) methods respectively. However, averaging the singles shift of eq. (7.124) gives \(\frac{4}{n}E_{\text{C} }\). The latter is also the leading term in the expansion of the AQCC shift for large \(n\). In view of the instability of ACPF in certain situations, Gdanitz has proposed to use the AQCC shift for the singles and the original ACPF shift for the doubles and called his new method ACPF/2 [292]. Based on what has been argued above, we feel that it would be most consistent with the ACPF approach to simply use \(\frac{4}{n}E_{\text{C} }\) as the appropriate singles shift. We refer to this as NACPF.
It is readily demonstrated that the averaged models may be obtained by a variation of the modified correlation energy functional:
with g\(_{S}\) and g\(_{D}\) being the statistical factors \(\frac{4}{n}\), \(\frac{2}{n}\), \(\frac{4n-6}{n\left({ n-1} \right)}\), as appropriate for the given method. Thus, unlike the CEPA models, the averaged models fulfill a stationarity principle and are unitarily invariant. However, if one thinks about localized internal MOs, it appears evident that the approximation of equal pair energies must be one of rather limited validity and that a more detailed treatment of the electron pairs is warranted. Maintaining a stationarity principle while providing a treatment of the pairs that closely resembles that of the CEPA methods was achieved by Ahlrichs and co-workers in an ingenious way with the development of the CPF method [16]. In this method, the correlation energy functional is written as:
with
The topological matrix for pairs \(P=\)(ij) and \(Q=(kl)\) is chosen as: [440]
with \(n_{i}\) being the number of electrons in orbital \(i\) in the reference determinant. The singles out of orbital \(i\) are formally equated with \(P=(ii)\). At the spin-orbital level, \(n_{i} =1\), for closed shells \(n_{i} =2\). Using the same topological matrix in \(\Delta_{P} =\sum\nolimits_Q { T_{PQ} \varepsilon_{Q} }\) one recovers the CEPA/1 shifts for the doubles in eq. (7.124). It is straightforward to obtain the CPF equivalents of the other CEPA models by adjusting the \(T_{PQ}\) matrix appropriately. In our program, we have done so and we refer below to these methods as CPF/1, CPF/2 and CPF/3 in analogy to the CEPA models (CPF/1 \(\equiv\)CPF). In fact, as discussed by Ahlrichs and co-workers, variation of the CPF-functional leads to equations that very closely resemble the CEPA equation and can be readily implemented along the same lines as a simple modification of a CISD program. Ahlrichs et al. argued that the energies of CEPA/1 and CPF/1 should be very close. We have independently confirmed that in the majority of cases, the total energies predicted by the two methods differ by less than 0.1 mEh.
An alternative to the CPF approach which is also based on variational optimization of an energy functional is the VCEPA method [455]. The equations resulting from application of the variational principle to the VCEPA functional are even closer to the CEPA equations than for CPF so that the resulting energies are practically indistinguishable from the corresponding CEPA values. The VCEPA variants are referred to as VCEPA/1, VCEPA/2, and VCEPA/3 in analogy to CEPA and CPF. A strictly size extensive energy functional (SEOI) which is invariant with respect to unitary transformations within the occupied and virtual orbital subspaces is also available [456] (an open-shell version is not implemented yet).
Again, a somewhat critical point concerns the single excitations. They do not account for a large fraction of the correlation energy. However, large coefficients of the single excitations lead to instability and deterioration of the results. Secondly, linear response properties are highly dependent on the effective energies of the singles and their balanced treatment is therefore important. Since the CEPA and CPF methods amount to shifting down the diagonal energies of the singles and doubles, instabilities are expected if the effective energy of an excitation approaches the reference energy of even falls below it. In the CPF method this would show up as denominators \(N_{P}\) that are too small. The argument that the CPF denominators are too small has led Chong and Langhoff to the proposal of the MCPF method which uses a slightly more elaborate averaging than (\(N_{P}N_{Q})^{1/2}\) [173].[1] However, their modification was solely based on numerical arguments rather than physical or mathematical reasoning. In the light of Eq. (7.124) and the performance of the NACPF, it appears to us that for the singles one should use twice the \(T_{PQ}\) proposed by Ahlrichs and co-workers. The topological matrix \(T_{PQ}\) is modified in the following way for the (very slightly) modified method to which we refer to as NCPF/1:
(note that \(T_{PQ} \ne T_{QP}\) for this choice). Thus, the effect of the singles on the doubles is set to zero based on the analysis of the CEPA approximations and the effect of the singles on the singles is also set to zero. This is a sensible choice since the product of two single excitations is a double excitation which is already included in the SD space and thus none of them can belong to the outer space. It is straightforward to adapt this reasoning about the single excitations to the CEPA versions as well as to NCPF/2 and NCPF/3.
The aforementioned ambiguities arising from the use of single excitations in coupled-pair methods can be avoided by using correlation-adapted orbitals instead of Hartree-Fock orbitals thus eliminating the single excitations. There are two alternatives: (a) Brueckner orbitals and (b) optimized orbitals obtained from the variational optimization of the electronic energy with respect to the orbitals. Both approaches have already been used for the coupled-cluster doubles (CCD) method [360, 776] and later been extended to coupled-pair methods [454]. In the case of CCD, orbital optimization requires the solution of so-called \(\Lambda\) (or Z vector) equations [744]. There is, however, a cheaper alternative approximating the Z vector by a simple analytical formula [457].
Furthermore, the parametrized coupled-cluster (pCCSD) method of
Huntington and Nooijen [405], which combines the
accuracy of coupled-pair type methods for (usually superior to CCSD, at
least for energies and energy differences) with the higher stability of
the coupled-cluster methods, is an attractive alternative. Comprehensive
numerical tests [404] indicate that particularly
pCCSD(-1,1,1) (or pCCSD/1a) and pCCSD (-1.5,1,1) (or pCCSD/2a) have
great potential for accurate computational thermochemistry. These
methods can be employed by adding the “simple” keywords pCCSD/1a or
pCCSD/2a to the first line of input. As mentioned in section
Local Coupled Pair and Coupled-Cluster Calculations, the LPNO
variants of the pCCSD methods are also available for RHF and UHF
references via usage of the simple keywords LPNO-pCCSD/1a and
LPNO-pCCSD/2a.
7.14.2. Closed-Shell Equations¶
Proceeding from spin-orbitals to the spatial orbitals of a closed-shell determinant leads to the actual working equations of this work. Saebo, Meyer and Pulay have exploited the generator state formalism to arrive at a set of highly efficient equations for the CISD problem [704]. A similar set of matrix formulated equations for the CCSD and QCISD cases has been discussed by Werner and co-workers [355] and the MOLPRO implementation is widely recognized to be particularly efficient. Equivalent explicit equations for the CISD and CCSD methods were published by Scuseria et al. [777][2] The doubles equations for the residual “vector” \(\mathbf{\sigma}\) are (\(i\leqslant j\), all \(a,b\)):
The singles equations are:
The following definitions apply:
The two-electron integrals are written in (11|12) notation and F is the closed-shell Fock operator with F \(^{V}\)being its virtual sub-block. We do not assume the validity of Brillouin’s theorem. The amplitudes \(C_{a}^{i} ,\,C_{ab}^{ij}\) have been collected in vectors \(\mathrm{\mathbf{C} }^{i}\) and matrices \(\mathrm{\mathbf{C} }^{ij}\) wherever appropriate. The shifts \(\Delta^{i}\) and \(\Delta^{ij}\) are dependent on the method used and are defined in Table Table 7.17 for each method implemented in ORCA.
Method |
Doubles Shift |
Singles Shift |
|---|---|---|
CISD |
\(E_{\text{C} }\) |
\(E_{\mathbf{C} }\) |
CEPA/0 |
0 |
0 |
CEPA/1 |
\(\frac{1}{2}(\varepsilon_{i} +\varepsilon_{j} )+\frac{1}{4}\sum\limits_k { (\varepsilon_{ik} +\varepsilon_{jk} ) }\) |
\(\frac{1}{2}\varepsilon_{ii} +\frac{1}{2}\sum\limits_k { \varepsilon_{ik} }\) |
CEPA/2 |
\(\delta_{ij} \varepsilon_{i} +\varepsilon_{ij}\) |
\(\varepsilon_{i} +\varepsilon_{ii}\) |
CEPA/3 |
\((\varepsilon_{i} +\varepsilon_{j} )-\delta_{ij} \varepsilon_{i} -\varepsilon_{ij} +\frac{1}{2}\sum\limits_k { (\varepsilon_{ik} +\varepsilon_{jk} ) }\) |
\(\varepsilon_{i} +\sum\limits_k { \varepsilon_{ik} }\) |
NCEPA/1 |
\(\frac{1}{4}\sum\limits_k { (\varepsilon_{ik} +\varepsilon_{jk} ) }\) |
\(\varepsilon_{ii} +\sum\limits_k { \varepsilon_{ik} }\) |
NCEPA/2 |
\(\varepsilon_{ij}\) |
\(2\varepsilon_{ii}\) |
NCEPA/3 |
\(-\varepsilon_{ij} +\frac{1}{2}\sum\limits_k { (\varepsilon_{ik} +\varepsilon_{jk} ) }\) |
\(2\sum\limits_k { \varepsilon_{ik} }\) |
CPF/1 |
\(N_{ij} \left\{{ \frac{1}{2}(\frac{\varepsilon_{i} }{N_{i} }+\frac{\varepsilon_{j} }{N_{j} })+\frac{1}{4}\sum\limits_k { (\frac{\varepsilon_{ik} }{N_{ik} }+\frac{\varepsilon_{jk} }{N_{jk} }) } } \right\}\) |
\(N_{i} \left\{{ \frac{1}{2}\frac{\varepsilon_{ii} }{N_{ii} }+\frac{1}{2}\sum\limits_k { \frac{\varepsilon_{ik} }{N_{ik} } } } \right\}\) |
CPF/2 |
\(N_{ij} \left\{{ \delta_{ij} \frac{\varepsilon_{i} }{N_{i} }+\frac{\varepsilon_{ij} }{N_{ij} } } \right\}\) |
\(N_{i} \left\{{ \frac{\varepsilon_{i} }{N_{i} }+\frac{\varepsilon_{ii} }{N_{ii} } } \right\}\) |
CPF/3 |
\(N_{ij} \left\{{ \frac{\varepsilon_{i} }{N_{i} }\left({ 1-\delta_{ij} } \right)+\frac{\varepsilon_{j} }{N_{j} }-\frac{\varepsilon_{ij} }{N_{ij} }+\frac{1}{2}\sum\limits_k { (\frac{\varepsilon_{ik} }{N_{ik} }+\frac{\varepsilon_{jk} }{N_{jk} }) } } \right\}\) |
\(N_{i} \left\{{ \frac{\varepsilon_{i} }{N_{i} }+\sum\limits_k { \frac{\varepsilon_{ik} }{N_{ik} } } } \right\}\) |
NCPF/1 |
\(\frac{1}{4}N_{ij} \sum\limits_k { (\frac{\varepsilon_{ik} }{N_{ik} }+\frac{\varepsilon_{jk} }{N_{jk} }) }\) |
\(N_{i} \left\{{ \frac{\varepsilon_{ii} }{N_{ii} }+\sum\limits_k { \frac{\varepsilon_{ik} }{N_{ik} } } } \right\}\) |
NCPF/2 |
\(N_{ij} \frac{\varepsilon_{ij} }{N_{ij} }\) |
\(2N_{i} \frac{\varepsilon_{ii} }{N_{ii} }\) |
NCPF/3 |
\(N_{ij} \left\{{ -\frac{\varepsilon_{ij} }{N_{ij} }+\frac{1}{2}\sum\limits_k { (\frac{\varepsilon_{ik} }{N_{ik} }+\frac{\varepsilon_{jk} }{N_{jk} }) } } \right\}\) |
\(2N_{i} \sum\limits_k { \frac{\varepsilon_{ik} }{N_{ik} } }\) |
ACPF |
\(\frac{2}{n}E_{\text{C} }\) |
\(\frac{2}{n}E_{\text{C} }\) |
ACPF/2 |
\(\frac{2}{n}E_{\text{C} }\) |
\(\left[{ 1-\frac{(n-3)(n-2) }{n(n-1) }} \right]E_{\text{C} }\) |
NACPF |
\(\frac{2}{n}E_{\text{C} }\) |
\(\frac{4}{n}E_{\text{C} }\) |
AQCC |
\(\left[{ 1-\frac{(n-3)(n-2) }{n(n-1) }} \right]E_{\text{C} }\) |
\(\left[{ 1-\frac{(n-3)(n-2) }{n(n-1) }} \right]E_{\text{C} }\) |
The QCISD method requires some slight modifications. We found it most convenient to think about the effect of the nonlinear terms as a “dressing” of the integrals occurring in equations (7.133) and (7.134). This attitude is close to the recent arguments of Heully and Malrieu and may even open interesting new routes towards the calculation of excited states and the incorporation of connected triple excitations.[391] The dressed integrals are given by:
The CCSD method can be written in a similar way but requires 15 additional terms that we do not document here. They may be taken conveniently from our paper about the LPNO-CCSD method [616].
A somewhat subtle point concerns the definition of the shifts in making the transition from spin-orbitals to spatial orbitals. For example, the CEPA/2 shift becomes in the generator state formalism:
(\(\tilde{{\Phi } }_{ij}^{ab}\) is a contravariant configuration state function, see Pulay et al. [744]. The parallel and antiparallel spin pair energies are given by:
This formulation would maintain the exact equivalence of an orbital and a spin-orbital based code. Only in the (unrealistic) case that the parallel and antiparallel pair correlation energies are equal the CEPA/2 shift of Table 7.17 arise. However, we have not found it possible to maintain the same equivalence for the CPF method since the electron pairs defined by the generator state formalism are a combination of parallel and antiparallel spin pairs. In order to maintain the maximum degree of internal consistency we have therefore decided to follow the proposal of Ahlrichs and co-workers and use the topological matrix \(T_{PQ}\) in equation (7.128) and the equivalents thereof in the CEPA and CPF methods that we have programmed.
7.14.3. Open-Shell Equations¶
We have used a non-redundant set of three spin cases (\(\alpha \alpha\), \(\beta \beta\), \(\alpha \beta\)) for which the doubles amplitudes are optimized separately. The equations in the spin-unrestricted formalism are straightforwardly obtained from the corresponding spin orbital equations by integrating out the spin. For implementing the unrestricted QCISD and CCSD method, we applied the same strategy (dressed integrals) as in the spin-restricted case. The resulting equations are quite cumbersome and will not be shown here explicitly [361].
Note that the definitions of the spin-unrestricted CEPA shifts differ from those of the spin-restricted formalism described above (see Kollmar et al. [361]). Therefore, except for CEPA/1 and VCEPA/1 (and of course CEPA/0), for which the spin-adaptation of the shift can be done in a consistent way, CEPA calculations of closed-shell molecules yield slightly different energies for the spin-restricted and spin-unrestricted versions. Since variant 1 is also the most accurate among the various CEPA variants [617], we recommend to use variant 1 for coupled-pair type calculations. For the variants 2 and 3, reaction energies of reactions involving closed-shell and open-shell molecules simultaneously should be calculated using the spin-unrestricted versions only.
A subtle point for open-shell correlation methods is the choice of the reference determinant [618]. Single reference correlation methods only yield reliable results if the reference determinant already provides a good description of the systems electronic structure. However, an UHF reference wavefunction suffers from spin-contamination which can spoil the results and lead to convergence problems. This can be avoided if quasi-restricted orbitals (QROs) are used [391, 611] since the corresponding zeroth-order wavefunction is an eigenfunction of the \(\hat{{S} }^{2}\) operator and thus, no severe spin-contamination will appear. The coupled-pair and coupled-cluster equations will be still solved in a spin-unrestricted formalism but the energy will be slightly higher compared to the results obtained with a spin-polarized UHF reference determinant. Furthermore, especially for more difficult systems like e.g. transition metal complexes, it is often advantageous to use Kohn-Sham (KS) orbitals instead of HF orbitals.
7.14.4. Local correlation¶
As described in previous sections of the manual, ORCA features the extremely powerful LPNO and DLPNO methods. “LPNO” stands for “local pair natural orbital” approximation and DLPNO for “domain based LPNO”. These methods are designed to provide results as close as possible to the canonical coupled-cluster results while gaining orders of magnitude of efficiency through a series of well-controlled approximations. In fact, typically about 99.8% to 99.9% of the canonical correlation energy is recovered in such calculations. Even higher accuracy is achievable but will ultimately come at much higher computational cost. The default cut-offs are set such that the vast majority of chemically meaningful energy differences are reproduced to better than 1 kcal/mol relative to the canonical results with the same basis set. Of the LPNO and DLPNO methods, the LPNO is the older one and will eventually be discarded from ORCA. It has some higher order scaling steps (up to \(N^5\)) while DLPNO is linear scaling and of similar accuracy. Amongst the DLPNO methods, the first generation implementation of the DLPNO methods (DLPNO2013) is only near-linear scaling, while the DLPNO implementation since ORCA 4.0 is linear scaling.
It is important to understand that the LPNO and DLPNO implementations are intimately tied to the RI approximation. Hence, in these calculations one must specify a fitting basis set. The same rules as for RI-MP2 apply: the auxiliary basis sets optimized for MP2 are just fine for PNO calculations. You can specify several aux bases for the same job and the program will sort it out correctly.
The theory of the LPNO methods has been thoroughly described in the literature in a number of original research papers.[616, 622, 720, 722]
Hence, it is sufficient to only point to a few significant design principles and features of these methods:
The correlation energy of any molecule can be written as a sum over the correlation energies of pairs of electrons, labelled by the internal indices (\(ij\)) of pairs of orbitals that are occupied in the reference determinant. If the orbitals (\(i\)) and (\(j\)) are localized, the pair correlation energy \(\epsilon_{ij}\) falls off very quickly with distance, quite typically by about an order of magnitude per chemical bond that is separating orbitals (\(i\)) and (\(j\)). Hence, by using a suitable cut-off for a reasonable pair correlation energy estimate many electron pairs can be removed from the high-level treatment and only a linear scaling number of electron pairs will make a significant contribution to the correlation energy.
The natural estimate for the pair correlation energy comes from second order many body perturbation theory (MP2). However, canonical MP2 is scaling with the fifth power of the molecular size and hence, is not really attractive from a theoretical nor computational point of view. Owing to the small pre-factor RI-MP2 goes a long way to provide reasonably cheap estimates for pair correlation energies. However, if one uses localized internal orbitals, then the MP2 energy expression must be cast in form of linear equations. On the other hand, if one uses canonical virtual orbitals together with localized internal orbitals and neglects the coupling terms coming from purely internal Fock matrix elements F(\(i\),\(k\)) and F(\(k\),\(j\)) then one ends up with a fair approximation that is termed “semi-canonical MP2” in ORCA. It serves as a guess in the older LPNO method. For DLPNO this method is also not attractive.
In DLPNO, the guess is more sophisticated. Here the virtual space is spanned by projected atomic orbitals (PAOs) that are assigned to domains of atoms that are associated with each local internal orbital (\(i\)) and with the union of two such domains (\(i\)) and (\(j\)) for the electron pair (\(ij\)). If one applies the semi-local approximation, one obtains an excellent approximation to the semi-canonical MP2 energy. This is called the “semi-local” approximation and it scales linearly with respect to computational effort. If one iterates these equations to self-consistency to eliminate the coupling terms F(\(i\),\(k\)) and F(\(k\),\(j\)) then one obtains the full local MP2 method (LMP2 or L-MP2). By making the domains large enough the results approach the canonical MP2 energy to arbitrary accuracy while still being linear scaling with respect to computational resources. This method is the default for the DLPNO method.
Basically, the high-spin open-shell version of the DLPNO uses the same strategy as the closed-shell variant to efficiently generate the open-shell PNOs in a consistent manner to the closed-shell formalism. Through the development of the UHF-LPNO-CCSD method, we have realized that use of the unrestricted MP2 (UMP2) approach to define the open-shell PNOs introduces a few complexities; (1) the PNOs for \(\beta\) spin orbitals cannot be defined for \(\alpha\)-\(\alpha\) electron pairs and vice versa, (2) the diagonal PNOs for singly occupied orbitals cannot be properly defined, and (3) the PNO space does not become identical to that in the closed-shell LPNO framework in the closed-shell limit. However, to program all the unrestricted CCSD terms in the LPNO basis, those PNOs are certainly necessary. Therefore, in the UHF-LPNO-CCSD implementation, several terms, which, in many cases, give rather minor contributions in the correlation energy are omitted. Due to these facts, the UHF-LPNO-CCSD does not give identical results to that of RHF-LPNO-CCSD in the closed-shell limit. In addition, screening of the weak pairs on the basis of the semi-canonical UMP2 pair-energy results in somewhat unbalanced treatment of the closed- and open-shell states in some cases, leading to rather larger errors in the reaction energies. To overcome those issues, in the high-spin open-shell DLPNO-CCSD method, the PNOs are generated in the framework of semi-canonical NEVPT2 which smoothly converges to the RHF-MP2 counterpart in the closed-shell limit. The open-shell DLPNO-CCSD, which is made available from ORCA 4.0, includes a full set of the open-shell CCSD equations and is designed as a natural extension to the RHF-DLPNO-CCSD.
Screening of the electron pairs according to a truncation parameter (in ORCA it is called \(T_\mathrm{CutPairs}\)) is not sufficient to obtain a highly performing local electronic structure method. The original work of Pulay suggested to limit excitations out of the internal orbitals (\(i\)) and (\(j\)) to the domain associated with the pair (\(ij\)). While this works well and has been implemented to perfection by Werner, Schütz and co-workers over the years,[754, 755, 774, 775] the pre-factor of such calculations is high, since the domains have to be chosen large in order to recover 99.9% or more of the canonical correlation energy.
The ORCA developers have therefore turned to an approach that has been used with a high degree of success in the early 1970s by Meyer, Kutzelnigg, Staemmler and their co-workers, namely the method of “pair natural orbitals” (PNOs).[10, 582, 583, 886]
As shown by Loewdin in his seminal paper from 1955, natural orbitals (the eigenfunctions of the one-particle density matrix) provide the fastest possible convergence of the correlated wavefunction with respect to the number of one-particle functions included in the virtual space. It has been amply established that approximate natural orbitals are almost as succesful as the true natural orbitals (which would require the knowledge of the exact wavefunction) in this respect. While the success of approximate correlation treatments of many particle systems that use approximate natural orbitals of the whole systems are somewhat limited, this is not the case for pair natural orbitals. The latter have first been suggested as a basis for correlation calculations by England and co-workers and, at the time, were given the name “pseudonatural orbitals”, a term that was used by Meyer throughout his pioneering work.The PNOs are approximate natural orbitals of a given electron pair. In order to generate them one requires a one particle pair density matrix \(D_{ij}\) the eigenfunctions of which are the PNOs themselves while the corresponding eigenvalues are the PNO occupation numbers. While there are many creative possibilities that can lead to slightly different PNO sets, a quite useful and natural approximation is to generate such a density from the MP2 amplitudes as an expectation value (the “unrelaxed” MP2 density. One then uses a second threshold (in ORCA\(T_\mathrm{CutPNO}\)) that controls the PNOs to be included in the calculation. PNOs with an occupation number \(< T_\mathrm{CutPNO}\) are neglected.
PNOs of a given electron pair form an orthonormal set. However, PNOs belonging to different electron pairs are not orthonormal and hence they overlap. This non-orthogonality leads to surprisingly few complications because the PNOs stay orthogonal to all occupied orbitals. In practice, the equations for PNO-based correlation calculations are hardly more complex than the canonical equations.
The nice feature of these pair densities is that they become small when the pair interaction becomes small. Hence weak pairs are correlated by very few PNOs. Therefore, the PNO expansion of the wavefunction is extremely compact and there only is a linear scaling number of significant excitation amplitudes that need to be considered.
A great feature of the PNOs is that they are “self-adapting” to the correlation situation that they are supposed to describe. Hence, they are as delocalized as required by the physics and there is no ad-hoc assumption about their location in space. However, it is clear that the PNOs are located in the same region of space as the internal orbitals that they correlate because otherwise they would not contribute to the correlation energy.
In the iterative local MP2, a set of PNOs is calculated for the MP2 calculation from the semicanonical amplitudes first using the cutoff
TCutPNO\(\times\)LMP2ScaleTCutPNO. The size of the resulting PNO space should be comparable with DLPNO-MP2. After the iterations have converged, a (smaller) set of PNOs for the CCSD double excitation amplitudes is generated from the iterated local MP2 amplitudes in the (larger) PNO basis. The PNOs for the single excitation amplitudes are always calculated using the semicanonical MP2 amplitudes, as they require a much more conservative PNO truncation threshold.Capitalizing on this feature one can define generous domains and expand the PNOs in terms of the PAOs and auxiliary fit functions (for the RI approximations) that are contained in these domains. In ORCA this is controlled by the third significant truncation parameter \(T_\mathrm{CutMKN}\). This is the basis of the DLPNO method. In the older LPNO method, the PNOs are expanded in terms of the canonical virtual orbitals and \(T_\mathrm{CutMKN}\) is only used for a local fit to the PNOs. In the linear-scaling DLPNO implementation the domains expanding the PNOs in terms of the PAOs are controlled via \(T_\mathrm{CutDO}\).
PNO expansions have the amazing advantage that the PNOs converge to a well-defined set as the basis set is approaching completeness. Hence, the increase in the number of PNOs per electron pair is very small upon enlargement of the orbital basis. In other words, correlation calculations with large basis sets are not that much more expensive than correlation calculations with small basis sets. Thus, the advantage of PNO over canonical calculations increases with the size of the basis set. This is a great feature as large and flexible basis sets are a requirement for meaningful correlation calculations.
In summary, DLPNO and LPNO calculations are controlled by only three cut-off parameters with well-defined meanings: a) \(T_\mathrm{CutPairs}\), the cut-off for the electron pairs to be included in the coupled-pair or coupled-cluster iterations, b) \(T_\mathrm{CutPNO}\) which controls how many PNOs are retained for a given electron pair and c) \(T_\mathrm{CutMKN}\) that controls how large the domains are that the PNOs expand over. For the linear-scaling DLPNO calculations the domain sizes are controlled via \(T_\mathrm{CutDO}\).
It is clear that owing to the truncations a certain amount of error is introduced in the results. However, having the LMP2 results available, one can compensate for the errors coming from \(T_\mathrm{CutPairs}\) and \(T_\mathrm{CutPNO}\). This is done in ORCA and the correction is added to the final correlation energy, thus bringing it very close (to mEh accuracy or better) to the canonical result. \(T_\mathrm{CutMKN}\) is best dealt with by making it conservative (at 1e\(-\)3 to 1e\(-\)4 the domains are about 20–30 atoms large, which is sufficient for an accurate treatment).
Note that the LPNO and DLPNO methods do not introduce any real space cut-offs. We refrain from doing so and insist in our method development by making all truncations based on wavefunction or energy parameters. We feel that this is the most unbiased approach and it involves no element of “chemical intuition” or “prejudice”.
In the DLPNO method a highly efficient screening mechanism is operative. In this method one first obtains a (quadratically scaling) multipole estimate for the pair correlation energy that is extremely fast to compute (a few seconds, even for entire proteins). Only if this estimate is large enough, a given electron pair is even considered for a LMP2 treatment. Quite typically pairs with energy contributions of 1e\(-\)6 Eh and smaller are very well described by the dipole-dipole estimate. Specifically, we drop pairs with estimated energies \(< 0.01 \times T_\mathrm{CutPairs}\) and add their multipole energy sum to the final correlation energy. These corrections tend to be extremely small, even for large molecules and are insignificant for the energy. However, importantly, the multipole estimate ensures linear scaling in the MP2 treatment. The pairs that then do not survive the pair-prescreening are called “weak pairs” in the ORCA or DLPNO sense. They still play a role in the calculation of the triple excitation correction.
The calculation of triple excitation contributions is more involved and one does not have a perturbative estimate available since the (T) contribution is perturbative itself. While the (T) contribution is much smaller than the CCSD correlation energy, the errors introduced by the various local and PNO approximations can be significant. We found that one has to include triples with at least one pair being a “weak” LMP2 pair (with its LMP2 amplitudes) in order to arrive at sufficiently accurate results.
Given these explanations the various cut-off parameters that can be controlled in LPNO and DLPNO calculations should be understandable and are listed below. We emphasize again that only the three thresholds \(T_\mathrm{CutPairs}\), \(T_\mathrm{CutPNO}\) and \(T_\mathrm{CutMKN}\) should be touched by the user, unless very specific questions are addressed. The recommended way to control the accuracy of calculations is to specify “TightPNO”, “NormalPNO” or “LoosePNO” keywords, rather than to change numeric values of cutoffs. Individual thresholds should normally not be changed, as the defaults are sensible and lead to good cost/performance ratios.
%mdci TCutPairs 1e-4 # cut-off for the pair truncation
TCutPNO 3.33e-7 # cut-off for the PNO truncation
TCutDO 1e-2 # cut-off for the DLPNO domain construction
TCutMKN 1e-3 # cut-off for the local fit
# for DLPNO2013: also domain construction
# remaining options, tied to the three main cut-offs,
EXPERTS ONLY!
Localize true # flag for using localized orbitals
LocMet AHFB # Localization method.
# Options: PM, FB, IAOIBO, IAOBOYS, NEWBOYS, AHFB
LocTol 1e-6 # Absolute threshold for the localization procedure
# Automatically adjusted by default.
LocTolRel 1e-8 # Relative threshold for the localization procedure
LocMaxIter 128 # Maximum number of localization iterations
LocRandom 1 # default, take random seed for any localization
# For internal orbitals: choose best of 32 localizations
# Switched off for AHFB
0 # take constant seed for any localization (for testing)
LocNAttempts np # number of localization attempts
# default: number of processes, minimum 8, if
# randomize true
# 1, if randomize false
# any number larger or equal np, if randomize true
PNONorm MP2Norm # default, old IEPANorm can also be used
(near identical results)
NrMP2Pairs_Trip 1 # number of MP2 pairs to be included in the triples
calculation
PAOOverlapThresh 1e-8 # generation of non-redundant PAOs from redundant ones
UseFullLMP2Guess false # Use iterative full LMP2 (for DLPNO)
SinglesFockUsePNOs true # compute the Singles Fock matrix (SFM) in PNOs.
# DLPNO2013: default for SinglesFockUsePNOs is false,
# by default RIJCOSX is used for the SFM, except when
# RIJK is given. In that case the RIJK-SFM is used.
LMP2MaxIter 50 # max no of iterations in the MP2 equations
LMP2TolE 1e-7 # LMP2 energy convergence tolerance
LMP2TolR 5e-7 # LMP2 residual convergence tolerance
LMP2ScaleTCutPNO # PNO cutoff for LMP2 is: TCutPNO*LMP2ScaleTCutPNO
# Default: TCutPNO(DLPNO-MP2)/TCutPNO(DLPNO-CCSD) with
# respective TCutPNOs specific to Loose/Normal/TightPNO
LMP2FCut 1e-5 # LMP2 neglect cut-off for
off-diagonal Fock matrix elements
TCutTNO 0 # Cut-off for triples natural orbitals (0=automatic)
TCutPNOSingles -1 # -1= use 0.03*TCutPNO
TCutPreScr -1 # -1= use 0.01*TCutPairs for
multipole estimate based screening
TCutDeloc 0.1 # delocalization threshold for specification
of extended domains.
Necessary because PAOs are not strictly localized
TCutOSV 1e-6 # orbital specific virtuals used for pre-screening.
No critical
TScaleMP2Pairs 0.1
TScaleMKNStrong 10
TScaleMKNWeak 100
For larger systems and tighter thresholds the disk I/O of a DLPNO calculation may become challenging. In this case, it might be usefull to keep some integrals in memory, if enough RAM is available. With the following flag
%mdci
StorageType Shared
end
ORCA will try to store certain much-used integrals in shared memory. If the amount of memory is not sufficient, ORCA will fall back to on-disk storage. NOTE: This flag will only work if all processes work on the same node.
IMPORTANT NOTE REGARDING ORBITAL LOCALIZATION
Localized orbitals for DLPNO are obtained via the Foster-Boys method with an augmented Hessian converger (
AHFB) by default.The localization tolerance (
LocTol) is coupled to the SCF gradient tolerance (TolG) with a constant factor by default. Selecting specific SCF convergence settings (such asTightSCF) therefore also ensures obtaining a set of appropriately well converged localized orbitals. This can be overridden by setting a different value forLocTol.An important feature of the augmented Hessian converger is that it systematically approaches a local maximum of the localization function (even though convergence to the global maximum cannot be guaranteed). As opposed to that, the conventional localization method (
FB) may stop, for example, at a saddle point. In bad cases, this can lead to deviations of several kJ/mol in the DLPNO energy. Likewise, it can contribute towards lack of reproducibility of results.No similar procedure has been implemented for the other localization methods (such as Pipek-Mezey) yet. The same problems as with the FB converger can occur in these cases.
No randomization is used for the
AHFBconverger.
The old default orbital localization settings of ORCA 4.0 can be reproduced with the following options:
%MDCI
LocMet FB
LocTol 1.0e-6
LocRandom 1
End
Regarding the methods that employ randomization (FB, PM, IAOIBO,
IAOBOYS) only, the following notes apply:
Generally, better DLPNO results are obtained when several runs of localization are undertaken using different initial guesses. The different initial guesses are obtained using randomization (LocRandom).
However, randomization of the initial guess can lead to differently localized MOs in different calculations. This can yield non-identical correlation energies, varying in the sub-kJ/mol range, for different runs on the same machine.
In order to yield identical correlation energy results, randomization can be switched off (LocRandom 0). However, switching off randomization only leads to identical results on the same machine, but can still lead to slightly different results (in the sub-kJ/mol range) on different machines.
Reproducibility of the correlation energy is expected to increase further if LocNAttempts is set to higher values.
The input below shows how to perform a DLPNO calculation with settings that exactly reproduce the canonical RI-MP2 result. They are not recommended for production use, but merely to show that if the local approximations are pushed, then the result coincides with the canonical one. If one would set \(T_\mathrm{CutPNO}\) to zero this would give canonical RI-CCSD. However, this is an absurdly inefficient calculation and hence not done.
#
! def2-SV(P) def2/JK def2-SVP/C RI-JK DLPNO-CCSD VeryTightSCF RI-MP2
# obtain a result that only contains errors from the PNO approximation
# but no others
%mdci TCutPairs 0
TCutMKN 0
UseFullLMP2Guess true
LMP2FCut 1e-9
LMP2MaxIter 25
LMP2TolE 1e-10
LMP2TolR 1e-11
PAOOVerlapThresh 1e-9
end
! Bohrs
* xyz 0 1
C -1.505246952209632 1.048213673267046 -3.005665895986369
C 1.289678561934891 0.246429688933291 -3.259735682020124
C 2.834670835163566 1.157307360133605 -0.990383454919828
C 1.924119415395082 -0.128330938291771 1.465070676514038
C -0.931529472233802 -0.722841293992075 1.397639867298547
C -2.347670084056626 1.213332291655600 -0.217984867773136
H 2.084955694093313 0.973408301535989 -5.037750251258102
H 1.426532559234904 -1.831017720289521 -3.371063003813707
H -1.795307501459984 2.891278294563413 -3.927855043896308
H -2.709613973668925 -0.308515546176734 -4.026627646697411
H -4.404246093821399 0.941639912907262 -0.071175054238094
H -1.962867323232915 3.122079490952855 0.528101313545138
H -1.245096579039474 -2.621186110634707 0.594784162223769
H -1.699155144887690 -0.782162821007662 3.328959985756973
H 2.347109421287126 1.104305785540561 3.087624818244846
H 2.990679065503112 -1.888017241218143 1.775287572161196
H 4.862301668284708 0.796425411350593 -1.279131939569907
H 2.634027658640572 3.226752635113244 -0.827936424652650
*
7.14.4.1. Including (semi)core orbitals in the correlation treatment¶
In some chemical applications some or all of the chemical (semi)core electrons must be included in the correlation treatment. In this case, it is necessary to tighten the TCutPNO thresholds for electron pairs in which chemical (semi)core electrons are involved. This is now the default in DLPNO calculations.
For instance, one can decide to switch off the frozen-core approximation and include all the electrons in the correlation treatment. In this case, the program will use tighter thresholds by default for all electron pairs and Singles that involve chemical core electrons. Note that, in this case, the use of properly optimized basis functions for correlating the inner electrons is highly recommened.
! DLPNO-CCSD(T) def2-SVP def2-SVP/C NoFrozenCore
%mdci TSCALEPNO_CORE 0.01 # scaling factor for TCutPNO for electron pairs and
# Singles involving chemical core electrons
end
* xyz 0 1
Ti 0.0001595288 0.0000775041 0.0000000000
F 1.7595996122 0.0000634675 -0.0000000011
F -0.5865076471 1.6586935196 0.0000000018
F -0.5866248292 -0.8294172469 -1.4362516915
F -0.5866248311 -0.8294172443 1.4362516907
*
Another option is to choose the involved chemical core electrons by using an energy window. In this way all electron pairs and Singles that involve chemical core electrons, which are in the defined energy window, are affected by TScalePNO_CORE.
! DLPNO-CCSD(T) def2-SVP def2-SVP/C
%method
FrozenCore FC_EWIN
end
%mdci
EWIN -40, 10000
end
* xyz 0 1
Ti 0.0001595288 0.0000775041 0.0000000000
F 1.7595996122 0.0000634675 -0.0000000011
F -0.5865076471 1.6586935196 0.0000000018
F -0.5866248292 -0.8294172469 -1.4362516915
F -0.5866248311 -0.8294172443 1.4362516907
*
A summary with the number of electrons affected by TScalePNO_Core for the examples just discussed is shown in Table Table 7.18.
Keyword |
Chemical Core Electrons |
Valence Electrons |
|
Frozen |
Included\(^{a}\) |
||
FrozenCore (default) |
18 |
0 |
40 |
NoFrozenCore |
0 |
18 |
40 |
EWIN -40, 10000 |
16 |
2 |
40 |
\(^{a}\) using TScalePNO_Core.
By default, ORCA provides a chemical meaningful definition for the number of electrons which belong to the chemical core of each element. As already discussed, these default values define which pairs are affected by TScalePNO_Core. However, the user can modify the number of chemical core electrons for a specific element via the NewNCore keyword.
! DLPNO-CCSD(T) def2-SVP def2-SVP/C NoFrozenCore
%method
NewNCore Ti 8 end
end
* xyz 0 1
Ti 0.0001595288 0.0000775041 0.0000000000
F 1.7595996122 0.0000634675 -0.0000000011
F -0.5865076471 1.6586935196 0.0000000018
F -0.5866248292 -0.8294172469 -1.4362516915
F -0.5866248311 -0.8294172443 1.4362516907
*
In the previous example, the number of chemical core electrons for Ti has been fixed to 8.
Starting from ORCA 6.0, in DLPNO calculations, tightened TCutPNO thresholds are used by default for “semicore” electron pairs involving the 3s and 3p orbitals of first-row transition metals.[32] This improves not only the accuracy of the results noticeably but also the efficiency of the computations as the system size grows (ref.). To reproduce results obtained with earlier versions, the number of semicore orbitals on first-row transition metals needs to be set to zero (the old default) instead of eight (the current default), as done below for a Zn atom with the NewNSemiCore keyword in the method block:
%method NewNSemiCore Zn 0 end end
NOTE
Of course, if electrons are replaced by ECPs, they are not included in the correlation treatment.
If ECPs are used, the number for NewNCore has to include the electrons represented by the ECPs as well. E.g. if an element is supposed to have 60 electrons in the ECP and additional 8 electrons should be frozen in the correlation calculation, NewNCore should be 68.
The different sets of orbitals (chemical core electrons included in the correlation treatment and valence electrons) are localized separately in order to avoid the mixing of core and valence orbitals.
7.14.4.2. Multi-Level Calculations¶
In many applications events are investigated that are located in a relatively small region of the system of interest. In these cases, combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have been proved to be extremely useful, especially in the modeling of enzymatic reactions. The basic idea of any QM/MM method is to treat a small region of the system at the QM level and to use an MM treatment for the remaining part of the system. Alternatively, QM/QM methods, where different parts of a system are studied at various quantum mechanical levels, are also available. Quantum mechanical methods are more computationally demanding than the molecular mechanics treatment, and this limits the applicability of all-QM approaches. Nevertheless, QM/QM methods retain some strong advantages over QM/MM schemes. For instance, force field parameters for the molecular mechanics part of the calculation are not necessary, and thus there are no restrictions on the type of chemical systems that can be treated. Moreover, problems usually caused by boundaries between QM and MM parts do not occur. Finally, the accuracy of an all-QM calculation is expected to be higher compared to the accuracy of QM/MM approaches, that is limited by the MM treatment.
In ORCA, the different methods that can be used in a QM/QM calculations are:
DLPNO-CCSD(T) with TightPNO thresholds
DLPNO-CCSD(T) with NormalPNO thresholds
DLPNO-CCSD(T) with LoosePNO thresholds
DLPNO-CCSD
DLPNO-MP2
HF
The user can define an arbitrary number of fragments in the input, the level of theory to be used for each fragment and for the interaction between different fragments. Localized molecular orbitals are then assigned to a given fragment on the basis of their Mulliken populations.
The following example shows the calculation of a benzene dimer, for which the individual monomers are calculated on MP2 level, and the interaction between the two monomers is calculated on TightPNO DLPNO-CCSD(T) level. More realistic use cases are discussed in ref. [811].
! DLPNO-CCSD(T) Def2-SVP Def2-SVP/C
%mdci
UseFullLmp2Guess True
TightPNOFragInter {1 2}
MP2FragInter {1 1} {2 2}
end
*xyz 0 1
C(1) 1.393 0.000 0.0
H(1) 2.474 0.000 0.0
C(1) 0.695 1.206 0.0
H(1) 1.238 2.143 0.0
C(1) -0.695 1.206 0.0
H(1) -1.238 2.143 0.0
C(1) -1.393 0.000 0.0
H(1) -2.474 0.000 0.0
C(1) -0.695 -1.206 0.0
H(1) -1.238 -2.143 0.0
C(1) 0.695 -1.206 0.0
H(1) 1.238 -2.143 0.0
C(2) 2.333 1.33 3.5
H(2) 3.414 1.33 3.5
C(2) 1.635 2.536 3.5
H(2) 2.178 3.473 3.5
C(2) 0.245 2.536 3.5
H(2) -0.298 3.473 3.5
C(2) -0.453 1.33 3.5
H(2) -1.534 1.33 3.5
C(2) 0.245 0.124 3.5
H(2) -0.298 -0.813 3.5
C(2) 1.635 0.124 3.5
H(2) 2.178 -0.813 3.5
*
For the calculation of the interaction energy, the energy of the individual benzene monomer should be calculated on the accuracy level of the monomer in the dimer calculation, i.e. using MP2 with full LMP2 guess for the above example.
All possible settings for the multi-level calculation are listed below.
# The one-keyword line defines the default method for the multi-level calculation.
# Options here are DLPNO-CCSD(T) or DLPNO-CCSD with the addition of the
# LoosePNO, NormalPNO and TightPNO keyword
!DLPNO-CCSD(T)
# The below given keywords define the changes with respect to the
# above given default method. The user should take care that each intra- or
# interfragment combination is defined only once (unlike in the example given below)
%mdci
LoosePNOFragInter {1 1} {2 2} # use LoosePNO settings for the intrafragment
# pair energies of fragments 1 and 2
NormalPNOFragInter {1 2} # use NormalPNO settings for the interfragment
# pair energies between fragment 1 and 2
TightPNOFragInter {1 3} # use TightPNO settings for the interfragment
# pair energies between fragment 1 and 3
NormalPNOTightPairFragInter {1 2} # use NormalPNO settings but with TCutPairs
# 1.e-5 for the interfragment pair energies
# between fragment 1 and 2
LoosePNOTightPairFragInter {1 3} # use LoosePNO settings but with TCutPairs 1.e-5
# for the interfragment pair energies between
# fragment 1 and 3
NoTriplesFragments 1, 2 # if all MOs of a triple are located on fragment
# 1 and / or 2, the triple is neglected
MP2FragInter {1 1} {2 2} # compute the intrafragment pair energies of
# fragments 1 and 2 on MP2 level
HFFragInter {1 1} {2 2} # compute the intrafragment energies on HF level
UseFullLmp2Guess false # default = false,
# if MP2FragInter is used: default = true
# The fragments need to be defined in the xyz section.
*xyz 0 1
C(1) 1.393 0.000 0.0
H(1) 2.474 0.000 0.0
...
7.14.4.3. Multi-Level Calculations for IP and EA-EOM-DLPNO-CCSD¶
The multi-layer method can be used to include the environmental effect in IP-and EA-EOM-DLPNO-CCSD method. A typical input file for the multi-layer IP-EOM-CCSD method will look like
! IP-EOM-DLPNO-CCSD TightSCF NORMALPNO def2-SVP def2-SVP/C def2/J RIJCOSX
%mdci
nroots 4
DTOl 1e-7
NORMALPNOFragInter { 1 1}
LOOSEPNOFragInter { 1 1}
HFFRAGMENTINTERACTION { 2 2}
end
*xyz 0 1
C(1) 2.782064 -1.456235 0.000000
C(1) 1.478695 -0.729491 0.000000
C(1) 0.274461 -1.343436 0.000000
N(1) -0.914790 -0.659079 0.000000
C(1) -0.988897 0.709718 0.000000
N(1) 0.241239 1.324758 0.000000
H(1) 0.224165 2.335424 0.000000
C(1) 1.507368 0.726274 0.000000
O(1) 2.518005 1.411594 0.000000
O(1) -2.043648 1.337449 0.000000
H(1) -1.808257 -1.143873 0.000000
H(1) 0.182931 -2.420317 0.000000
H(1) 2.626386 -2.532891 0.000000
H(1) 3.370859 -1.183830 -0.874506
H(1) 3.370859 -1.183830 0.874506
O(2) -3.661424 -0.883408 0.000000
H(2) -3.462053 0.068032 0.000000
H(2) -4.615649 -0.964529 0.000000
*
Here the example is a mono-hydrated thymine molecule, where the thymine is treated at the main fragment and water is treated at the environment. It will result in the following output
----------------------
EOM-CCSD RESULTS (RHS)
----------------------
IROOT= 1: 0.322826 au 8.785 eV 70852.1 cm**-1
Amplitude Excitation
-0.689911 37 -> x
Percentage singles character= 96.89
IROOT= 2: 0.364959 au 9.931 eV 80099.2 cm**-1
Amplitude Excitation
-0.689956 35 -> x
Percentage singles character= 95.49
IROOT= 3: 0.378175 au 10.291 eV 82999.9 cm**-1
Amplitude Excitation
-0.691827 36 -> x
Percentage singles character= 93.93
IROOT= 4: 0.403845 au 10.989 eV 88633.7 cm**-1
Amplitude Excitation
-0.690254 34 -> x
Percentage singles character= 96.55
The result of a full a IP-EOM-DLPNO-CCSD calculation with NORMALPNO setting would have looked like
IROOT= 1: 0.322576 au 8.778 eV 70797.3 cm**-1
Amplitude Excitation
0.689734 37 -> x
Percentage singles character= 96.88
IROOT= 2: 0.364691 au 9.924 eV 80040.3 cm**-1
Amplitude Excitation
-0.689947 35 -> x
Percentage singles character= 95.50
IROOT= 3: 0.377966 au 10.285 eV 82954.0 cm**-1
Amplitude Excitation
-0.691801 36 -> x
Percentage singles character= 93.94
IROOT= 4: 0.402497 au 10.953 eV 88337.9 cm**-1
Amplitude Excitation
-0.690138 34 -> x
Percentage singles character= 96.50
The results in multi-layer IP-EOM-DLPNO-CCSD method has been found to be in excellent agreement with standard variant. The \(MP2FragInter\) treatment is not available for the EOM method. To get a reasonable accuracy one need to treat the fragment from where the ionization is happening (thymine in the above example) at the highest possible level. The interaction between the main fragment (thymine) and environment (water) should be treated at the intermediate level accuracy. The environment can safely be treated with \(HFFragInter\) for almost all the cases. One can decide the size of the main fragment by looking at HF occupied orbitals as the Koopmans’ approximation is a very good zeroth order guess for the IP values. The electron attached states are much less localized as compared to the ionization problem. Consequently, the multi-layer EA-EOM-DLPNO-CCSD requires much more tighter thresholds than the IP variant. An typical input file for multi-layer EA-EOM-DLPNO-CCSD will look as follows.
! EA-EOM-DLPNO-CCSD NORMALPNO ma-def2-SVP RIJCOSX aug-cc-pVDZ/C def2/J
%mdci
NRoots 4
FollowCIS true
TCutPNOSingles 1e-12
MaxIter 2000
DTol 1e-7
NDAV 10
NormalPNOFragInter { 1 1 }
LoosePNOFragInter { 1 2 } { 2 2 }
end
* xyz 0 1
N(1) -1.114 -0.934 -3.554
C(1) -0.343 -0.202 -4.483
H(1) -0.668 -0.311 -5.520
C(1) 0.635 1.107 -2.633
O(1) 1.241 2.018 -2.013
N(1) 0.022 0.050 -1.776
H(1) -0.069 0.339 -0.800
C(1) -0.975 -0.782 -2.233
O(1) -1.697 -1.422 -1.333
C(1) 1.087 1.852 -4.986
H(1) 1.673 2.594 -4.418
H(1) 1.771 1.340 -5.697
H(1) 0.348 2.403 -5.609
H(1) -1.823 -1.635 -3.888
N(2) -4.904 -4.773 -4.958
H(2) -4.634 -3.987 -5.572
C(2) -4.875 -4.201 -3.599
C(2) -3.704 -3.226 -3.310
H(2) -5.818 -3.646 -3.423
H(2) -4.859 -5.012 -2.850
O(2) -3.026 -2.826 -4.301
H(2) -4.078 -5.386 -5.029
O(2) -3.559 -2.899 -2.077
H(2) -2.494 -2.057 -1.754
C(2) 0.440 0.898 -4.018
end
It is a thymine-glycine complex where the thymine is treated as the main fragment and glycine as the environment. One needs to use a more tighter value of \(TCutPNOSingles\) for EA as in the case of standard EA-EOM-DLPNO-CCSD. The \(TCutPNOSingles\) for the respective fragments automatically gets adjusted based on their respective \(TCutPNO\) values. The output will be
----------------------
EOM-CCSD RESULTS (RHS)
----------------------
IROOT= 1: -0.039822 au -1.084 eV -8739.9 cm**-1
Amplitude Excitation
0.689012 x -> 53
Percentage singles character= 93.33
IROOT= 2: 0.025156 au 0.685 eV 5521.0 cm**-1
Amplitude Excitation
-0.614385 x -> 54
0.139410 x -> 55
0.297903 x -> 59
Percentage singles character= 96.86
IROOT= 3: 0.044569 au 1.213 eV 9781.7 cm**-1
Amplitude Excitation
0.116867 x -> 54
0.684240 x -> 55
Percentage singles character= 98.16
IROOT= 4: 0.053677 au 1.461 eV 11780.7 cm**-1
Amplitude Excitation
0.695211 x -> 56
Percentage singles character= 98.12
The results are in excellent agreement with the standard EA-EOM-DLPNO-CCSD method.
----------------------
EOM-CCSD RESULTS (RHS)
----------------------
IROOT= 1: -0.038862 au -1.057 eV -8529.3 cm**-1
Amplitude Excitation
-0.689412 x -> 53
Percentage singles character= 93.47
IROOT= 2: 0.025448 au 0.692 eV 5585.2 cm**-1
Amplitude Excitation
-0.654101 x -> 54
0.131404 x -> 55
0.207630 x -> 59
Percentage singles character= 97.47
IROOT= 3: 0.044651 au 1.215 eV 9799.8 cm**-1
Amplitude Excitation
0.112183 x -> 54
0.684562 x -> 55
Percentage singles character= 98.17
IROOT= 4: 0.053780 au 1.463 eV 11803.3 cm**-1
Amplitude Excitation
0.695635 x -> 56
Percentage singles character= 98.16
To get the reasonable accuracy in multi-layer EA-EOM-CCSD one need to treat the environment and inter-fragment interaction atleast at \(LoosePNOFragInter\) level.
7.14.5. The singles Fock term¶
In most MDCI calculations, there is an intermediate, which resembles closely to the SCF Fock matrix, and similar methods are available to efficiently calculate it. In the followings, a short discussion will be given of the so-called singles Fock term, which in the closed shell case has the form
The singles Fock matrix can be obtained via transformation from its counterpart (\(G(\mathbf{t}_1)_{\mu\nu}\)) in the atomic orbital (AO) basis
where
is the analogue of the SCF density matrix for the singles Fock case. For the singles Coulomb (\(J(\mathbf{t}_1)_{\mu\nu}\)) case, the density may be symmetrized (\(\tilde{P}(\mathbf{t}_1)_{\kappa\tau}=P(\mathbf{t}_1)_{\kappa\tau}+P(\mathbf{t}_1)_{\tau\kappa}\)), and one may use the resolution of identity approximation
where \(A, B\) are elements of the RI/DF auxiliary fitting basis. Note that the factor of 2 in ((7.148)) is taken care of by symmetrization. Since we are using a symmetric density, we may use the same efficient routine to evaluate the singles Coulomb term as in the SCF case, see Using the RI-J Approximation to the Coulomb Part and The Split-RI-J Coulomb Approximation.
For the exchange case (\(K(\mathbf{t}_1)_{\mu\nu}\)), one possibility is to use the COSX approximation (see Using the RI Approximation for Hartree-Fock and Hybrid DFT (RIJCOSX))
The COSX routine is able to deal with asymmetric densities as well, and thus, it can be used here similar to the SCF case.
The other possibility is to use RI for exchange (RIK),
where
and the \(\tilde{j}\) is an “orbital” transformed using \(C(\mathbf{t}_1)\).
Using these approximations, there are two approximations for the total singles Fock term, RIJCOSX called by the simple keyword RCSinglesFock and RIJK called by RIJKSinglesFock, see Basics. For canonical and LPNO methods, by default the program chooses the same approximation used in the SCF calculation. DLPNO2013 uses RIJCOSX by default, while in DLPNO, the singles Fock term is also evaluated using PNOs via SinglesFockUsePNOs, see Local correlation. This behavior can also be changed by keywords in the method block.
%method RIJKSinglesFock 1 # 0 false, 1 true
RCSinglesFock 0 # 0 false, 1 true
end
7.14.6. Use of the MDCI Module¶
The MDCI module is fairly easy to use. The flags for the “simple” input lines have been described in section Keyword Lines. The detailed listing of options is found below:
%mdci citype CISD # CI singles+doubles
QCISD # quadratic CI (singles+doubles)
CCSD # coupled-cluster singles+doubles
CEPA_1 # coupled-electron pair approximation ''1''
CEPA_2 #
CEPA_3 #
NCEPA_1 # our slightly modified versions of CEPA
NCEPA_2 # and CPF
NCEPA_3 #
NCPF_1 #
NCPF_2 #
NCPF_3 #
ACPF # averaged coupled-pair functional
ACPF_2 # Gdanitz modification of it
NACPF # our modification of it
AQCC # Szalay + Bartlett
SEOI # a strictly size extensive energy functional
# maintaining unitary invariance (not yet
# available for UHF)
ewin -3,1e3 # orbital energy window to determine which
# MOs are included in the treatment
# (respects settings in %method block)
Singles true # include single excitations in the
# treatment (default true)
Triples 0 # (T) correction in CCSD(T)/QCISD(T)
# default is no triples
1 # Algorithm 1 (lots of memory, fast)
2 # Algorithm 2 (less memory, about 2x slower,
# not yet available for UHF)
Brueckner true # use Brueckner orbitals
# (default false)
Denmat none # no evaluation of density matrices
linearized # density matrix obtained by retaining
# only CEPA_0-like terms, i.e., those
# linear in the excitation amplitudes
unrelaxed # unrelaxed density matrices, i.e.,
# density matrices without orbital
# relaxation
orbopt # perform orbital optimization yielding
# fully relaxed density matrices (if
# citype chosen as CCSD or QCISD this option
# implies evaluation of the Z vector).
# (default: linearized)
ZSimple true # simplified evaluation of the Z vector
# in case of orbital optimized CCD
# (citype chosen as CCSD or QCISD and
# Denmat as orbopt) by using an
# analytical formula
false # explicit solution of Z vector
# equations
# in case of orbital optimized CCD
# (default: false)
UseQROs # use of quasi-restricted orbitals
# (default false)
Localize 0 # use localized MOs. Presently very little
# use is made of locality. It may help
# for interpretations. Localization is
# incompatible with the (T) correction
PM # Use Pipek-Mezey localized MOs
FB # use Foster-Boys localized MOs
NatOrbIters 0 # Perform natural orbital iterations.
# default is none. Not possible for CCSD
# and QCISD
pCCSDAB # the three parameters for parametrized
pCCSDCD # coupled-cluster (default is 1.0 which
pCCSDEF # corresponds to normal CCSD
# this defines how the rate limiting step is handled
# MO and AOX need lots of disk and I/O but if they
# can be done they are fast
KCOpt KC_MO # Perform full 4-index transformation
KC_AOBLAS# AO direct with BLAS (preferred)
# (not yet available for UHF, switch to KC_AOX)
KC_AO # AO direct handling of 3,4 externals
# (not yet available for UHF, switch to KC_AOX)
KC_RI # RI approximation of 3,4 externals
# (not yet available for UHF)
KC_RI2 # Alternative RI (not recommended)
# (not yet available for UHF)
KC_AOX # Do it from stored AO exchange integrals
PrintLevel 2 # Control the amount of output. For 3 and
# higher things like pair correlation
# energies are printed.
MaxIter 35 # Max. number of iterations
# How the integral transformation is done.
# Note that it is fine to do AOX or AO or AOBLAS
# together with trafo_ri
TrafoType trafo_jk # Partial trafo to J+K operators
trafo_ri # RI transformation of all
# integrals up to 2-externals
# (3-ext for (T))and rest on the
# fly
trafo_full # Full four index transformation.
# Automatically chosen for
# KCOpt=KC_MO
MaxCore 350 # Memory in MB - used for integral
# trafos and batching and for storage of
# integrals and amplitudes
# don't be too generous
STol 1e-5 # Max. element of the residual vector
# for convergence check
LShift 0.3 # Level shift to be used in update of
# coefficients
MaxDIIS 7 # Max number of DIIS vectors to be stored
# this lets you control how much and what is residing
# in central memory. May speed up things. Note that
# MaxCore is not respected here
InCore 0 # nothing in core
1 # + sigma-vector and amplitudes (default)
2 # + Jij(a,b) Kij(a,b) operators
3 # + DIIS vectors
4 # + 3-exernal integral Kia(b,c)
5 # + 4-external integrals Kab(c,d)
# this is identical to ALL
# the default is AUTO which means that incore
# is chosen based on MaxCore
end