7.51. Calculation of Properties¶
7.51.1. Electric Properties¶
Calculation of first order (electric dipole and quadrupole
moments) and second order (polarizabilities) electric properties can be requested
in the %elprop input block. The second order properties can be
calculated through the solution of the CP-SCF (see CP-SCF Options)
or CP-CASSCF (see CASSCF Linear Response) equations. Details are
shown below:
%elprop
Dipole true
Quadrupole true
Polar true
PolarVelocity true # polarizability w.r.t. velocity perturbations
PolarDipQuad true # dipole-quadrupole polarizability
PolarQuadQuad true # quadrupole-quadrupole polarizability
freq_r 0.00 # purely real frequency (default: static calculation)
freq_i 0.00 # purely imaginary frequency (default: static calculation)
Solver CG # CG(conjugate gradient)
# other options: DIIS or POPLE(default)
MaxDIIS 5 # max. dimension of DIIS method
Shift 0.2 # level shift used in DIIS solver
Tol 1e-3 # Convergence of the CP-SCF equations
# (norm of the residual)
MaxIter 64 # max. number of iterations in CPSCF
PrintLevel 2
Origin CenterOfElCharge # center of electronic charge
CenterOfNucCharge # center of nuclear charge
CenterOfSpinDens # center of spin density
CenterOfMass # center of mass (default)
N # position of atom N (starting at 0)
X,Y,Z # explicit position of the origin
# in coordinate input units (Angstrom by default)
end
The most efficient and accurate way to calculate the polarizability analytically is to use the coupled-perturbed SCF method. The most time consuming and least accurate way is the numerical second derivative of the total energy. Note that the numerical differentiation requires: (a) tightly or even very tightly converged SCF calculations and (b) carefully chosen field increments. If the field increment is too large then the truncation error will be large and the values will be unreliable. On the other hand, if the field increment is too small the numerical error associated with the finite difference differentiation will get unacceptably large up to the point where the whole calculation becomes useless.
7.51.2. The Spin-Orbit Coupling Operator¶
Several variants of spin-orbit coupling operators can be used for
property calculations [610]. They are based on effective
potential and mean-field approaches, and have various parameters that
can be selected via the %rel block. Note that the SOMF operator
depends on the density matrix, so the operator itself can differ for
example between a CASSCF and an MRCI calculation.
Note: The defaults have slightly changed in ORCA 5.0, see SOCFlags
in the following.
%rel
# ---------------------------------------------------
# SPIN ORBIT COUPLING OPERATORS
# ---------------------------------------------------
SOCType 0 # none
1 # effective nuclear charge
2 # mean-field with atomic densities read from
# disk; similar to SOCType=4
3 # mean-field/effective potential (default)
4 # mean-field with atomic densities generated
# on the fly; see bellow
# ---------------------------------------------------
# Flags for construction of potential; operative
# only for SOCType 3.
# ---------------------------------------------------
SOCFlags 1,4,3,0 # default if nothing is specified
# 1,3,3,0 default if RI is applied and AuxJ available
# e.g. when using !RIJCOSX (default for DFT) or !RIJONX
# Flag 1 = 0 - do not include 1-electron terms
# = 1 - do include 1-electron terms
# Flag 2 = 0 - do not include Coulomb terms
# = 1 - compute Coulomb terms fully numeric
# = 2 - compute Coulomb term seminumeric
# = 3 - compute Coulomb term with RI approx
# = 4 - compute Coulomb term exactly
# Flag 3 = 0 - do not include exchange terms
# = 1 - do include local X-alpha exchange
# the X-Alpha parameter can be chosen via
# % rel Xalpha 0.7 (default)
# = 2 - same as 1 but with sign reversed
# = 3 - exchange via one-center exact
# integrals including the spin-other
# orbit interaction
# = 4 - all exchange terms full analytic
# (this is expensive)
# Flag 4 = 0 - do not include DFT local correlation
# terms
# = 1 - do include local DFT correlation (here
# this is done with VWN5)
#
SOCMaxCenter 4 # max. number of centers to include in
# the integrals (not fully consistently
# implemented yet; better leave equal to 4)
# The following simple input equivalents can also be used:
# SOMF(1X) = SOCType 3, SOCFlags 1,2,3,0 and SOCMaxCenter 4
# RI-SOMF(1X) = SOCType 3, SOCFlags 1,3,3,0 and SOCMaxCenter 4
# VEFF-SOC = SOCType 3, SOCFlags 1,3,1,1 and SOCMaxCenter 4
# VEFF(-2X)-SOC = SOCType 3, SOCFlags 1,3,2,1 and SOCMaxCenter 4
# AMFI = SOCType 3, SOCFlags 1,4,3,0 and SOCMaxCenter 1
# AMFI-A = SOCType 4, SOCFlags 1,4,3,0 and SOCMaxCenter 1
# (AMFI-like approach that uses pre-calculated atomic densities;
# this can be combined with the SOCOff option
# to exclude contributions from specific atoms)
# NOTE: If you choose the RI option you need to specify an auxiliary basis set
# even if the underlying SCF calculation does not make use of any form
# of RI!
# -----------------------------------------------
# For the effective nuclear charge SOC operator
# the nuclear charges can be adjusted.
# -----------------------------------------------
Zeff[26] 0.0 # set the effective nuclear charge
# of iron (Z = 26) to zero
# -----------------------------------------------
# Neglecting SOC contributions from particular
# atoms
# -----------------------------------------------
SOCOff 0,5 # turn off the SOC for atoms 0 and 5
# this makes sense if the SOC operator
# has only one center contributions
# (e.g. effective nuclear charge)
Simple input equivalents are described in more detail in [610]. More details on the AMFI-A approach which uses pre-calculated atomic densities can be found in [284].
The Breit-Pauli spin-orbit coupling operator is given by:
with the one- and two-electron contributions
This operator would be hard to handle exactly; therefore it is common to introduce mean field and/or effective potential approaches in which the operator is written as an effective one-electron operator:
The simplest approximation is to simply use the one-electron part and regard the nuclear charges as adjustable parameters. Reducing their values from the exact nuclear charge is supposed to account in an average way for the screening of the nuclear charge by the electrons. In our code we use the effective nuclear charges of Koseki et al. This approximation introduces errors which are usually smaller than 10% but sometimes are larger and may approach 20% in some cases. The approximation is best for first row main group elements and the first transition row (2p and 3d elements). For heavier elements it becomes unreliable.
A much better approximation is to take the two-electron terms into account precisely. Without going into details here – the situation is as in Hartree-Fock (or density functional) theory and one gets Coulomb, exchange and correlation terms. The correlation terms (evaluated in a local DFT fashion) are negligible and can be safely neglected. They are optionally included and are not expensive computationally. The Coulomb terms is (after the one-electron term) the second largest contribution and is expensive to evaluate exactly. The situation is such that in the Coulomb-part the spin-other orbit interaction (the second term in the two-electron part) does not contribute and one only has to deal with the spin-own-orbit contribution. The exact evaluation is usually too expensive to evaluate. The RI and seminumeric approximation are much more efficient and introduce only minimal errors (on the order of usually not more than 1 ppm in g-tensor calculations for example) and are therefore recommended. The RI approximation is computationally more efficient. Please note that you have to specify an auxiliary basis set to take advantage of the RI approximation, even if the preceding SCF calculation does not make use of any form of RI. The one-center approximation to the Coulomb term introduces much larger errors. The fully numeric method is both slower and less accurate and is not recommended.
The exchange term has contributions from both the spin-own-orbit and
spin-other-orbit interaction. These are taken both into account in the
mean-field approximation which is accessed by Flag 3 = 3. Here a
one-center approximation is much better than for the Coulomb term since
both the integrals and the density matrix elements are short ranged.
Together with the Coulomb term this gives a very accurate SOC operator
which is recommended. The DFT-Veff operator suffers from not treating
the spin-other-orbit part in the exchange which gives significant errors
(also, local DFT underestimates the exchange contributions from the
spin-same-orbit interaction by some 10% relative to HF but this is not a
major source of error). However, it is interesting to observe that in
the precise analytical evaluation of the SOMF operator, the
spin-other-orbit interaction is exactly -2 times the spin-own-orbit
interation. Thus, in the DFT framework one gets a much better SOC
operator if the sign of the DFT exchange term is simply reversed! This is accessed by Flag 3 = 2.
7.51.2.1. Exclusion of Atomic Centers¶
In ORCA it is possible to change the spin-orbit coupling operator in order to exclude contributions from user-defined atoms. This approach can be useful, for example, in quantifying the contribution of the ligands to the zero-field splitting (ZFS); for an application of this method see Ref. [834].
This is illustrated for the calculation of the SOC contribution to the ZFS of the triplet oxygen molecule. Using the input below we start by a normal calculation of the ZFS, including both oxygen atoms. Note that we use here the effective nuclear charge operator. This is required as not all implemented SOC operators are compatible with the decomposition in terms of individual centers contributions.
! def2-TZVP def2-TZVP/c
%casscf nel 8
norb 6
mult 3,1
nroots 1,3
rel dosoc true
end
end
%rel
SOCType 1
end
*xyz 0 3
O 0 0 0
O 0 0 1.207
*
The calculated value of the D parameter is approximately 2.573
cm\(^{-1}\). In a second calculation we exclude the contribution from the
first oxygen atom. For this we change the %rel block to the one below.
%rel
SOCType 1
SOCOff 0
end
Now the D parameter is calculated to be approximately 0.643 cm\(^{-1}\), a result that deviates quite significantly from half of the value calculated previously, implying that non-additive effects are important. In addition to the effective nuclear charge operator, the AMFI-A operator described previously can be used. Given that this is based on precalculated atomic densities, it might be preferred for heavier elements where the effective nuclear charge operator becomes unreliable. The method is not limited to CASSCF calculations as described above, and can be used in DFT, MRCI and ROCIS calculations.
7.51.3. EPR and NMR properties¶
Calculation of EPR and NMR response properties can be requested in
the %eprnmr input block. The individual flags are given below.
%eprnmr
# Calculate the g-tensor using CP-KS theory
gtensor true
# Calculate and print one- and two-electron contributions to the g-tensor
gtensor_1el2el true
# Calculate the D-tensor
DTensor so # spin orbit part
ss # spin-spin part
ssandso # both parts
DSOC qro # quasi-restricted method; must be done with the keyword !uno
pk # Pederson-Khanna method.
# NOTE: both approximations are only valid for
# pure functionals (no HF exchange)
cp # coupled-perturbed method (default)
cvw # van W\"ullen method
DSS direct # directly use the canonical orbitals for the spin density
uno # use spin density from UNOs
PrintLevel n # Amount of output (default 2)
# whether to calculate and print the Euler angles via `orca_euler` if the
# calculation of the g-tensor or the D-tensor is requested
PrintEuler false
# For the solution of the CP-SCF equations:
Solver Pople # Pople solver (default)
CG # conjugate gradient
DIIS # DIIS type solver
MaxIter 64 # maximum number of iterations
MaxDIIS 10 # max. number of DIIS vectors (only DIIS)
Tol 1e-3 # convergence tolerance
LevelShift 0.05 # level shift for DIIS (ignored otherwise)
Ori CenterOfElCharge # center of electronic charge
CenterOfNucCharge # center of nuclear charge
CenterOfSpinDens # center of spin density
CenterOfMass # center of mass
GIAO # use the GIAO formalism (default)
N # number of the atom to put the origin
X,Y,Z # explicit position of the origin
# in coordinate input units (Angstrom by default)
# Calculate the NMR shielding tensor
NMRShielding 1 # for chosen nuclei - specified with the Nuclei keyword
2 # for all nuclei - equivalent to the 'NMR' simple input keyword
# treatment of 1-electron integrals in the RHS of the CPSCF equations
giao_1el = giao_1el_analytic # analytical, default
giao_1el_numeric # numerical - for testing only
# treatment of 2-electron integrals in the RHS of the CPSCF equations
# various options for combination of approximations in Coulomb (J) and
# exact (HF) exchange (K) part. The default is the same as used for the SCF.
giao_2el = giao_2el_rijk # RIJK approximation
giao_2el_same_as_scf # use same scheme as in SCF
giao_2el_analytic # fully analytical, for testing only
giao_2el_rijonx # RIJ approximation with analytical K
giao_2el_cosjonx # COSJ approximation with analytical K
giao_2el_rijcosx # RIJ approximation with COSX approximation
giao_2el_cosjx # COSJ approximation with COSX apprxomation
giao_2el_exactjcosx # exact J with COSX approximation
giao_2el_exactjrik # exact J with RIK approximation
# for g-tensor calculations using the SOMF-operator for the SOC
# treatment, the 2-electron contribution to the GIAO terms can be
# computed as well, but they take much more time and usually do not
# contribute significantly and therefore are disabled by default
do_giao_soc2el false
# treatment of tau in meta-GGA DFT - see below
Tau = Dobson # (default) Other options: 0, MS, GI
# use effective nuclear charges for the gauge correction to the A-tensor
# (this makes sense if an effective 1-electron SOC operator is used)
hfcgaugecorrection_zeff true
# calculate diamagnetic spin-orbit (DSO) integrals needed for the gauge
# correction to the A-tensor numerically (faster than the analytical way)
hfcgaugecorrection_numeric true
# Grid settings for the above: <0 means to use the DFT grid setting
hfcgaugecorrection_angulargrid -1
hfcgaugecorrection_intacc -1
hfcgaugecorrection_prunegrid -1
hfcgaugecorrection_bfcutoff -1
hfcgaugecorrection_wcutoff -1
Nuclei = all type { flags }
# Calculate nuclear properties. Here the properties
# for all nuclei of "type" are calculated ("type" is
# an element name, for example Cu.)
# Flags can be the following:
# aiso - calculate the isotropic part of the HFC
# adip - calculate the dipolar part of the HFC
# aorb - 2nd order contribution to the HFC from SOC
# fgrad - calculate the electric field gradient
# rho - calculate the electron density at the
# nucleus
# shift - NMR shielding tensor (orbital contribution)
# srot - spin-rotation tensor
# ssdso - spin-spin coupling constants, diamagnetic spin orbit term
# sspso - spin-spin coupling constants, paramagnetic spin orbit term
# ssfc - spin-spin coupling constants, Fermi contact term
# sssd - spin-spin coupling constants, spin dipole term
# ssall - spin-spin coupling constants, calculate all above contributions
# In addition you can change several parameters
# e.g. for a different isotope.
Nuclei = all N { PPP=39.1, QQQ=0.5, III=1.0 };
# PPP - the HFC proportionality factor for this nucleus
# = ge*gN*betaE*betaN
# QQQ - the quadrupole moment for this nucleus
# III - the spin for this nucleus
# ist - isotope
# ssgn - nuclear g-factor for spin-spin coupling
# and spin-rotation constants (overrides ist)
# For example:
# calculates the hyperfine coupling for all nitrogen atoms
Nuclei = all N { aiso, adip, fgrad, rho};
# calculates the spin-spin coupling constants for all carbon atoms
# assuming all carbon atoms are 13C
Nuclei = all C { ssall, ist = 13};
# You can also calculate properties for individual atoms
# as in the following example:
Nuclei = 1,5,8 { aiso, adip};
# NOTE: Counting starts with atom 1!
# WARNING: All the nuclei, mentioned in one line
# as above will be assigned the same parameters !
# For spin-spin coupling constants, a distance threshold is
# applied in the eprnmr block to restrict the number of couplings
# to be computed, given in Angstroms:
SpinSpinRThresh 5.0 # default
# Coupling can be restricted to certain element pairs
# (if they are added to the "Nuclei" list and are within RThresh).
# The syntax accepts multiple pairs of element symbols or
# atomic numbers or "*" as a wildcard
SpinSpinElemPairs {C C} {H *} {6 7} # default is {* *}, i.e. all
# Similarly, coupling can be restricted to certain atom pairs
# (if they are added to the "Nuclei" list and are within RThresh).
# The syntax accepts multiple pairs of indices (starting at 0!)
# or "*" as a wildcard
SpinSpinAtomPairs {1 0} {5 *} # default is {* *}, i.e. all
# whether to print reduced spin-spin coupling constants
PrintReducedCoupling false
end
7.51.3.1. Hyperfine and Quadrupole Couplings¶
The hyperfine coupling has four contributions:
(a) The isotropic Fermi contact term that arises from the finite spin density on the nucleus under investigation. It is calculated for nucleus \(N\) from:
Here, \(\left\langle { S_{z} } \right\rangle\) is the expectation value of the z-component of the total spin, \(g_{e}\) and \(g_{N}\) are the electron and nuclear g-factors and \(\beta_{e}\) and \(\beta_{N}\) are the electron and nuclear magnetons respectively. \(\rho \left({ \vec{{R} }_{N} } \right)\) is the spin density at the nucleus. The proportionality factor \(P_{N} = g_{e} g_{N} \beta_{e} \beta_{N}\) is commonly used and has the dimensions MHz bohr\(^{3}\) in ORCA.
(b) The spin dipole part that arises from the magnetic dipole interaction of the magnetic nucleus with the magnetic moment of the electron. It is also calculated as an expectation value over the spin density as:
where \(\mathrm{\mathbf{\rho } }\) is the spin-density matrix and \(\vec{{r} }_{N}\) is a vector of magnitude \(r_{N}\) that points from the nucleus in question to the electron (\(\left\{\phi \right\}\) is the set of basis functions).
(c) The second order contribution that arises from spin-orbit coupling. Presently ORCA can calculate all these contributions. The first two are calculated as simple expectation values of the appropriate operators over the self-consistent spin density, but the second order contribution requires the solution of the coupled-perturbed SCF equations and is consequently computationally more demanding. The contribution can be written:
The derivative of the spin density is computed from solving the coupled-perturbed SCF equations with respect to the nucleus-orbit coupling as perturbation. The nucleus-orbit coupling is represented by the operator
where the sum is over electrons and \(A\) is the nucleus in question.
(d) The gauge-correction contribution.[492] This term is often small. However, it is needed in order to get exactly gauge-invariant results. We recently showed that the gauge correction can become crucial in the long-distance limit between the nuclear spin and the electron spin. This is relevant for pseudocontact NMR chemical shifts (PCS).[492]
The field gradient tensor is closely related to the dipole contribution to the hyperfine coupling. The main differences are that the electron instead of the spin density enters its calculation and that it contains a nuclear contribution due to the surrounding nuclei. It is calculated from
with \(Z_{A}\) as the nuclear charge of nucleus \(A\) and \(\vec{{R} }_{AN}\) as a vector of magnitude \(R_{AN}\) that points from nucleus \(A\) to nucleus \(N\). \(\mathrm{\mathbf{P} }\) is the first order density matrix.
NOTE:
Hyperfine and quadrupole couplings are properties where the standard basis sets that have been designed for geometry optimization and the like may not be entirely satisfactory (especially for atoms heavier than Ne). You should probably look into tailoring the basis set according to your needs. While it is likely that a later release will provide one or two special basis sets for “core-property” calculations at this time you have to make sure yourself that the basis set has enough flexibility in the core region, for example by uncontracting core basis functions and adding s-primitives with large exponents (or using the “decontraction feature”, section Choice of Basis Set). If you add these tight functions and use DFT make sure that the numerical integration is still satisfactory. Use the “SpecialGrid” feature to enlarge grids for individual atoms without increasing the computational effort too drastically.
For heavy nuclei you may want to consider the possibility of relativistic effects. Scalar relativistic effects can be handled with several quasi-relativistic Hamiltonians in ORCA, an overview of the possibilities and some recommendations can be found in Relativistic Calculations. Note that relativistic calculations may have special requirements on basis sets, and in the context of property calculations, you should be especially aware of the importance of using picture change corrections (see Relativistic Calculations and Picture-Change Effects). In quasi-relativistic calculations with DFT, one should also be very cautious about accuracy of the numerical integration, especially for heavier (transition metal) nuclei.
Second order HFCs require the calculation of the spin-orbit coupling contributions which in turn requires the calculation of the coupled perturbed SCF equations. These effects can be quite significant for heavier nuclei and should definitely be included for transition metal complexes. The spin-orbit coupling treatment used is the same as described under The Spin-Orbit Coupling Operator.
7.51.3.2. The g-Tensor¶
The EPR g-tensor is a property that can be calculated as a second derivative of the energy and it is implemented as such in ORCA for the SCF methods—e.g. HF and DFT—, CASSCF, as well as all-electron MP2 (or RI-MP2) and double-hybrid DFT. The following four contributions arise for the g-tensor (SZ = spin Zeeman, RMC = relativistic mass correction, DSO = diamagnetic spin-orbit correction, PSO = orbital Zeeman/SOC term):
Here, \(g_{e}\) is the free-electron g-value (\(=\)2.002319…), \(S\) is the total spin, \(\alpha\) the fine structure constant, \(P^{\alpha -\beta }\) is the spin density matrix, \(\phi\)is the basis set, \(\hat{T}\) is the kinetic energy operator, \(\xi \left({ r_{A} } \right)\) an approximate radial operator, \(h^{\text{SOC} }\) the spatial part of an effective one-electron spin-orbit operator and \(B_{\mu }\) is a component of the magnetic field. The calculation of the derivative of the spin-density depends on the chosen level of theory. For the SCF-level it is done based on the coupled-perturbed SCF theory with respect to a magnetic field perturbation.
Accuracy. g-tensor calculations at the SCF level are not highly demanding in terms of basis set size. Basis sets that give reliable SCF results (at least valence double-zeta plus polarization) usually also give reliable g-tensor results. For many molecules the Hartree-Fock approximation will give reasonable predictions. In a number of cases, however, it breaks down completely. DFT is more robust in this respect and the number of molecules where it fails is much smaller. Among the density functionals, the hybrid functionals seem to be the most accurate. In my hands PBE0 is perhaps the best although PWP1 and B3LYP are not much worse. The GGA functionals such as BP, PW91, BLYP or PBE are equally good for small radicals but are significantly inferior to their hybrid counterparts for transition metal complexes.
Gauge dependence. Unfortunately, the g-tensor is a gauge dependent property, i.e. the results depend on where the origin is chosen within the molecule. Unless fully invariant procedures (such as GIAOs) are used, this undesirable aspect is always present in the calculations. GIAOs are now available for calculations on the SCF-level in ORCA. However, if the choice of gauge origin is not outrageously poor, the gauge dependence is usually so small that it can be ignored for all practical purposes, especially if large basis sets are used. ORCA gives you considerable freedom in the choice of gauge origin. It can either be the center of mass, the center of nuclear charge, the center of electronic charge, GIAOs (recommended if available), a special atom or a user-defined point in space. It is wise to check the sensitivity of the results with respect to the choice of origin, especially when small g-shifts on the order of only a few hundred ppm are calculated.
Spin-orbit coupling operator. In previous versions of the code, the g-tensor module used the parameterization of Koseki et al. [462, 463, 464] for the spin-orbit operator. This is expected to be a reasonable approximation for the 2p and 3d elements and less satisfactory for heavier main group or transition metal containing systems. Thus, the main target molecules with the simple operators are radicals made of light atoms and first row transition metal complexes. More accurate SOC operators (at only moderately increased computational cost) have now been implemented and are described in section The Spin-Orbit Coupling Operator. With these operators there are fewer restrictions. However, for very heavy elements they will suffer from the shortcomings of the Breit-Pauli approximation and future releases will modify these operators to take into account the ZORA or DKH corrections to the SOC.
7.51.3.3. Zero-Field-Splitting¶
It is well known that the ZFS consists of a first order term arising from the direct spin-spin interaction[364]:
(\(K\),\(L=\)x,y,z). Here \(\alpha\) is the fine structure constant (\(\approx\) 1/137 in atomic units), \(\mathrm{\mathbf{r} }_{ij}\) is the electronic distance vector with magnitude \(r_{ij}\) and \(\hat{{s} }_{i}\) is the spin-vector operator for the \(i\)’th electron. \(\left|0SS \right\rangle\) is the exact ground state eigenfunction of the Born-Oppenheimer Hamiltonian with total spin \(S\) and projection quantum number \(M_{S} = S\). Since the spin-spin interaction is of first order, it presents no particular difficulties. The more complicated contribution to the D-tensor arises from the spin-orbit interaction, which gives a second order contribution. Under the assumption that the spin-orbit coupling (SOC) operator can to a good approximation be represented by an effective one-electron operator (\(\hat{{H} }_{\text{SOC} } =\sum\nolimits_i { \mathrm{\mathbf{\hat{{h} }} }_{i}^{\text{SOC} } { \mathrm{\mathbf{\hat{s} }} }_{i} } )\), ref [621] has derived the following sum-over-states (SOS) equations for the SOC contribution to the ZFS tensor:
Here the one-electron spin-operator for electron \(i\) has been written in terms of spherical vector operator components \(s_{i,m}\) with \(m=0,\pm 1\) and \(\Delta_{b} =E_{b} -E_{0}\) is the excitation energy to the excited state multiplet \(\left| b^{SS} \right\rangle\) (all \(M_{S}\) components are degenerate at the level of the BO Hamiltonian).
One attractive possibility is to represent the SOC by the spin-orbit mean-field (SOMF) method developed by Hess et al.,[390] widely used in the AMFI program by Schimmelpfennig [761] and discussed in detail by Berning et al.[97] as well as in ref. [610]. In terms of an (orthonormal) one-electron basis, the matrix elements of the SOMF operator are:
and:
\(\mathrm{\mathbf{\hat{l} }}_{iA} =(\mathrm{\mathbf{\hat{r} }}_{i} -\mathrm{\mathbf{R} }_{A} )\times \mathrm{\mathbf{\hat{p} }}_{i}\) is the angular momentum of the \(i\)’th electron relative to nucleus \(A\). The vector \(\mathrm{\mathbf{\hat{r} }}_{iA} ={\mathrm{\mathbf{\hat{r} }} }_{i} -\mathrm{\mathbf{R} }_{A}\) of magnitude \(r_{iA}\) is the position of the \(i\)’th electron relative to atom \(A\). Likewise, the vector \({\mathrm{\mathbf{\hat{r} }} }_{ij} =\mathrm{\mathbf{\hat{r} }}_{i} -\mathrm{\mathbf{\hat{r} }}_{j}\) of magnitude \(r_{ij}\) is the position of the \(i\)th electron relative to electron \(j\) and \(\mathrm{\mathbf{\hat{l} }}_{ij} =(\mathrm{\mathbf{\hat{r} }}_{i} -\mathrm{\mathbf{\hat{r} }}_{j} )\times \mathrm{\mathbf{\hat{p} }}_{i}\) is its angular momentum relative to this electron. P is the charge density matrix of the electron ground state (\(P_{pq} =\left\langle { 0SS\left|{E_{q}^{p} } \right|0SS} \right\rangle\)with \(E_{q}^{p} =a_{p\beta }^{+} a_{q\beta } +a_{p\alpha }^{+} a_{q\alpha }\) where \(a_{p\sigma }^{+}\) and \(a_{q\sigma }\) are the usual Fermion creation and annihilation operators).
7.51.3.4. General Treatment of ZFS¶
The zero-field splitting (ZFS) is typically the leading term in the Spin-Hamiltonian (SH) for transition metal complexes with a total ground state spin \(S\)>1/2 (for reviews and references see chapter Publications Related to ORCA). Its net effect is to introduce a splitting of the \(2S+1\) \(M_{S}\) levels (which are exactly degenerate at the level of the Born-Oppenheimer Hamiltonian), even in the absence of an external magnetic field. Thus, an analysis and interpretation of the ZFS is imperative if the information content of the various physical methods that are sensitive to ZFS effects.
In 2007, we have developed a procedure that makes the ZFS calculation compatible with the language of analytic derivatives.[613] Perhaps the most transparent route is to start from the exact solutions of the Born-Oppenheimer Hamiltonian. To this end, we look at the second derivative of the ground state energy (\(E=\left\langle { 0^{SS}\left|{\hat{{H} }} \right|0^{SS} } \right\rangle)\) with respect to a spin-dependent one-electron operator of the general form:
Where \(h_{pq}^{K}\) is the matrix of the \(K\)’th component of the spatial part of the operator (assumed to be imaginary Hermitian as is the case for the spatial components of the SOC operator) and \(\hat{{S} }_{pq}^{\left( m \right)}\) is the second quantized form of the spin vector operator (\(m=0,\pm 1\)). The quantity \(x_{K}^{\left( m \right)}\) is a formal perturbation parameter. Using the exact eigenfunctions of the BO operator, the first derivative is:
With the components of the spin density:
The second derivative with respect to a second component for \(m'=-m\) is:
The derivative of the spin density may be written as:
Expanding the perturbed wavefunction in terms of the unperturbed states gives to first order:
Where \(\left| n \right\rangle\) is any of the \(\left| b^{S'M'} \right\rangle\). Thus, one gets:
The equality holds for exact states. For approximate electronic structure treatments, the analytic derivative approach is more attractive since an infinite sum over states can never be performed in practice and the calculation of analytic derivative is computationally less demanding than the calculation of excited many electron states.
Using eq. (7.372), the components of the SOC-contribution to the D-tensor are reformulated as
These are general equations that can be applied together with any non-relativistic or scalar relativistic electronic structure method that can be cast in second quantized form. Below, the formalism is applied to the case of a self-consistent field (HF, DFT) reference state.
For DFT or HF ground states, the equations are further developed as follows:
The SCF energy is:
Here \(V_{\text{NN} }\) is the nuclear repulsion energy and \(h_{\mu \nu }\) is a matrix element of the one-electron operator which contains the kinetic energy and electron-nuclear attraction terms (\(\left\langle { \mathrm{\mathbf{ab} }} \right\rangle\) denotes the trace of the matrix product \(\mathrm{\mathbf{ab} }\)). As usual, the molecular spin-orbitals \(\psi_{p}^{\sigma }\) are expanded in atom centered basis functions (\(\sigma =\alpha ,\beta )\):
with MO coefficients \(c_{\mu p}^{\sigma }\). The two-electron integrals are defined as:
The mixing parameter \(a_{\text{X} }\) controls the fraction of Hartree-Fock exchange and is of a semi-empirical nature. \(E_{\text{XC} } \left[{ \rho_{\alpha } ,\rho _{\beta } } \right]\) represent the exchange-correlation energy. The parameter \(c_{\text{DF} }\) is an overall scaling factor that allows one to proceed from Hartree-Fock theory (\(a_{\text{X} } =1, c_{\text{DF} } =0)\) to pure DFT (\(a_{\text{X} } =0, c_{\text{DF} } =1)\) to hybrid DFT (\(0<a_{\text{X} } <1, c_{\text{DF} } =1\)). The orbitals satisfy the spin-unrestricted SCF equations:
With \(V_{XC}^{\sigma } =\frac{\delta E_{XC} }{\delta \rho_{\sigma } \left({\mathrm{\mathbf{r} }} \right)}\) and \(P_{\mu \nu } =P_{\mu \nu }^{\alpha } +P_{\mu \nu }^{\beta }\) being the total electron density. For the SOC perturbation it is customary to regard the basis set as perturbation independent. In a spin-unrestricted treatment, the first derivative is:
For the second derivative, the perturbed orbitals are required. However, in the presence of a spin-dependent perturbation they can no longer be taken as pure spin-up or spin-down orbitals. With respect to the \(L\)’th component of the perturbation for spin-component \(m\), the orbitals are expanded as:
Since the matrix elements of the spin-vector operator components are purely real and the spatial part of the SOC operator has purely imaginary matrix elements, it follows that the first order coefficients are purely imaginary. The second derivative of the total SCF energy becomes:
Examination of the three cases \(m=0, \pm 1\) leads to the following equations for the D-tensor components:
Where a special form of the perturbed densities has been chosen. They are given in the atomic orbital basis as:
The special form of the coupled perturbed equations are implemented in ORCArun as follows: The perturbation is of the general form \(h^{K}\hat{{s} }_{m}\). The equations (7.385)-(7.390) and (7.391)-(7.396) below have been written in such a way that the spin integration has been performed but that the spin-dependent factors have been dropped from the right-hand sides and included in the prefactors of eqs. (7.385)-(7.387). The explicit forms of the linear equations to be solved are:
\(m=0\):
\(m=+1\):
\(m=-1\):
Note that these coupled-perturbed (CP) equations contain no contribution from the Coulomb potential or any other local potential such as the exchange-correlation potential in DFT. Hence, in the absence of HF exchange, the equations are trivially solved:
It is interesting that the “reverse spin flip coefficients” \(U_{a_{\beta } i_{\alpha } }^{K\left({ +1} \right)}\) and \(U_{a_{\alpha } i_{\beta } }^{K\left({ -1} \right)}\) are only nonzero in the presence of HF exchange. In a perturbation expansion of the CP equations they arise at second order (V\(^{2}\)/\(\Delta \varepsilon^{2})\) while the other coefficients are of first order (V/\(\Delta \varepsilon\); V represents the matrix elements of the perturbation). Hence, these contributions are of the order of \(\alpha^{4}\) and one could conceive dropping them from the treatment in order to stay consistently at the level of \(\alpha^{2}\). These terms were nevertheless kept in the present treatment.
Equations (7.391)-(7.396) are referred to as CP-SOC (coupled-perturbed spin-orbit coupling) equations. They can be solved by standard techniques and represent the desired analogue of the CP-SCF magnetic response equations solved for the determination of the g-tensor and discussed in detail earlier [663]. It is readily confirmed that in the absence of HF exchange, eqs. (7.397)-(7.402) inserted into eqs. (7.385)-(7.390) lead back to a modified Pederson-Khanna type treatment of the SOC contributions to the D-tensor [661]. In the framework of the formalism developed above, the Pederson-Khanna formula can be re-written in the form:
This equation was derived from second-order non-self-consistent perturbation theory without recourse to spin-coupling. For the special case of no Hartree-Fock exchange, the main difference to the treatment presented here is that the correct prefactors from eqs. (7.374)-(7.376) occur in front of the spin-flip contributions rather than \(\pm\) 1/(4\(S^{2})\) in eq. (7.403). In the presence of HF exchange it is suggested that the consistent generalization of the PK method are eqs. (7.385)-(7.387) with the \(\pm\) 1/(4\(S^{2})\) prefactors and this way the method has been implemented as an option into the ORCA program.
For completeness, the evaluation of the spin-spin term in the SCF case proceeds conveniently through:
as derived by McWeeny and Mizuno and discussed in some detail by Sinnecker and Neese.[799] In this reference it was found that DFT methods tend to overestimate the spin-spin contribution if the calculations are based on a spin-unrestricted SCF treatment. A much better correlation with experiment was found for open-shell spin restricted calculations. The origin of this effect proved to be difficult to understand but it was shown in ref [354] that in the case of small spin-contamination, the results of ROKS calculations and of those that are obtained on the basis of the spin-unrestricted natural orbital (UNO) determinant are virtually indistinguishable. It is therefore optionally possible in the ORCA program to calculate the spin-spin term on the basis of the UNO determinant.
7.51.3.5. Spin-rotation constants¶
Spin-rotation constant calculations are implemented using perturbation-dependent atomic orbitals following [290]. As given in eq. 34 of that reference, the spin-rotation tensor of nucleus \(K\), \(\mathbf{M}_K\), is related to the nuclear shielding tensor computed with GIAOs, \(\boldsymbol{\sigma}_K^\text{GIAO}\), and the diamagnetic part of the shielding tensor with the gauge origin set at that nucleus, \(\boldsymbol{\sigma}_K^\text{dia}(\mathbf{R}_K)\):
where \(\mathbf{M}_K^\text{nuc}\) is the nuclear component (given in eq. 14 of the reference), \(\mathbf{I}\) is the inertia tensor, and \(\gamma_K\) is the nuclear magnetogyric ratio. Accordingly, upon requesting spin-rotation constants, ORCA automatically computes the NMR shieldings with GIAOs as well. The following input shows an example calculation of \(\mathbf{M}(^{17}\text{O})\) in H212C17O:
! HF pcSseg-1 Mass2016 Bohrs
*xyz 0 1
O -0.00000000 -0.00000000 1.13863731 M=16.999131
C -0.00000000 -0.00000000 -1.14131773 M=12.0
H -0.00000000 1.76770755 -2.24076285
H 0.00000000 -1.76770755 -2.24076285
*
%eprnmr
Nuclei = all O {srot, ist=17}
end
Note
The magnetogyric ratio used can be changed either by choosing the correct isotope via
ist, or by providing the nuclear g-factor directly viassgn.The masses used to compute the inertia tensor are independent of the chosen isotopes! The example above requests atomic masses of the most abundant isotope (via the
Mass2016keyword) and explicitly specifies those of 12C (which is the default) and 17O.
7.51.3.6. Cartesian Index Conventions for EPR and NMR Tensors¶
The NMR shielding tensor \(\boldsymbol{\sigma}\) and the EPR \(g\) and \(A\) tensors are in general nonsymmetric matrices. It is therefore important to know the conventions used with regard to their cartesian indices. These conventions stipulate the order of the vector–matrix–vector multiplications in the respective spin Hamiltonians. Unless stated otherwise, ORCA adopts the following conventions:
For the NMR shielding tensor the nuclear Zeeman Hamiltonian assumes the form:
where \(\mathbf{B}\) is the applied magnetic field vector.
For the EPR \(g\) and \(A\) tensors the EPR spin Hamiltonian assumes the form:
7.51.3.7. MP2 level magnetic properties¶
Presently, hyperfine couplings (excluding the \(A_\text{orb}\) term),
g-tensors, and shielding tensors without GIAOs can be calculated for
both canonical and RI-MP2 and double-hybrid DFT without the frozen core
approximation. The \(A_\text{orb}\) term of the hyperfine couplings is
available only for RI-MP2 and double-hybrid DFT with and without frozen
core approximation. In case the RIJCOSX approximation is used, the
keywords Z_GridX, Z_GridX_RHS, KCOpt, KC_GridX and KC_IntAccX
are relevant – see
sections RIJCOSX-RI-MP2 Gradients and
MP2 and RI-MP2 Second Derivatives. NMR shielding and g-tensor
calculations with GIAOs are available for RI-MP2 and double-hybrid DFT
with or without a frozen core. The implementation is described in detail
in refs [828, 849] and the available
options are shown in
section RI-MP2 and Double-Hybrid DFT Response Properties. Note that for double-hybrid
DFT the correct properties are printed after the density heading containing
“Method : MP2” and “Level : Relaxed density”.
DLPNO-MP2 (and double-hybrid) shielding tensors are also available - see
section Local MP2 Response Properties.
7.51.3.8. Nucleus-independent chemical shielding¶
Aromaticity is a fundamental concept in chemistry and much attention has
been payed to its analysis in quantum chemistry. One possibility to gain
insight into aromaticity is to sample the ring current effect in NMR
close to the \(\pi\)-system. As this is not done by inspecting the
chemical shielding of any of the atoms, this quantity is called nucleus
independent chemical shielding (NICS). Usually, a “dummy” atom is placed
in the center of the ring and/or at some distance away from it. In ORCA,
one needs to use a ghost atom, e.g. “H:” to ensure that the program
generates DFT or COSX grid points in the region on interest. For
technical reasons, this atom must also have at least one basis function,
which can be set with NewGTO. An s-function with a sufficiently large
exponent will not overlap with any other basis function in the molecule
and will thus have no effect on the results, but only satisfy the
technical requirement (note that the extra grid points may change some
of the calculation results by increasing the accuracy of the numerical
integration). If RI is used, a dummy fitting function must also be added
to the AuxJ, AuxJK, and/or AuxC basis set. A typical input for benzene
looks like this:
! TightSCF NMR PBE def2-TZVP RI def2/J
* gzmt 0 1
H: NewGTO S 1 1 1e6 1 end NewAuxJGTO S 1 1 2e6 1 end
H: 1 1 NewGTO S 1 1 1e6 1 end NewAuxJGTO S 1 1 2e6 1 end
C 1 1.39 2 90
C 1 1.39 2 90 3 60
C 1 1.39 2 90 4 60
C 1 1.39 2 90 5 60
C 1 1.39 2 90 3 -60
C 1 1.39 2 90 6 60
H 3 1.09 4 120 1 180
H 4 1.09 3 120 1 180
H 5 1.09 4 120 1 180
H 6 1.09 5 120 1 180
H 7 1.09 3 120 1 180
H 8 1.09 7 120 1 180
*
7.51.3.9. Shielding tensor orbital decomposition¶
It is possible to decompose the NMR shielding tensor into orbital- or orbital pair contributions. One such option is the Natural Chemical Shielding analysis (see Section Natural Chemical Shielding Analysis (NCS)), while another is presented here. The shielding tensor for nucleus \(A\) can be exactly decomposed as follows:
Note that for SCF methods (HF or DFT),
\(\boldsymbol{\sigma}^{A,\text{para/dia} }_a=0\). To request the analysis,
a valid GBW file must be given with the keyword LocOrbGBW, which can
contain any orthonormal MOs in the same basis, e.g. canonical or
localized, although the decomposition above assumes that the Brillouin
condition is fulfilled and may be misleading if performed over NBOs for
example. Orbital contributions over 0.01 ppm are printed – this is
currently not user-configurable. The separate orbital pair contributions
can also be printed using Printlevel=3. The following example input
calculates HF and RI-MP2 shieldings for formaldehyde and decomposes them
over Pipek–Mezey localized orbitals (note that the virtual orbitals are
likely not well localized - AHFB would be better suited there). The
second sub-calculation just prints the LMOs for convenient visualization
with Avogadro.
! RI-MP2 NMR TightSCF RIJK pcSseg-1 cc-pwCVDZ/C def2/JK
%base "H2CO_NMR_PM"
%loc
LocMet PM # Pipek-Mezey localization
Occ true # localize both occupied
Virt true # and virtual orbitals (better use AHFB for these!)
T_core -1e5 # including the core
end
%eprnmr
PrintLevel 3 # print orbital pair contributions (lots of output!)
LocOrbGBW "H2CO_NMR_PM.loc" # use these orbitals
end
!xyzfile # store the coordinates
*gzmt 0 1
O
C 1 1.2078
H 2 1.1161 1 121.74
H 2 1.1161 1 121.74 3 180
*
# read the localized orbitals and print them in the output
# for easy visualization with Avogadro
$new_job
! HF pcSseg-1 NoRI MORead NoIter PrintBasis PrintMOs
%moinp "H2CO_NMR_PM.loc"
*xyzfile 0 1 H2CO_NMR_PM.xyz
7.51.3.10. Treatment of Tau in Meta-GGA Functionals¶
For GIAO-based calculations with meta-GGAs, different options are available for the kinetic energy density \(\tau\). The current-independent \(\tau_0\) is not gauge-invariant. Ignoring the terms, which produce the gauge-dependence, leads to an ad-hoc gauge-invariant treatment (this was the default up to ORCA 5). A gauge-invariant definition \(\tau_\text{MS}\), containing an explicit dependence on the magnetic field, was proposed by Maximoff and Scuseria.[566] However, this ansatz produces unphysical paramagnetic contributions to the shielding tensor.[757] The last option, introduced by Dobson,[216] is gauge-invariant but requires the solution of the CP-SCF equations, even for pure density functionals. For a discussion and comparison of these alternatives see refs [79] (in the context of TDDFT) and [713, 757] (in the context of NMR shielding). Note that the calculated shieldings can differ substantially between the different approaches! Some other electronic structure programs use the MS ansatz by default, so be careful when comparing results between different codes. In ORCA the treatment of \(\tau\) in GIAO-based calculations is chosed as follows
%eprnmr
Tau 0 # gauge-variant
GI # ad-hoc gauge-invariant
MS # field-dependent, gauge-invariant version of Maximoff and Scuseria
Dobson # (default) current density-dependent, gauge-invariant version of Dobson
end
7.51.4. Paramagnetic NMR shielding tensors¶
For systems with spin \(S>0\), the nuclear shielding contains a contribution which arises from the paramagnetism of the unpaired electrons.[1] This contribution is temperature-dependent and is called the “paramagnetic shielding” (\(\boldsymbol{\sigma}^\mathrm{p}\)). It adds to the temperature-indendent contribution to the shielding, also called the “orbital” contribution:
ORCA currently supports the calculation of \(\boldsymbol{\sigma}^\mathrm{p}\) for systems whose paramagnetism can be described by the effective EPR spin Hamiltonian
The theoretical background can be found in Refs. [809, 862]. We reproduce here the main equations.
For a spin state described by Eq. (7.406), the paramagnetic shielding tensor is given by
where \(\mathbf{Z}\) is a dimensionless \(3\times3\) matrix which is determined by the ZFS and the temperature, as follows: Diagonalization of the ZFS Hamiltonian \(\mathbf{S}\mathbf{D}\mathbf{S}\) yields energy levels \(E_\lambda\) and eigenstates \(|S\lambda a\rangle\), where \(a\) labels degenerate states if \(E_\lambda\) is degenerate. Then \(\mathbf{Z}\) is defined as (\(i,j=x,y,z\))
where \(Q_0 = \sum_{\lambda',a} e^{-E_\lambda/kT}\) denotes the partition function. An important property of the \(\mathbf{Z}\) matrix as defined above is that it goes to the identity matrix as \(\mathbf{D}/kT\) goes to zero.
The orbital part of the shielding, \(\boldsymbol{\sigma}^\mathrm{orb}\), is calculated in the same manner as for closed-shell molecules. It is available in ORCA for the unrestricted HF and DFT methods and for MP2 (see Section MP2 level magnetic properties for more information on the latter).
The orca_pnmr tool uses Eq.
(7.407) to calculate
\(\boldsymbol{\sigma}^\mathrm{p}\). Usage of orca_pnmr is described in
Section
orca_pnmr.
7.51.5. Calculating properties from existing densities¶
Occasionally, one may calculate a density matrix using an expensive correlated method such as CCSD and realize afterwards that a certain property such as the quadrupole moment or a hyperfine coupling constant (HFCC) is also required. Rather than start the whole calculation from scratch, one may wish to use the existing density matrix to calculate the properties. For this purpose, we have experimentally introduced a “properties only” calculation mode, whereby the MOs are read from an existing BaseName.gbw file and the densities are read from an existing BaseName.densities file and only the property calculations are performed. Note however, that this presents many possibilities for error, so only use it as a last resort and be very careful with the results!
Take, for example, this CCSD calculation:
! ccsd def2-svp
%base "BO-CCSD"
%mdci density unrelaxed end
*xyz 0 2
B 0 0 0
O 0 0 1.2049
*
This produces the files BO-CCSD.gbw and BO-CCSD.densities. To obtain
the CCSD quadrupole moment and HFCCs without repeating the whole
calculation, we can copy these two files into a new directory (highly
recommended!) and start a second job with the !PropertiesOnly keyword.
Note that the basename of the second job must be identical!
! ccsd def2-svp
%base "BO-CCSD"
%mdci density unrelaxed end
*xyz 0 2
B 0 0 0
O 0 0 1.2049
*
# Everything above must be the same as in the first job!
# Request the property calculations
! PropertiesOnly
%elprop
quadrupole true
end
%eprnmr
Nuclei = all B { Aiso, Adip }
Nuclei = all O { Aiso, Adip }
end
7.51.6. Local Energy Decomposition¶
The DLPNO-CCSD(T) method provides very accurate relative energies and allows to successfully predict many chemical phenomena. In order to facilitate the interpretation of coupled cluster results, we have developed the Local Energy Decomposition (LED) analysis [33, 105, 767], which permit to divide the total DLPNO-CCSD(T) energy (including the reference energy) into physically meaningful contributions. A practical guide to the LED scheme is reported in Section Local Energy Decomposition. Examples of applications of this scheme can be found in Ref. [84, 106, 294, 536, 537, 903]
As a word of caution we emphasize that only the total energy is an observable and its decomposition is, to some extent, arbitrary. Nevertheless, the LED analysis appears to be physically well grounded and logical to us, it is straightforward to apply and comes typically at a negligible computational cost compared to DLPNO-CCSD(T) calculations. Starting from ORCA 4.1, the LED scheme is available for both closed shell and open shell calculations. The code has also been made parallel and more efficient.
The LED scheme makes no assumption about the strength of the intermolecular interaction and hence it remains valid and consistent over the entire potential energy surface. Alternative schemes, such as the popular symmetry adapted perturbation theory, are perturbative in nature and hence are best applied to weakly interacting systems.
The idea of the LED analysis is rather simple. In local correlation methods occupied orbitals are localized and can be readily assigned to pre-defined fragments in the molecule. The same can be done for the correlation energy in terms of pair correlation energies that refer to pairs of occupied orbitals. In this way, both the correlation and the reference energy can be decomposed into intra- and interfragment contributions. The fragmentation is user defined. An arbitrary number of fragments can be defined. In the case that more than 2 fragments are defined, the interfragment interaction is printed for each fragment pair.
A very important feature of the LED scheme is the possibility to distinguish between dispersive and non-dispersive part of the DLPNO-CCSD(T) correlation energy. In brief, we exploit the fact that each CCSD pair correlation energy contribution can be expressed as a sum of double excitations from pairs of occupied orbitals into the virtual space. As in the DLPNO-CCSD(T) scheme the virtual space is spanned by Pair Natural Orbitals(PNOs) that are essentially local, the entire correlation energy can be decomposed into double excitations types, depending on where occupied and virtual orbitals are localized. For each pair of fragments, the sum of all excitations corresponding to the interaction of instantaneous local dipoles located on different fragments defines the so called “London dispersion” attraction between the two fragments in the LED framework.
For a system of two fragments, one can use as a reference point the geometrically and electronically relaxed fragments that constitute the interacting super-molecule. Relative to this reference state, the binding energy between the fragments can be written as:
where \(\Delta E_{geo-prep}\) is the energy needed to distort the fragments from their equilibrium configuration to the interacting geometry (also called “strain” in other energy decomposition schemes). The \(\Delta E_{el-prep}^{ref.}\) term represents the electronic preparation energy and describes how much energy is necessary to bring the fragments into the electronic structure that is optimal for interaction. \(E_{exch}^{ref.}\) is the inter-fragment exchange interaction (it always gives a binding contribution in our formalism) and \(E_{elstat}^{ref.}\) is the electrostatic energy interaction between the distorted electronic clouds of the fragments. The sum of these terms gives the Hatree-Fock energy in the closed shell case and the energy of the QRO determinant in the open shell case. Finally, the correlation energy is decomposed into dispersive \(E^{C-CCSD}_{dispersion}\) and non-dispersive \(\Delta E^{C-CCSD}_{non-dispersion}\) contributions plus a triples correction term to the interaction energy \(\Delta E^{C-(T) }_{int}\).
The \(E^{C-CCSD}_{dispersion}\) term contains the London dispersion contribution from the strong pairs described above plus the interfragment component of the weak pairs, which is essentially dispersive in nature. The \(\Delta E^{C-CCSD}_{non-dispersion}\) correlation term serves to correct the contributions to the binding energy approximately included in the reference energy, e.g it counteracts the overpolarization typical of the HF method. It contains the so called charge transfer excitations from the strong pairs \(E_{C-SP}^{CT(X\rightarrow Y) } + E_{C-SP}^{CT(X\leftarrow Y) }\), which represent the sum of all double excitation contributions that do not conserve the charge within each fragment. Moreover, the non-dispersive term also includes the electronic preparation from strong (\(\Delta E_{el-prep}^{C-SP}\)) and weak (\(\Delta E_{el-prep}^{C-WP}\)) pairs. Finally, \(\Delta E^{C-(T) }_{int}\) represents the triples correction contribution to the interaction energy between the fragments. In the LED scheme, this term can be further decomposed into intra- and interfragment components. This can be useful, for example, to estimate the London dispersion contribution from the triples correction term, as suggested in ref.[536].