7.4. Choice of Computational Model

7.4.1. Features Common to All Calculations

The computational model is specified in the block %method. The following choices exist:

%method
  Method  HFGTO  (or HF)   # Hartree-Fock with GTOs
          DFGTO  (or DFT)  # Density Functional Theory with GTOs
          MP2              # Second-order Møller-Plesset Perturbation Theory
          CNDO             # Complete neglect of differential overlap
          INDO             # Intermediate neglect of differential overlap
          NDDO             # Neglect of diatomic differential overlap
end

In the case of Hartree-Fock calculations [840], nothing else is required in this block. Density functional calculations [451, 647] need slightly more attention, as will be seen in the next section.

The type of calculation to be performed can be chosen by the RunTyp flag as follows:

%method
  RunTyp  Energy (or SinglePoint)     # Single point calculation (default)
          Scan   (or Trajectory)      # Geometric scan
          Gradient                    # Single point energy and gradient
          Opt    (or GeometryOpt)     # Geometry optimization
          MD                          # Molecular dynamics calculation
end

7.4.2. Density Functional Calculations

7.4.2.1. Choice of Functional

Basic Choice of Density Functional. If you are doing a DFT calculation [451, 647], the following choices for local, gradient-corrected, and hybrid density functionals are available. See also the simple input keywords in section Density Functional Methods for a more complete list of density functionals (including double-hybrids, which cannot be called from the %method block for technical reasons).

%method

  Functional
      #***************************************
      # Local functionals
      #***************************************
        HFS        # Hartree-Fock Slater (Slater exchange only)
        LSD        # Local spin density (VWN-5A form)
        VWN5       # Local spin density (VWN-5)
        VWN3       # Local spin density (VWN-3)
        PWLDA      # Local spin density (PW-LDA)
      #***************************************
      # ``Pure'' GGA functionals
      #***************************************
        BNULL      # Becke '88 exchange, no correlation
        BVWN       # Becke '88 exchange, VWN-5 correlation
        BP         # Becke '88 exchange, Perdew '86 correlation
        PW91       # Perdew-Wang (PW) GGA-II '91 functional
        mPWPW      # Modified PW exchange, PW correlation
        mPWLYP     # Modified PW exchange, Lee-Yang-Parr (LYP) correlation
        BLYP       # Becke '88 exchange, LYP correlation
        GP         # Gill '96 exchange, Perdew '86 correlation
        GLYP       # Gill '96 exchange, LYP correlation
        PBE        # Perdew-Burke-Ernzerhof (PBE) functional
        revPBE     # Revised PBE (exchange scaling)
        RPBE       # Revised PBE (modified exchange functional)
        PWP        # PW '91 exchange, Perdew '86 correlation
        OLYP       # Hoe/Cohen/Handy's optimized exchange, LYP correlation
        OPBE       # Hoe/Cohen/Handy's optimized exchange, PBE correlation
        XLYP       # Xu/Goddard exchange, LYP correlation
        B97D       # Grimme's GGA including D2 dispersion correction
        PW86PBE    # PW '86 exchange, PBE correlation (as used for vdw-DF and related)
        RPW86PBE   # Revised PW '86 exchange, PBE correlation
      #***************************************
      # Meta-GGA functionals
      #***************************************
        M06L       # Truhlar's semi-local functional
        TPSS       # TPSS functional
        revTPSS    # Revised TPSS functional
        SCANfunc   # Perdew's SCAN functional
        RSCAN      # Regularized SCAN functional
        R2SCAN     # Regularized and restored SCAN functional
      #***************************************
      # GGA Hybrid functionals
      #***************************************
        B1LYP      # 1-parameter hybrid of BLYP (25% HF exchange)
        B1P        # Similar with Perdew '86 correlation
        G1LYP      # 1-parameter analog with Gill '96 exchange
        G1P        # Similar with Perdew '86 correlation
        B3LYP      # 3-parameter Hybrid of BLYP (20% HF exchange)
        B3P        # Similar with Perdew '86 correlation
        G3LYP      # 3-parameter analog with Gill '96 exchange
        G3P        # Similar with Perdew '86 correlation
        PBE0       # 1-parameter version of PBE (25% HF exchange)
        PWP1       # 1-parameter version of PWP (analog of PBE0)
        mPW1PW     # 1-parameter version of mPWPW (analog of PBE0)
        mPW1LYP    # 2-parameter version of mPWLYP (analog of PBE0)
        PW91_0     # 1-parameter version of PW91 (analog of PBE0)
        O3LYP      # 3-parameter version of OLYP 
        X3LYP      # 3-parameter version of XLYP 
        B97        # Becke's original hybrid functional
        BHANDHLYP  # Becke's half-and-half hybrid functional (50% HF exchange)
      #***************************************
      # Meta-GGA Hybrid functionals
      #***************************************
        TPSSh      # TPSS hybrid with 10% HF exchange
        TPSS0      # TPSS hybrid with 25% HF exchange
        PW6B95     # Truhlar's 6-parameter hybrid functional 
        M06        # Truhlar's 2006 low-HF hybrid (27% HF exchange)
        M062X      # Truhlar's 2006 high-HF hybrid (54% HF exchange)
        r2SCANh    # r2SCAN global hybrid with 10% HF exchange
        r2SCAN0    # r2SCAN global hybrid with 25% HF exchange 
        r2SCAN50   # r2SCAN global hybrid with 50% HF exchange
      #***************************************
      # Range-Separated Hybrid functionals
      #***************************************
        wB97       # Head-Gordon's fully variable DF 
        wB97X      # Head-Gordon's DF with minimal HF exchange 
        CAMB3LYP   # Handy's fit 
        LC_BLYP    # Hirao's original application 
        LC_PBE     # Hirao's PBE-based range-separated hybrid 
        wr2SCAN    # r2SCAN range-separated hybrid 
end

Note that Functional is a compound keyword. It chooses specific values for the variables Exchange, Correlation, and ACM described below. If given as a simple input keyword, in some cases, it will also activate a dispersion correction. You can explicitly give these variables instead or in addition to Functional. However, make sure that you specify these variables after you have assigned a value to Functional or the values of Exchange, Correlation and ACM will be reset to the values chosen by Functional.

Empirical Parameters in Density Functionals. Some functionals incorporate empirical parameters that can be changed to improve agreement with experiment. Currently, there are a few parameters that can be changed (other than the parameters used in the hybrid functionals, which are described later).

The first of these parameters is \(\alpha\) of Slater’s X\(_{\alpha }\) method. Theoretically, it has a value of 2/3 and this is used in the HFS and LSD functionals. However, the exchange contribution is underestimated by about 10% by this approximation (quite significant!) and a value around 0.70-0.75 is recommended for molecules.

The second parameter is the \(\beta\) for Becke’s gradient-corrected exchange functional. Becke determined the value 0.0042 by fitting the exchange energies for rare gas atoms. There is some evidence that with smaller basis sets, a slightly smaller value such as 0.0039 gives improved results for molecules.

The next parameter is the value \(\kappa\), which occurs in the PBE exchange functional. It has been given the value 0.804 by Perdew et al. in order to satisfy the Lieb-Oxford bound. Since then, other workers have argued that a larger value for this parameter (around 1.2) gives better energetics, which is explored in the revPBE functional. Note, however, that while revPBE gives slightly better energetics, it also gives slightly poorer geometries.

The last two parameters are also related to PBE. Within the PBE correlation functional, there is \(\beta_C\) (not to be confused with the \(\beta\) exchange parameter in Becke’s exchange functional), whose original value is \(\beta_C=0.066725\). Modified variants exist with different \(\beta_C\) values, e.g., the PBEsol functional and the PBEh-3c compound method. Furthermore, the \(\mu\) parameter in the PBE exchange functional may be modified. In the original formulation, it is related to \(\beta_C\) via \(\mu = \beta_C \frac{\pi^{2} }{3}\), but has been modified in later variants as well.

%method
  XAlpha   0.75       # Slater's alpha parameter         (default 2/3)
  XBeta    0.0039     # Becke's beta parameter           (default 0.0042)
  XKappa   0.804      # PBE(exchange) kappa parameter    (default 0.804)
  XMuePBE  0.21952    # PBE(exchange) mue parameter      (default 0.21952)
  CBetaPBE 0.066725   # PBE(correlation) beta parameter  (default 0.066725)
end

Specifying Exchange and Correlation approximations individually. The following keywords are available for specifying the exchange and correlation approximations individually. In addition, scaling parameters can be defined to construct user-defined hybrid or “extended” hybrid functionals:

%method

  Exchange    
      X_NOX      # No exchange
      X_SLATER   # Slater's local exchange 
      X_B88      # Becke '88 gradient exchange 
      X_wB88     # Short-range Becke '88 exchange for range-separated functionals
      X_G96      # Gill '96 gradient exchange 
      X_PW91     # Perdew-Wang (PW) '91 gradient exchange  
      X_mPW      # Adamo-Barone modification of PW exchange
      X_PBE      # Perdew-Burke-Ernzerhof (PBE) exchange
      X_RPBE     # Revised PBE exchange
      X_OPTX     # Hoe/Cohen/Handy's optimized exchange 
      X_X        # Xu/Goddard exchange
      X_TPSS     # TPSS meta-GGA exchange 
      X_B97D     # Grimme's modified exchange for the B97-D GGA 
      X_B97BECKE # Becke's original exchange for the B97 hybrid 
      X_SCAN     # Perdew's constrained exchange for the SCAN meta-GGA 
      X_RSCAN    # Perdew's constrained exchange for the rSCAN meta-GGA 
      X_R2SCAN   # Perdew's constrained exchange for the r2SCAN meta-GGA 
  
  Correlation 
      C_NOC      # No correlation
      C_VWN5     # Local Vosko-Wilk-Nusair correlation (parameter set "V")
      C_VWN3     # Local Vosko-Wilk-Nusair correlation (parameter set "III") 
      C_PWLDA    # Local PW correlation
      C_P86      # Perdew '86 correlation 
      C_PW91     # PW '91 correlation 
      C_PBE      # PBE correlation 
      C_LYP      # LYP correlation 
      C_TPSS     # TPSS meta-GGA correlation 
      C_B97D     # Grimme's modified correlation for the B97-D GGA 
      C_B97BECKE # Becke's original correlation for the B97 hybrid 
      C_SCAN     # Perdew's constrained correlation for the SCAN meta-GGA 
      C_RSCAN    # Perdew's constrained correlation for the rSCAN meta-GGA 
      C_R2SCAN   # Perdew's constrained correlation for the r2SCAN meta-GGA 
  
  # Parameters for hybrid functionals
  ACM   ACM-A, ACM-B, ACM-C
                 # ACM-A: fraction of HF exchange in hybrid DFT
                 # ACM-B: scaling of GGA exchange part of DFT
                 # ACM-C: scaling of GGA correlation part of DFT
  
  # Parameters for "extended" hybrid functionals
  ScalLDAC 1.0   # scaling of the LDA correlation part
  ScalMP2C 0.0   # fraction of MP2 correlation mixed into the density functional

end

Hybrid Density Functionals. The hybrid DFs [87, 88] are invoked by choosing a nonzero value for the variable ACM. (ACM stands for “adiabatic connection model”). Specifically, these functionals have the following form:

(7.1)\[E_{\text{XC} } = a E_{\text{HF} }^{\text{X} } +\left( 1 - a \right)E_{\text{LSD} }^{\text{X} } +bE_{\text{GGA} }^{\text{X} } +E_{\text{LSD} }^{\text{C} } +cE_{\text{GGA} }^{\text{C} } \]

Here, \(E_{\text{XC} }\) is the total exchange/correlation energy, \(E_{\text{HF} }^{\text{X} }\) is the Hartree-Fock exchange, \(E_{\text{LSD} }^{\text{X} }\) is the local (Slater) exchange, \(E_{\text{GGA} }^{\text{X} }\) is the gradient correction to the exchange, \(E_{\text{LSD} }^{\text{C} }\) is the local, spin-density based part of the correlation energy, and \(E_{\text{GGA} }^{\text{C} }\) is the gradient correction to the correlation energy.

This brings us to a slightly awkward subject: several hybrid functionals with the same name give different values in different programs. The reason for this is that they choose slightly different default values for the parameters \(a\), \(b\), and \(c\) and/or differ in the way they treat the local part of the correlation energy.

Different parametrizations exist. The most popular example of this is due to Vosko, Wilk, and Nusair (VWN, [874]). VWN in their classic paper give two sets of parameters — one in the main body (parametrization of RPA results; known as VWN-III) and one in their table 5 of correlation energies (parametrization of the Ceperley/Alder Monte Carlo results; known as VWN-V). Some programs use VWN-III, while others use VWN-V. Additionally, a slightly better fit to the uniform electron gas has been produced by Perdew and Wang [668]. The results from this fit are very similar to those using VWN5, whereas VWN3 yields quite different values.

In ORCA, almost all functionals choose PWLDA as the underlying LDA functional. A special situation arises if LYP is the correlation functional [502] since LYP is not a correction to the correlation, but rather includes the full correlation. It is therefore used in the B3LYP functional as:

(7.2)\[E_{\text{B3LYP} }^{\text{C} } =E_{\text{LSD} }^{\text{C} } +c\left({ E_{\text{LYP} }^{\text{C} } -E_{\text{LSD} }^{\text{C} } } \right)\]

In ORCA, VWN5 is chosen for the local correlation part of B3LYP. This choice is consistent with the TurboMole program [12, 13, 14], but not with the Gaussian program [278]. However, the user has full control over this setting. The underlying local part of any correlation functional can be set in the input with the variable LDAOpt:

%method
  LDAOpt  C_PWLDA
          C_VWN5
          C_VWN3
end

Specifying C_VWN3 for LDAOpt together with Functional=B3LYP should give results very close to the B3LYP functional as implemented in the Gaussian series of programs[1]. Due to the popularity of the B3LYP functional, the following aliases are defined in order to facilitate comparisons with other major electronic structure packages:

%method
  Functional B3LYP     # consistent with TurboMole
             B3LYP_TM  #    := Functional B3LYP
                       #       LDAOpt     C_VWN5
             B3LYP_G   # consistent with Gaussian
                       #    := Functional B3LYP
                       #       LDAOpt     C_VWN3
end

One-Parameter Hybrid Density Functionals. The underlying LDA-dependence of the three-parameter (i.e. ACM) hybrids causes slightly different results from programs which use a different choice of LDA. It has been argued from a theoretical viewpoint that the optimal mixing of HF exchange is 25% [248]. It has been further shown that use of this fixed ratio and not scaling the GGA correlation or exchange contributions gives results that are as good as the original three-parameter hybrids [9]. The one-parameter hybrid PBE0 has been advertised as a hybrid functional of overall well-balanced accuracy [8].

As such, the one-parameter hybrids are based more in theory and remove some arbitrariness from the hybrid procedures. The slightly higher HF-exchange (0.25 in favor of 0.20 used in the original three-parameter hybrids) is also reasonable for many systems. The one-parameter hybrids have the simple form:

(7.3)\[E_{\text{XC}} = E_{\text{DFT}}^{\text{X}} + a' \left( E_{\text{HF}}^{\text{X}} - E_{\text{DFT}}^{\text{X}} \right) + E_{\text{DFT}}^{\text{C}} \]

with \(a' = \frac{1}{4}\). This is the same as the three-parameter hybrids (equation (7.1)) with \(a = a'\), \(b = 1 - a'\), and \(c = 1\) and is how it is implemented.

Extended “double-hybrid” functionals. In addition to mixing the HF-exchange into a given DF, Grimme has proposed to mix in a fraction of the MP2 correlation energy as calculated with hybrid DFT orbitals [320]. Such functionals may be referred to as “extended” or “double” hybrid functionals. Grimme’s expression is:

(7.4)\[E_{\text{XC}} = aE_{\text{X}}^{\text{HF}} +\left({1-a} \right)E_{\text{X}}^{\text{DFT}} +\left({1-c} \right)E_{\text{C}}^{\text{DFT}} +cE_{\text{C}}^{\text{MP2}} \]

Such functionals can be defined by the user in ORCA as follows:

%method
  ScalHFX  = a
  ScalDFX  = 1-a
  ScalGGAC = 1-c
  ScalLDAC = 1-c
  ScalMP2C = c
end

Grimme recommends the B88 exchange functional, the LYP correlation functional and the parameters \(a = 0.53\) and \(c = 0.27\). This gives the B2PLYP functional which appears to be a fair bit better than B3LYP based on Grimme’s detailed evaluation study.

Presently, this methodology covers single point calculations, analytic gradients (and thus also geometry optimizations, relaxed scans, and transition state searches), and frequencies and other second-derivative properties (without the frozen core approximation in the MP2 part). Note that %mp2 density relaxed end should be specified in order to get the response density that is consistent with first-order properties from analytic first-derivatives. By default, this response density is not calculated since its construction adds significant overhead to the calculation. Therefore, you have to specifically request it to look at the consistent density. The considerably less costly unrelaxed (expectation value-like) density may instead requested by %mp2 density unrelaxed end. However, this is not recommended since the changes to the relaxed density are considerable in our experience and the unrelaxed density has a much weaker theoretical status than its relaxed counterpart.

Range-separated hybrid functionals. ORCA supports functionals based on the error function splitting of the two-electron operator used for exchange as first realized by Hirao and coworkers [410]:

(7.5)\[r_{12}^{-1} = \underbrace{\text{erfc}(\mu \cdot r_{12}) \cdot r_{12}^{-1} }_{\text{SR} } + \underbrace{\text{erf}(\mu \cdot r_{12}) \cdot r_{12}^{-1} }_{\text{LR} }\]

where \(\text{erf}(x) = \frac{2}{\sqrt{\pi} } \int_{0}^{x} \exp(-t^{2}) dt\) and \(\text{erfc}(x) = 1 - \text{erf}(x)\). Note that the splitting is only applied to the exchange contribution; all other contributions (one-electron parts of the Hamiltonian, the electron-electron Coulomb interaction and the approximation for the DFT correlation) are not affected. Later, Handy and coworkers generalized the ansatz to [900]:

(7.6)\[r_{12}^{-1} = \underbrace{\frac{1 - [ \alpha + \beta \cdot \text{erf}(\mu \cdot r_{12})]}{ r_{12} }}_{\text{SR} } + \underbrace{\frac{\alpha + \beta \cdot \text{erf}(\mu \cdot r_{12}) }{r_{12} }}_{\text{LR} }\]

This form of the splitting used in ORCA is shown visually (according to Handy and coworkers) in Fig. 7.1.

Fig. 7.1 Graphical description of the Range-Separation ansatz. The gray area corresponds to Hartree-Fock exchange. \(\alpha\) and \(\beta\) follow Handy’s terminology [900].

The splitting has been used to define the \(\omega\)B97 family of functionals, wherein the short-range part (SR) is described by DFT exchange and the long-range part by exact exchange, i.e. Hartree-Fock exchange. The same is true for CAM-B3LYP, LC-BLYP, and LC-PBE. It is possible to use a fixed amount of Hartree-Fock exchange (EXX) and/or a fixed amount of DFT exchange in this ansatz. Here are some examples of range-separated hybrid functionals available in ORCA and their parameters:

Functional	Keyword	fixed EXX	variable part	\(\mu\)/\(\text{bohr}^{-1}\)	fixed DFT-X	Reference
\(\omega\)B97	WB97	—	100%	0.40	—	[151]
\(\omega\)B97X	WB97X	15.7706%	84.2294%	0.30	—	[151]
\(\omega\)B97X-D3	WB97X-D3	19.5728%	80.4272%	0.25	—	[525]
\(\omega\)B97X-V	WB97X-V	16.7%	83.3%	0.30	—	[555]
\(\omega\)B97X-D3BJ	WB97X-D3BJ	16.7%	83.3%	0.30	—	[603]
CAM-B3LYP	CAM-B3LYP	19%	46%	0.33	35%	[900]
LC-BLYP	LC-BLYP	—	100%	0.33	—	[846]
LC-PBE	LC-PBE	—	100%	0.47	—	[410]

The currently available speed-up options are RIJONX and RIJCOSX. Integral-direct single-point calculations, calculations involving the first nuclear gradient (i.e. geometry optimizations), frequency calculations, TDDFT, TDDFT nuclear gradient, and EPR/NMR calculations are the supported job types with range-separated hybrid functionals thus far.

In principle, it is possible to change the amount of fixed Hartree-Fock exchange (ACM-A), the amount of variable exchange (RangeSepScal), and \(\mu\), though this is not recommended. The amount of fixed DFT Exchange (ACM-B) can only be changed for CAM-B3LYP and LC-BLYP. In other words, ACM-B is ignored by the \(\omega\)B97 approaches as no corresponding \(\mu\)-independent exchange functional has been defined.

! RKS CAM-B3LYP DEF2-SVP

# default parameters for CAM-B3LYP:
%method
  RangeSepEXX  True    # must be set 
  RangeSepMu   0.33    # should not be set to 0.0 or below
  RangeSepScal 0.46    # should sum to 1 with ACM-A and ACM-B
  ACM 0.19, 0.35, 0.81 # ACM-A, ACM-B, ACM-C(same as B3LYP)
end

* xyz 0 1
H 0.0 0.0 0.0
H 0.0 0.0 0.7
*

Note

For information on the ACM formalism, see preceding section called “Specifying Exchange and Correlation approximations individually”. While it is technically possible to choose an exchange functional that has no \(\mu\)-dependence, this makes no sense conceptually.

7.4.2.2. LibXC Functionals

Since ORCA 4.2, it is possible to use the functionals provided by LibXC[2] within the ORCA framework. The LibXC version used by ORCA is printed at the beginning of the output. The LibXC reference should be cited when LibXC is used as part of your calculations. For reference, see [506].

The complete list of functionals available via the LibXC interface can always be inspected by typing at the command line

orca -libxcfunctionals

The list of functionals has the following form:

Functionals available via LibXC:
  No.:     ID / Keyword                    - Name
    0:     18 / lda_c_1d_csc               - Casula, Sorella & Senatore
    1:     26 / lda_c_1d_loos              - P-F Loos correlation LDA
    2:     15 / lda_c_2d_amgb              - AMGB (for 2D systems)
    3:     16 / lda_c_2d_prm               - PRM (for 2D systems)
    4:    552 / lda_c_br78                 - Brual & Rothstein 78

   ...

The list is grouped by “level” of functional (LDA, GGA, meta-GGA, etc.) and then by part of the energy it models (correlation, exchange, exchange-correlation). Correlation functionals carry a ’_c_’ in their names, exchange functionals an ’_x_’, and combined exchange-correlation functionals an ’_xc_’. Additional information for a specific functional can be requested using

orca -libxcinfo [functional name]

where the functional name is the keyword in the above list.

Specification of LibXC functionals in the ORCA input follows the standard format:

%method 
  method     dft
  functional hyb_gga_xc_b3lyp
end

or in the case of separate exchange and correlation specifications:

%method 
  method      dft
  exchange    mgga_x_m06_l
  correlation mgga_c_m06_l
end

CAM-type range-separated functionals are supported through the LibXC interface since ORCA 5.0. So are functionals that include non-local (VV10) correlation (see DFT Calculations with the Non-Local, Density Dependent Dispersion Correction (VV10): DFT-NL). The VV10 part is computed internally by ORCA. Other non-local functionals, such as BEEF-vdW, are not supported. Meta-GGA functionals that depend on the kinetic energy density \(\tau\) are supported, but not those that depend on the Laplacian of the density \(\nabla^2\rho\).

Double-hybrid functionals can be constructed with LibXC components as described in section DFT Calculations with Second Order Perturbative Correction (Double-Hybrid Functionals), but only with separate exchange and correlation components. So exchange=gga_x_pbe and correlation=gga_c_pbe can be used, but functional=hyb_gga_xc_pbeh cannot be used in a double-hybrid formulation. Beware that the exchange and correlation contributions calculated by LibXC are simply scaled by ScalDFX and ScalGGAC, respectively, and no care is taken to separately scale LDA components or alter other internal parameters!

7.4.2.2.1. Screening Thresholds

When the density is smaller than a certain threshold, LibXC skips the evaluation of the functional and instead sets all the output quantities to zero in order to avoid under- and/or overflows. The default thresholds for different functionals are set by the LibXC developers, but may sometimes be too low. We have not performed extensive testing, but allow the user to set the threshold. Similarly, there are thresholds for minimum values of \(\zeta\) and \(\tau\) in order to avoid division by zero. The default values are functional-independent and can be changed using the following keywords.

%method 
  LibXCDensityThreshold 1e-12    # seems to be reasonable
  LibXCZetaThreshold    2e-16    # default value in LibXC
  LibXCTauThreshold     1e-20    # default value in LibXC
end

7.4.2.2.2. Modifying LibXC Functional Parameters

Starting with ORCA 5.1 it is possible to modify the “external parameters” of a functional that the LibXC interface provides. The available parameters and their default values can be seen in the output or with the orca -libxcinfo command above. Please exercise caution when using this interface and if using a published functional reparametrization, make sure you can reproduce results from the original publication. The syntax requires one of the ExtParamX, ExtParamC, or ExtParamXC keywords (depending on which functional is modified), followed by the parameter name in quotation marks, and finally the new value. For example, here is an input for the PWPB95 double-hybrid functional, where the simple input keyword is used to set most parameters (such as the MP2 scaling factors) and only the exchange and correlation parameters of the mPW91 and B95 LibXC functionals, respectively, are modified.

! Opt PWPB95 def2-SVP def2-SVP/C

%method
  Exchange    gga_x_mpw91
  Correlation mgga_c_bc95
  ExtParamX   "_bt"     0.00444
  ExtParamX   "_alpha" 19.823391
  ExtParamX   "_expo"   3.7868
  ExtParamC   "_css"    0.03241
  ExtParamC   "_copp"   0.00250
end

*xyz 0 1
  H 0 0 0
  F 0 0 0.9
*

Note that the variable definitions in LibXC may be different from the ones used internally in ORCA or in the original publication, due to various constant factors, etc. When in doubt, please consult the LibXC documentation or source code repository — we simply provide access to the interface.

7.4.2.2.3. Simple Input of LibXC Functionals

Some LibXC functionals are accessible via the simple input keyword !LibXC(Keyword), e.g. !LibXC(TPSS). The keywords match those in section Density Functional Methods for functionals that are also implemented natively in ORCA. Using this syntax, parameters for the DFT-D dispersion corrections are also set automatically (if implemented). Table Table 7.4 lists the available functionals and their LibXC names.

Table 7.4 LibXC functionals available via the simple input !LibXC(Keyword)

Keyword	LibXC name(s)	Notes
GGA functions
B97-D	gga_xc_b97_d	Uses D2
B97-D3	gga_xc_b97_d	Uses D3(0)
B97-D4	gga_xc_b97_d	Uses D4
BLYP	gga_x_b88 + gga_c_lyp
BP86	gga_x_b88 + gga_c_p86
KT2	gga_xc_kt2
KT3	gga_xc_kt3
PBE	gga_x_pbe + gga_c_pbe
Meta-GGA functionals
B97M-V	mgga_xc_b97m_v	Uses VV10
B97M-D3BJ	mgga_xc_b97m_v	Uses D3(BJ)
B97M-D4	mgga_xc_b97m_v	Uses D4
M06L	mgga_x_m06_l + mgga_c_m06_l
MN15L	mgga_x_mn15_l + mgga_c_mn15_l
r2SCAN	mgga_x_r2scan + mgga_c_r2scan
rSCAN	mgga_x_rscan + mgga_c_rscan
SCAN	mgga_x_scan + mgga_c_scan
TPSS	mgga_x_tpss + mgga_c_tpss
Hybrid GGA functionals
B3LYP	hyb_gga_xc_b3lyp5
B3LYP/G	hyb_gga_xc_b3lyp
B97	hyb_gga_xc_b97
PBE0	hyb_gga_xc_pbeh
wB97	hyb_gga_xc_wb97
wB97X	hyb_gga_xc_wb97x
wB97X-D3	hyb_gga_xc_wb97x_d3	Uses D3(0)
wB97X-V	hyb_gga_xc_wb97x_v	Uses VV10
wB97X-D3BJ	hyb_gga_xc_wb97x_v	Uses D3(BJ)
wB97X-D4	hyb_gga_xc_wb97x_v	Uses D4
Hybrid meta-GGA functionals
M05	hyb_mgga_x_m05 + mgga_c_m05
M052X	hyb_mgga_x_m05_2x + mgga_c_m05_2x
M06	hyb_mgga_x_m06 + mgga_c_m06
M062X	hyb_mgga_x_m06_2x + mgga_c_m06_2x
MN15	hyb_mgga_x_mn15 + mgga_c_mn15
PW6B95	hyb_mgga_xc_pw6b95
SCAN0	hyb_mgga_x_scan0 + mgga_c_scan
TPSS0	hyb_mgga_xc_tpss0
TPSSH	hyb_mgga_xc_tpssh
wB97M-V	hyb_mgga_xc_wb97m_v	Uses VV10
wB97M-D3BJ	hyb_mgga_xc_wb97m_v	Uses D3(BJ)
wB97M-D4	hyb_mgga_xc_wb97m_v	Uses D4
Double-hybrid functionals
DSD-PBEB95-D3	gga_x_pbe + mgga_c_bc95	Custom mixing, uses D3(BJ)
DSD-PBEB95-D4	gga_x_pbe + mgga_c_bc95	Custom mixing, uses D4
PWPB95	gga_x_mpw91 + mgga_c_bc95	Custom mixing and ExtParam’s

7.4.2.3. Using the RI-J Approximation to the Coulomb Part

Note

This is the default for non-hybrid DFT! Can be turned off by using !NORI.

A very useful approximation that greatly speeds up DFT calculations unless the molecule gets very large is the so called “RI-approximation” [65, 225, 245, 246, 440, 862, 888]. RI stands for “Resolution of the identity”. In short, charge distributions arising from products of basis functions are approximated by a linear combination of auxiliary basis functions.

(7.7)\[\phi_{i} \left({ \vec{{r} }} \right)\phi_{j} \left({ \vec{{r} }} \right)\approx \sum\limits_k { c_{k}^{ij} \eta_{k} (\mathrm{\mathbf{r} }) } \]

There are a variety of different possibilities to determine the expansion coefficients \(c_{k}^{ij}\). A while ago, Almlöf and coworkers [860] have shown that for the approximation of electron repulsion integrals, the best choice is to minimize the residual repulsion[3].

Define:

(7.8)\[R_{ij} \equiv \phi_{i} \left({ \vec{{r} }} \right)\phi_{j} \left({\vec{{r} }} \right)-\sum\limits_k { c_{k}^{ij} \eta_{k} \left({ \vec{{r} }} \right)} \]

and

(7.9)\[T_{ij} =\int{ \int{ R_{ij} \left({ \vec{{r} }} \right)\frac{1}{\left|{\vec{{r} }-{\vec{{r} }}'} \right|}R_{ij} \left({ \vec{{r} }} \right)d^{3}rd^{3}{r}'} } \]

Determining the coefficients that minimize \(T_{ij}\) leads to

(7.10)\[\mathrm{\mathbf{c} }^{ij} =\mathrm{\mathbf{V} }^{-1}\mathrm{\mathbf{t} }^{ij} \]

where:

(7.11)\[t_{k}^{ij} =\left\langle { \phi_{i} \phi_{j} \left|{ r_{12}^{-1} } \right|\eta_{k} } \right\rangle\]

(7.12)\[V_{ij} =\left\langle { \eta_{i} \left|{ r_{12}^{-1} } \right|\eta_{j} }\right\rangle\]

Thus, an ordinary two-electron integral becomes

(7.13)\[\begin{split}\begin{aligned} \left\langle { \phi_{i} \phi_{j} \left|{ r_{12}^{-1} } \right|\phi_{k} \phi_{l} } \right\rangle &\approx \sum\limits_{p,q} { c_{p}^{ij} c_{q}^{kl} V_{pq} } \\ &=\sum\limits_{p,q} {V_{pq} \sum\limits_r { \left({ \mathrm{\mathbf{V} }^{-1} } \right)_{pr} t_{r}^{ij} } \sum\limits_s { \left({ \mathrm{\mathbf{V} }^{-1} } \right)_{qs} t_{s}^{kl} } } \\ &=\sum\limits_{r,s} {\left({ \mathrm{\mathbf{V} }^{-1} } \right)_{rs} t_{r}^{ij} t_{s}^{kl} } \end{aligned} \end{split}\]

and the total Coulomb energy becomes

(7.14)\[\begin{split}\begin{aligned} E_{J} &=\sum\limits_{i,j} { \sum\limits_{k,l} { P_{ij} P_{kl} \left\langle {\phi_{i} \phi_{j} \left|{ r_{12}^{-1} } \right|\phi_{k} \phi_{l} }\right\rangle} } \\ &\approx \sum\limits_{i,j} { \sum\limits_{k,l} { P_{ij} P_{kl} \sum\limits_{r,s} { \left({ \mathrm{\mathbf{V} }^{-1} } \right)_{rs} t_{r}^{ij} t_{s}^{kl} } } } \\ &=\sum\limits_{r,s} { \left({ \mathrm{\mathbf{V} }^{-1} } \right)_{rs} } \underbrace{ \sum\limits_{i,j} { P_{ij} t_{r}^{ij} } }_{\mathrm{\mathbf{X} }_{r}} \underbrace{ \sum\limits_{k,l} { P_{kl} t_{s}^{kl} } }_{\mathrm{\mathbf{X} }_{s}} \end{aligned} \end{split}\]

where \(\mathrm{\mathbf{P} }\) is the total density matrix.

In a similar way, the Coulomb contribution to the Kohn-Sham matrix is calculated. There are substantial advantages from this approximation: the quantities to be stored are the matrix \(\mathrm{\mathbf{V} }^{-1}\) — which depends only on two indices — and the three-index auxiliary integrals \(t_{r}^{ij}\). This leads to a tremendous reduction of storage requirements compared to a four-index list of repulsion integrals. In addition, the two- and three-index electron repulsion integrals are easier to compute than the four-index integrals, leading to further reductions in processing time. Furthermore, the Coulomb energy and the Kohn-Sham matrix contributions can be quickly assembled by simple vector/matrix operations, leading to large time savings. This arises because each auxiliary basis function \(\eta_{k} \left({ \vec{{r} }} \right)\) appears in the expansion of many charge distributions \(\phi_{i} \left({\vec{{r} }} \right)\phi_{j} \left({ \vec{{r} }} \right)\). Unfortunately, a similar strategy is less easily applied (or, at least, with less benefit) to the Hartree-Fock exchange term.

If the auxiliary basis set \(\left\{\eta \right\}\) is large enough, the approximation is also highly accurate. Since any DFT procedure already has a certain, sometimes sizable, error from the noise in the numerical integration of the XC part, it might be argued that a similarly large error in the Coulomb part is perfectly acceptable without affecting the overall accuracy of the calculation much. Furthermore, the errors introduced by the RI method are usually much smaller than the errors in the calculation due to basis set incompleteness. It is therefore recommended to use the RI procedure for pure DFs. However, one should probably not directly mix absolute total energies obtained from RI and non-RI calculations as the error in the total energy accumulates and will rise with increased molecular size, while the errors in the relative energies will tend to cancel.

There are several choices for auxiliary basis sets described in the next section, which depend on the choice of the primary GTO basis set used to expand the molecular orbitals[4].

In ORCA, the RI approximation is toggled by the input

%method
  RI  on   # do use the RI-J approximation
      off  # do not use the RI-J approximation
end

Note

If you use RI, you must specify an auxiliary basis set (in the %basis section or using the appropriate simple keyword) or use the !AutoAux simple keyword.

7.4.2.4. The Split-RI-J Coulomb Approximation

There is an improved version of the RI-algorithm that has been implemented since ORCA 2.2.09. This Split-RI-J algorithm yields the same Coulomb energy as the standard RI-algorithm, but is significantly faster if the basis set contains many high angular momentum functions (d-, f-, g-functions). For small basis sets, there is virtually no difference between the two algorithms, except that Split-RI-J uses more memory than standard RI. However, calculations with ca. 2000 basis functions only need about an extra 13 MB for Split-RI-J, which is a trivial requirement on present-day hardware.

The Split-RI-J algorithm is invoked with the !Split-RI-J simple keyword. Split-RI-J is presently only available for SCF and gradient calculations.

Note

• The Split-RI-J algorithm is the default if RI is turned on via !RI. If you do not want to use Split-RI-J, please also use the keyword !NoSplit-RI-J

7.4.2.5. Using the RI Approximation for Hartree-Fock and Hybrid DFT (RIJONX)

The RI approximation can be used, although with less benefit, for hybrid DFT and Hartree-Fock (RHF and UHF) calculations. In this case, a different algorithm[5] is used that allows a fair approximation to the Hartree-Fock exchange matrix (RI-JK). It can be difficult to make this approximation highly accurate. It is, however, usefully fast compared to direct SCF if the molecule is “dense” enough. There are special auxiliary basis sets for this purpose (see section Basis Sets).

%method
  RI  on
end

%basis
  Aux "def2/JK" 
end

Note

There has been little experimentation with this feature. It is provided on an experimental basis here.

Alternatively, the RI method can be used for the Coulomb term and the standard treatment for the exchange term. This method is called RIJONX since the exchange term should tend towards linear scaling for large molecules. This feature can be used for Hartree-Fock and hybrid DFT calculations by using:

%method
  RI      on      # do use the RI approximation
  RIFlags 1       # ...but treat exchange exactly
end

# Equivalently, use the following simple keyword:
! RIJONX

# Remember to assign an auxiliary basis set!

The requirements for the auxiliary basis are the same as for the normal RI-J method.

7.4.2.6. Using the RI Approximation for Hartree-Fock and Hybrid DFT (RIJCOSX)

Note

This is the default for hybrid DFT! Can be turned off by using !NOCOSX.

The aim of this approximation is to efficiently compute the elements of exchange-type matrices [6]:

(7.15)\[K_{\mu \nu } =\sum\limits_{\kappa \tau } { P_{\kappa \tau } } (\mu \kappa \vert \nu \tau ) \]

where \(\mathrm{\mathbf{P} }\) is some kind of density-type matrix (not necessarily symmetric) and the two-electron integrals are defined over the basis set \(\{\varphi \}\) by:

(7.16)\[(\mu \kappa \vert \nu \tau )=\int{ \mu (\mathrm{\mathbf{r} }_{1} )\kappa (\mathrm{\mathbf{r} }_{1} )\nu (\mathrm{\mathbf{r} }_{2} )\tau (\mathrm{\mathbf{r} }_{2} )r_{12}^{-1} } d\mathrm{\mathbf{r} }_{1} d\mathrm{\mathbf{r} }_{2} \]

The approximation pursued here can be written as follows:

(7.17)\[K_{\mu \nu } \approx \sum\limits_g { X_{\mu g} } \sum\limits_\tau{A_{\upsilon \tau } (\mathrm{\mathbf{r} }_{g} ) } \sum\limits_\kappa{ P_{\kappa \tau } X_{\kappa g} } \]

Here, the index \(g\) refers to grid points \(\mathrm{\mathbf{r} }_{g}\) and:

(7.18)\[X_{\kappa g} =w_{g}^{1/2} \kappa (\mathrm{\mathbf{r} }_{g} ) \]

(7.19)\[A_{\upsilon \tau } (\mathrm{\mathbf{r} }_{g} )=\int{ \frac{\nu (\mathrm{\mathbf{r} })\tau (\mathrm{\mathbf{r} }) }{\vert \mathrm{\mathbf{r} }-\mathrm{\mathbf{r} }_{g} \vert }d\mathrm{\mathbf{r} }} \]

where \(w_{g}\) denotes the grid weights. Thus, the first integration is carried out numerically and the second one analytically. Note that this destroys the Hermitian character of the two-electron integrals.

Equation (7.17) is perhaps best evaluated in three steps:

(7.20)\[F_{\tau g} =(\mathrm{\mathbf{PX} })_{\tau g} \]

(7.21)\[G_{\nu g} =\sum\limits_\tau{ A_{\nu \tau } (\mathrm{\mathbf{r} }_{g} )F_{\tau g} } \]

(7.22)\[K_{\mu \nu } =(\mathrm{\mathbf{XG} }^{+})_{\mu \nu } \]

As such, the equations are very similar to the pseudo-spectral method extensively developed and discussed by Friesner and co-workers since the mid 1980s and commercially available in the Jaguar quantum chemistry package. The main difference at this point is that instead of \(X_{\kappa g}\) there appears a least-square fitting operator \(Q_{\kappa g}\) in Friesner’s formulation. Note that an analogue of the fitting procedure has also been implemented in ORCA and — in contrast to Friesner’s pseudo-spectral method — does not need specially optimized grids. The basic idea is to remove the grid errors within the basis set by “fitting” the numerical overlap to the analytical one. Due to its nature, overlap fitting is supposed to work better with larger basis sets.

7.4.2.6.1. The RIJCOSX gradient

Given the exchange matrix, the exchange energy is given by (a sum over spin cases is left out here for simplicity):

(7.23)\[E_{\text{X} } =\frac{1}{2}\sum\limits_{\mu \nu } { P_{\mu \nu } K_{\mu \nu } (\mathrm{\mathbf{P} }) } \]

Previous to ORCA6, the gradient of the COSX contribution to the energy was taken as an approximation:

(7.24)\[\frac{\partial E_{X} }{\partial \lambda }\approx 2\sum\limits_g {\sum\limits_{\mu \nu } { \frac{\partial F_{\mu g} }{\partial \lambda } } } G_{\nu g} \]

with

(7.25)\[\frac{\partial F_{\mu g} }{\partial \lambda }=w_{g}^{1/2} \sum\limits_\kappa{P_{\kappa \mu } \frac{\partial X_{\mu g} }{\partial \lambda } } \]

as published in the original implementation paper [624].

Starting from ORCA6, this was updated to the full derivative of the energy component, including the derivative of all terms: grid weight derivatives, derivative of the SFitting matrices and all derivatives of \(F\) and \(G\) [7]. The gradient is thus now more accurate and less noisy. In case one wants to revert to the previous approximated version (not recommended), just set:

%method
  cosxgradtype grad_n
  useqgradfit false
end

7.4.2.6.2. Working with the COSX Numerical Grids

For expert users, the grid parameters for the exchange grids can be even more finely controlled:

%method
  IntAccX  Acc1, Acc2, Acc3
  GridX    Ang1, Ang2, Ang3
end

There are three grids involved: the smallest grid (Acc1, Ang1) that is used for the initial SCF iterations, the medium grid (Acc2, Ang2) that is used until the end of the SCF and the largest grid (Acc3, Ang3) that is used for the final energy and the gradient evaluations. UseFinalGridX can be used to turn this last grid on or off, though changing this is not generally recommended. More details about the grid constructions can be found in Details on the numerical integration grids.

7.4.2.6.3. Some SFitting Parameters

To modify the overlap fitting parameters, the following input can be specified:

%method
  UseSFitting false   # Same as NoSFitting in the simple input
                      # (Default is true)
  UseQGradFit false   # Uses the SCF fitting matrix for gradient calculations
                      # (Default is true)
end

Note that overlap fitting also works for HF and MP2 gradients without specifying any additional keyword. The UseQGradFit parameter uses the same fitting matrix for the gradients as for the energy calculation and is the default behavior since ORCA6.

7.4.2.6.4. Use of Partial Contraction Scheme

Since ORCA 5.0, generally-contracted basis sets can be handled efficiently by using an intermediate partially contracted (pc) atomic-orbital basis for the exchange-matrix computation without affecting the results [383]. Depending on the basis set and element type, computational speedups by many orders of magnitude are possible. For testing or benchmark purposes, the K matrix computation can be done in the original basis by using the flag

%method
  COSX_PartialContraction false   # No intermediate basis for generally contracted shells
                                  # (Default is true)
end

7.4.2.6.5. Restoring Full Symmetry

The semi-numerical integration scheme in the default COSX algorithm breaks the permutational symmetry of the two-electron integrals. We have observed that this flaw is often the cause of convergence problems for iterative algorithms, in particular for multi-reference theories [382]. An input option is available since ORCA 6.0 to preserve the full eight-fold permutational symmetry of the two-electron integrals:

%method
  COSX_IntSymmetry Full      # Fully symmetrized integrals
                   Standard  # Original COSX algorithm
end

The full-symmetrization algorithm often improves the numerical stability and is enabled by default for TRAH-CASSCF and CASSCF linear-response calculations. However, the full-symmetrization algorithm may come with additional costs that depend on the number of \(\mathbf{F}\) intermediates. The number of \(\mathbf{F}\) intermediates depends on the symmetry of the density matrix (symmetric (S), anti-symmetric (A), and non-symmetric (N)) and whether overlap fitting is employed, as summarized in the table below.

Table 7.5 Number of COSX \(\mathbf{F}\) intermediates per density matrix

S fitting	P symmetry (S, A, N)	Number of \(\mathbf{F}\) intermediates
No	S / A	1
No	N	2
Yes	S / A	2
Yes	N	4

Note that for symmetric (S) and anti-symmetric (A) densities, we symmetrize and anti-symmetrize, respectively, exchange matrices at the end, which reduces the number of \(\mathbf{F}\) intermediates by a factor of 2. The actual computational costs usually do not increase linearly with the number of \(\mathbf{F}\) intermediates since we compute the costly analytic integrals (\(A_{\nu\tau}(\mathbf{r}_g)\)) only once and then contract them with the additional \(\mathbf{F}\) intermediates. From our experience the overhead of the full-symmetrization algorithm is roughly between 1.2 and 1.5 times that of the original algorithm.

7.4.2.7. COSX Grid and Convergence Issues

Symptoms of convergence issues: Erratic convergence behavior, often starting from the first SCF step or possibly setting in at a later stage, which produce crazy energy values with “megahartree” jumps. If overlap fitting is on, the following error message may also be encountered: “Error in Cholesky inversion of numerical overlap”.

Convergence issues may arise when the chosen grid has difficulties in representing the basis set. This is the “grid equivalent” of a linear dependence problem, as discussed in Linear Dependence. It should also be mentioned that the grid-related problem discussed here often goes hand in hand with basis set linear dependence, although not necessarily. The most straightforward way of dealing with this is to increase the size of the integration grid. This, however, is not always desirable or possible.

One way to avoid the Cholesky inversion issue is to turn overlap fitting off (!NoSFitting). However, this means that the numerical problems are still present, but are ignored. Due to the fact that overlap fitting operates with the numerical overlap and its inverse, it is more sensitive to linear dependence issues, so turning off the fitting procedure may lead to convergence. This may be a pragmatic — but by no means clean — solution, since it relies on the assumption that the numerical errors are small.

On the other hand, overlap fitting also gives a similar tool to deal with linear dependence issues as the one discussed for basis sets. The eigenvalues of the numerical overlap can be similarly inspected and small values screened out. There is unfortunately no universal way to determine this screening parameter, but see Linear Dependence for typical values.

The parameters influencing the method used for inversion and obtaining the fitting matrix are:

%method
  SFitInvertType Cholesky    # Cholesky inversion (default)
                 Cholesky_Q  # Cholesky + full Q matrix
                 Diag        # Inversion via diagonalization
                 Diag_Q      # Diag + full Q matrix
  SInvThresh     1e-8        # Inversion threshold for Diag and Diag_Q (default 1e-8)
end

By default, the inversion procedure proceeds through Cholesky decomposition as it is the fastest option. Ideally, the overlap matrix is non-singular if the basis set is not linearly dependent. For singular matrices, the Cholesky procedure will fail. It should be noted at this point that the numerical overlap can become linearly dependent even if the overlap of basis functions is not, and so a separate parameter will be needed to take care of grid-related issues. To achieve this, a diagonalization procedure (Diag) can be used instead of Cholesky with the corresponding parameter to screen out eigenvectors belonging to eigenvalues below a threshold (SInvThresh). For both Cholesky and diagonalization procedures, a “full Q” approach is also available (Cholesky_Q and Diag_Q), which corresponds to the use of a more accurate (untruncated) fitting matrix.

7.4.2.8. Treatment of Dispersion Interactions with DFT-D3

7.4.2.8.1. Introduction

DFT-D3 is an atom-pairwise (atom-triplewise) dispersion correction which can be added to the KS-DFT energies and gradient [324]:

(7.26)\[E_{\text{DFT-D3}} = E_{\text{KS-DFT}} + E_{\text{disp}} \]

\(E_{\text{disp}}\) is then the sum of the two- and three-body contributions to the dispersion energy, \(E_{\text{disp} } =E^{(2) }+E^{(3) }\). Most important is the two-body term, which is given at long range by:

(7.27)\[E_{\text{disp} } = - \sum_{A < B} \sum_{n=6,8} s_{n} \frac{C_{n}^{AB} }{r_{AB}^{n} } \]

where \(C_{n}^{AB}\) denotes the averaged (isotropic) \(n^{\text{th}}\)-order dispersion coefficient for atom pair AB and \(r_{AB}\) is their internuclear distance. \(s_{n}\) is a functional-dependent scaling factor (see below). In the general case, an adequate damping function must be employed.

7.4.2.8.2. Damping Functions

In order to avoid near-singularities for small \(r_{AB}\), the dispersion contribution needs to be damped at short distances. One possible way is to use rational damping as proposed by Becke and Johnson [90, 424, 425]:

(7.28)\[E^{(2) }= - \sum_{A < B} \sum_{n=6,8} s_{n} \frac{C_{n}^{AB} }{r_{AB}^{n} +f(R_{0}^{AB} )^{n} } \]

with [425]

(7.29)\[R_{0}^{AB} =\sqrt{ \frac{C_{8}^{AB} }{C_{6}^{AB} } } \]

and

(7.30)\[f(R_{0}^{AB} )=a_{1} R_{0}^{AB} +a_{2} . \]

Damping the dispersion contribution to zero for short ranges (as in Ref. [324]) is also possible:

(7.31)\[E^{(2) }= - \sum_{A < B} \sum_{n=6,8} s_{n} \frac{C_{n}^{AB} }{r_{AB}^{n} }f_{d,n} (r_{AB} ) \]

with

(7.32)\[f_{d,n} =\frac{1}{1+6(\frac{r_{AB} }{s_{r,n} R_{0}^{AB} })^{-\alpha_{n} } } \]

Note that the \(R_{0}^{AB}\) used with this damping are from Ref. [324]. For more information on the supported damping functions, see Ref [326]. The recommended variant is that with Becke-Johnson damping and is invoked by the keyword !D3BJ or simply !D3. The dispersion correction with zero damping is invoked by the keyword !D3ZERO.

7.4.2.8.3. Three-body term

It is possible to calculate the three-body dispersion contributions with DFT-D3, according to

(7.33)\[E^{(3) }= - \sum_{A < B < C} \frac{C_{9}^{ABC} (3\cos{\theta_{a} } \cos{\theta_{b} } \cos{\theta_{c} } +1) }{(r_{AB} r_{BC} r_{CA} )^{3} } f_{d,(3) }(\overline{r}_{ABC}) \]

where \(\theta_{a}\), \(\theta_{b}\) and \(\theta_{c}\) are the internal angles of the triangle formed by \(r_{AB}\), \(r_{BC}\) and \(r_{CA}\). The \(C_{9}\) coefficient is approximated by

(7.34)\[C_{9}^{ABC} \approx -\sqrt{ C_{6}^{AB} C_{6}^{AC} C_{6}^{BC} } \]

The three-body contribution has a small effect on medium-sized molecules and is damped according to equation (7.33). The damping function \(f_{d,(3) }(\overline{r}_{ABC})\) is similar to the one shown in equation (7.32) with \(\overline{r}_{ABC}\) being the geometric mean of \(r_{AB}\), \(r_{BC}\) and \(r_{CA}\). It can be used with both variants of the \(E^{(2) }\) term, although the three-body term itself will always be calculated using the zero damping scheme. Adding the three-body correction has proven to give quite accurate results for the thermochemistry of supramolecular systems[322]. The three-body term is invoked by the keyword !ABC, along with either the !D3ZERO keyword for zero damping of \(E^{(2) }\) or !D3BJ for Becke-Johnson damping of \(E^{(2) }\).

7.4.2.8.4. Options

Note that correcting Hartree-Fock (HF) is only parametrized with BJ-damping. For a constantly updated list of supported functionals, check the website of the Grimme group [329]. If there is a functional on this website that is parametrized, but the parameters are not yet implemented into ORCA, you can specify the parameters manually as shown below (using the respective parameters from [329]). In the same fashion, one could also use one’s own parameters, but this is not recommended.

Important: GGA and hybrid functionals should only be used with \(s_{6} =1.0\) to ensure asymptotically correct behavior. Double-hybrid functionals already account for parts of the dispersion interaction and should therefore be used with \(s_{6} < 1.0\). In the %method block, it is possible to change the \(s_6\), \(a_1\), \(s_8\), and \(a_2\) parameters for the variant with Becke-Johnson damping:

! d3bj b2plyp

%method
  D3S6 0.64
  D3A1 0.3065
  D3S8 0.9147
  D3A2 5.0570
end

The variant with zero damping offers the parameters \(s_6\), \(rs_6\), \(s_8\), and \(\alpha_{6}\).

! d3zero blyp

%method
  D3S6     1.0
  D3RS6    1.094
  D3S8     1.682
  D3alpha6 14
end

If a geometry optimization is performed (!OPT), then the program automatically calls the DFT-D3 gradient. There are also special functional parameters, which were optimized for triple-zeta basis sets. This option is only available with zero damping and can be invoked by the keywords !D3ZERO D3TZ. Preliminary results in the SI of Ref. [324] indicate that results are only slightly worse than with quadruple-zeta basis sets using the default parameters. This option should be carefully tested for future use in large computations.

7.4.2.8.5. Example input files

Following are some example input files. In this first example, a computation using the DFT-D3 dispersion correction with BJ-damping is run. The use of !D3BJ is identical to !D3 as the BJ-damping is on by default.

! pbe svp d3bj

* xyz 0 1
  C     0.000000     0.000000     0.000000
  O     0.000000     0.000000     1.400000
  O     0.000000     0.000000    -1.400000
*

The output for the dispersion correction in the ORCA output will look like this:

-------------------------------------------------------------------------------
                          DFT DISPERSION CORRECTION                            
                                                                               
                              DFTD3 V3.1  Rev 1                                
                          USING Becke-Johnson damping                          
-------------------------------------------------------------------------------
The PBE functional is recognized
Active option DFTDOPT                   ...         4   

molecular C6(AA) [au] = 156.562480 


            DFT-D V3
 parameters
 s6 scaling factor         :     1.0000 
 a1 scaling factor         :     0.4289 
 s8 scaling factor         :     0.7875 
 a2 scaling factor         :     4.4407 
 ad hoc parameters k1-k3   :    16.0000     1.3333    -4.0000 

 Edisp/kcal,au: -0.563071585638  -0.000897311593 
 E6   /kcal   :  -0.390909076 
 E8   /kcal   :  -0.172162510 
 % E8         :  30.575598941 

-------------------------   ----------------
Dispersion correction           -0.000897312
-------------------------   ----------------


-------------------------   --------------------
FINAL SINGLE POINT ENERGY      -188.136908447288
-------------------------   --------------------

\(E_{\text{disp} }\) is given as the “Dispersion correction” and is automatically included in the final single point energy. As discussed above, the parameters \(s_6\), \(a_1\), \(s_8\), and \(a_2\) may be manually defined by:

! pbe svp d3bj

%method 
  D3S6 1.0
  D3A1 0.4289
  D3S8 0.7875
  D3A2 4.4407
end

*xyz 0 1
  C     0.000000     0.000000     0.000000
  O     0.000000     0.000000     1.400000
  O     0.000000     0.000000    -1.400000
*

This results in the same output as above, but with additional messages that user inputs were found for the parameters:

A user input s6-coefficient scaling factor has been recognized
A user input a1-coefficient scaling factor has been recognized
A user input s8-coefficient scaling factor has been recognized
A user input a2-coefficient scaling factor has been recognized

The calculation of the same system using zero damping is run by the input:

! pbe svp d3zero

*xyz 0 1
  C     0.000000     0.000000     0.000000
  O     0.000000     0.000000     1.400000
  O     0.000000     0.000000    -1.400000
*

As stated above, the parameters \(s_6\), \(rs_6\), \(s_8\) and \(\alpha_{6}\) can be defined by the user. The next example shows this along with the call for the three-body correction (!ABC):

! pbe svp d3zero abc

%method 
  D3S6 1.0000
  D3RS6 1.2170
  D3S8 0.7220
  D3ALPHA6 14.0
end

*xyz 0 1
  C     0.000000     0.000000     0.000000
  O     0.000000     0.000000     1.400000
  O     0.000000     0.000000    -1.400000
*

7.4.2.9. DFT Calculations with the Non-Local, Density Dependent Dispersion Correction (VV10): DFT-NL

Accounting for the missing van der Waals (vdW, dispersion) forces in standard Kohn-Sham Density Functional Theory (DFT) has become essential in many studies of chemical and physical electronic structure problems. Common approaches use atom pair-wise additive schemes such as the popular DFT-D3 [324, 326] method, which is also available in ORCA by invoking the keyword !D3 (see section Treatment of Dispersion Interactions with DFT-D3), but these are not self-consistent, and thus do not correct the MOs or any computed property besides the energies.

A different route is followed by the van der Waals Density Functional (vdW-DF) as pioneered by Langreth and Lundquist [503]. These methods use only the electron density to include such dispersion/correlation effects. The vdW correlation functional VV10 of Vydrov and Van Voorhis [876] currently seems to be the most promising candidate for a general and accurate electronic structure method.

We use the term DFT-NL for any density functional in combination with the non-local part of the VV10 functional with an optimized parameter \(b\), which will be defined below. The performance of these methods has been evaluated in Ref. [400] using the GMTKN30 [307, 308, 309] database and the S66 set [6]. The performance of weak hydrogen bonds were evaluated in Ref. [401].

DFT-NL and DFT-D3 perform very similarly, but NL is to be preferred for metallic systems or when the basic electronic structure changes significantly (e.g. oxidations or ionizations). As a recent example, Janes and Iron showed that for functionals such as wB97X-V, including VV10 correlation results in very high quality reaction barriers when metals are involved [411].

The total exchange-correlation (XC) energy of VV10-type functionals is defined in eq. (7.35). It is composed of standard exchange (X) and correlation (C) parts and the non-local (NL) term, which covers (mainly) long-range dispersive energy:

(7.35)\[E_{\text{XC} }^{\text{DFT-NL} } = E_{\text{X} }^{\text{(hybrid)GGA} } +E_{\text{C} }^{\text{GGA} } +E_{\text{C-NL} }^{\text{VV10} } \]

The NL term is given by the following double integral:

(7.36)\[E_{\text{C-NL} }^{\text{VV10} } = \int dr\rho (r)\left[{ \beta +\frac{1}{2}\mathop \int\nolimits dr'\rho ( { r'} ) \varphi (r,r') } \right]\]

where \(\rho\) is the total electron density, and the definitions of the kernel \(\varphi (r,r')\) and \(\beta\) are as follows (in a.u.):

(7.37)\[\varphi \left({ r,r'} \right) = -\frac{3}{2gg'\left({ g+g'} \right)} \]

(7.38)\[\beta = \frac{1}{32}\left[{ \frac{3}{b^{2} }} \right]^{3/4} \]

with

\[\begin{split}\begin{aligned} g\left( r \right) &= \omega_{0} \left( r \right)R^{2}+\kappa (r) \\ R &= \left|{ r-r'} \right| \\ \omega_{0} \left( r \right) &= \sqrt{ C\left|{ \frac{\nabla \rho (r) }{\rho (r) }} \right|^{4}+\frac{4\pi }{3}\rho (r) } \\ \kappa \left( r \right) &= b\frac{3\pi }{2}\left[{ \frac{\rho (r) }{9\pi } } \right]^{1/6} \end{aligned} \end{split}\]

In the original definition, the short-range attenuation parameter \(b\) appearing in \(\kappa\) and \(\beta\) was fitted to the S22 set [429] of non-covalent interactions (\(b = 5.9\) for the rPW86PBE GGA). The other parameter \(C = 0.0093\), appearing in \(\omega_{0}\), determines the long-range behavior, and was set to its original value. Other DFT-NL functionals are constructed analogously. For a detailed discussion of the derivation of the formulas and their physical meaning and basis, see the references given above.

The defined energy of the non-local DFT-NL exchange-correlation functional can be computed in two ways: as a post-SCF addition based on a converged SCF density or in a self-consistent treatment. We take B3LYP as an example.

In our implementation of the non-self-consistent B3LYP-NL functional, a self-consistent B3LYP computation is performed as the first step. In the second step, the optimized electron density from the B3LYP computation is taken as input for the energy calculation of the non-local part. This procedure is invoked by the combination of the keywords !B3LYP NL. The use of the keywords !B3LYP SCNL would request a self-consistent treatment in which orbitals and density are optimized in the presence of the full B3LYP + VV10 exchange-correlation potential. According to many test calculations, an SCNL treatment is rarely necessary for normal energy evaluations.

The computation of the double integral given in eq. (7.36) requires using an integration grid, just like for normal exchange-correlation functionals. This grid is built similarly to the regular grids available in the ORCA, see Sec. Details on the numerical integration grids.

In the following, we compute the energy of the argon dimer at a distance of 3.76 Å with the def2-TZVP basis set and using the B3LYP hybrid functional as an example. Here, we choose the non-self-consistent variant of the DFT-NL dispersion correction.

! B3LYP NL
! def2-TZVP def2/JK RIJK

*xyz 0 1
  Ar  0.0 0.0 0.0
  Ar  0.0 0.0 3.76
*

The DFT-NL output for this example is shown below:

-------------------------------------------------------------------------------
                    post-SCF DFT-NL dispersion correction                      
-------------------------------------------------------------------------------

SCF Energy:    -1054.960547980 
NL    Energy:      0.209449625 
SC+NL Energy:  -1054.751098355 
NL done in  :              0.1 sec
-------------------------------------------------------------------------------

----------------
TOTAL SCF ENERGY
----------------

(...)

DFT components:
N(Alpha)           :       17.999996327923 electrons
N(Beta)            :       17.999996327923 electrons
N(Total)           :       35.999992655845 electrons
E(X)               :      -47.880990295917 Eh       
E(C)               :       -1.761924720763 Eh       
NL Energy, E(C,NL) :        0.209449624689 Eh       
E(XC)              :      -49.433465391991 Eh       

Here, we find the B3LYP total energy (”SCF Energy”), the non-local contribution (”NL Energy”), and their sum (”SC+NL Energy”), which is the final total energy. In the “DFT components” section, the non-local contribution is listed separately (”NL Energy, E(C,NL)”) in order to be consistent with the !SCNL output.

In the current version of ORCA, there are several functionals with pre-fitted \(b\) parameters (the parameter \(C\) is not changed), available for use by the keyword !DF NL or !DF SCNL, where DF stands for one of the following density functional keywords:

DF keyword
(meta-)GGA
BLYP
PBE
REVPBE
RPBE
RPW86PBE
SCANfunc
TPSS
hybrid
B3LYP
B3LYP/G
B3P86
B3PW91
mPW1PW
PBE0
PW1PW
PW6B95
REVPBE0
REVPBE38
TPSSh
TPSS0
R2SCANh (see [128])
R2SCAN0 (see [128])
R2SCAN50 (see [128])
range-separated hybrid
WR2SCAN50 (see [891])
double-hybrid
B2PLYP (see [51])
DSD-BLYP (see [906]). Same \(b\) used for DSD-BLYP/2013
DSD-PBEP86 (see [906]). Same \(b\) used for DSD-PBEP86/2013
PWPB95 (see [906])
PR2SCAN50 (see [891])
KPR2SCAN50 (see [891])
PR2SCAN69 (see [891])
range-separated double-hybrid
WPR2SCAN50 (see [891])

Additionally, one can include the non-local term in Hartree-Fock (HF) with a parameter of \(b = 3.9\) using the keywords “!HF NL”.

7.4.2.9.1. The B97-V Family

Head-Gordon’s \(\omega\)B97X-V functional [555] is a reparametrized version of the range-separated \(\omega\)B97X, which makes use of the non-local VV10 kernel to capture London-dispersion effects (\(b = 6.0\) and \(C = 0.01\); note that \(C\) is, unlike for the other functionals, changed for \(\omega\)B97X-V). The keyword !wB97X-V evaluates the VV10 kernel in a post-SCF way by default (i.e. only the semi-local exchange-correlation part is converged self-consistently and the resulting density is then used to assess the VV10-type energy contribution). A recent study showed that this saves computer time and does not have a significant effect on the resulting relative energies [603]. The keyword “!NL” does not have to be specified in this case as the VV10 kernel is invoked automatically. If a user wishes to carry out fully self-consistent calculations with \(\omega\)B97X-V, the “!SCNL” keyword must be specified.

The range-separated meta-GGA hybrid \(\omega\)B97M-V [557] and the meta-GGA B97M-V [556] are also available (via keywords !wB97M-V and !B97M-V, respectively). In the spirit of \(\omega\)B97X-V, the VV10 (NL) correction is called automatically in the post-SCF way by default.

7.4.2.9.2. Changing the \(b\), \(C\) and Scaling Parameters

All density functionals that are available in ORCA, but for which no \(b\) parameter is available, can be used with the DFT-NL method by providing a value for the parameter \(b\) as shown here:

%method
  NLb 5.0
end

The other parameter \(C\), appearing in \(\omega_{0}\), may also be changed, as shown in the following example.

! B3LYP NL 
! def2-TZVP def2/JK RIJK nososcf nopop

%method
  NLC    0.0083
end

*xyz 0 1
  Ar  0.0 0.0 0.0
  Ar  0.0 0.0 3.76
*

Of course, both parameters may be changed by using the NLb and NLC keywords in the %method block at the same time.

Note

In order to improve the scaling for larger systems, a distance cutoff was also introduced, controlled by the vdWdistTCUT flag in the %method block. The default value is 10 Å, so two grid points more than 10 Å away from each other do not correlate via the VV10 potential. This is already very conservative and has practically zero effect on the final energy, but can be altered if needed.

For methods where the long-range correlation is already partly covered, e.g. in double-hybrids, the NL energy can be scaled down using the NLScal parameter as shown below. By default, the scaling factor is \(1.0\) if not otherwise specified by the employed functional.

%method
  NLScal 0.5
end

7.4.2.9.3. Self-consistent Computations with the DFT-NL Dispersion Correction

Self-consistent calculations with the DFT-NL dispersion correction are possible by using the keyword !SCNL in combination with one of the available density functionals. Note that due to technical reasons, self-consistent calculations are not possible with the Hartree-Fock method.

For example, we can use the B3LYP hybrid functional with the self-consistent DFT-NL variant by the following input:

! B3LYP SCNL 
! def2-TZVP def2/JK RIJK

*xyz 0 1
  Ar  0.0 0.0 0.0
  Ar  0.0 0.0 3.76
*

The output is the same as a normal SCF calculation, but after convergence, a line with NL Energy, E(C,NL) is printed under “DFT components”, as it was for the post-SCF DFT-NL variant.

• Analytic gradients are already available, thus geometry optimizations with numerical frequencies can be computed using these functionals.

• TD-DFT calculations are not yet available.

• Any calculation that requires second derivatives of the NL functional is not yet possible. These are needed for real type perturbations in the CP-SCF solutions, e.g. for analytic Hessians, dipole polarizabilities, and double-hybrid gradients.

• Strictly imaginary perturbations such as NMR shielding and EPR g-tensors (both also with GIAOs), and hyperfine couplings are available.

7.4.2.10. DFT and HF Calculations with the Geometrical Counterpoise Correction: gCP

The central idea of the gCP correction [476] is to add a semi-empirical correction \(\Delta E_{\text{gCP} }\) to the energies of molecular systems that removes artificial overbinding effects from BSSE (see section Counterpoise Correction). gCP uses atomic corrections and therefore also has the ability to correct for intramolecular BSSE. The parametrization is constructed such that it approximates the Boys and Bernadi counterpoise (CP) correction \(\Delta E_{CP}\) in the intermolecular case. That is,

(7.39)\[\Delta E_{CP}\approx \Delta E_{\text{gCP} }\]

where e.g. for a complexation reaction \(A+B\to C\), our correction is given by

(7.40)\[\Delta E_{\text{gCP} }=E_{\text{gCP} }(C)-E_{\text{gCP} }(A)-E_{\text{gCP} }(B)\]

In practice, \(E_{\text{gCP} }\) can be simply added to the HF/DFT energy

(7.41)\[E_{\text{total}} = E_{\text{HF/DFT}} + E_{\text{gCP} }\]

which is done automatically in ORCA (the FINAL SINGLE POINT ENERGY includes the gCP correction).

The central equation of the gCP correction over all atoms \(N\) reads:

(7.42)\[E_{\text{gCP} } = \sigma \cdot \sum_{a}^{N} \sum_{b \neq a}^{N} e_a^{\text{miss} } \cdot f_{\text{dec} }(R_{ab})\]

where the energy \(e_a^{\text{miss} }\) is a measure of the incompleteness of the chosen target basis set (which is typically small) and \(f_{\text{dec} }(R_{ab})\) is a decay function that depends on the interatomic distance \(R_{ab}\). The scaling factor \(\sigma\) is one of 4 parameters needed for every method/basis set combination. More details on this can be found in the original gCP paper [476].

Analytical gradients with gCP are available for geometry optimization. Due to its semi-empirical nature, the correction is calculated within seconds, even for very large systems.

The gCP correction can be invoked by using the !GCP(level) keyword, where level is a compound of the method (HF or DFT) and the basis set keyword. See Table 7.6 for the available basis sets and the corresponding keywords. For example, gCP can be invoked in a B3LYP calculation with the def2-SV(P) basis set using the input:

! B3LYP def2-SV(P) GCP(DFT/SV(P))
*xyzfile 0 1 example.xyz

Table 7.6 Overview of parametrized basis sets. The level keyword in !GCP(level) is a compound of HF or DFT and the basis set keyword. Valid inputs are, for example: !GCP(HF/MINIS), !GCP(DFT/LANL), !GCP(HF/TZ), !GCP(DFT/631GD), \(\ldots\)

parametrized basis set	HF	DFT	basis set
MINIS	yes	yes	MINIS
SV	yes	yes	SV
6-31G(d)	yes	yes	631GD
6-31G(d) + LANL2DZ (Sc-Zn)	no	yes	LANL
def2-SV(P)	no	yes	SV(P)
def2-SV(P/h,c)	no	yes	SVX or SV(P/H,C)
def2-SVP	yes	yes	SVP
def2-TZVP	yes	yes	TZ

At all print levels, warnings from the gCP program are printed. Using the default print level, the only additional output is the final gCP correction before FINAL SINGLE POINT ENERGY. Using !LargePrint or %output Print[P_gCP] 2 end states the gCP level, the 4 parameters mentioned above, and the computed correction. A larger print level can be specified to get more details (more information on this below).

Several warnings and notices may be issued. Elements of the 5th and higher periods are treated as their 4th period analogs — e.g. Sn is treated with the same parameters as Ge. If this is the case, a note is printed. Another note is printed if there is a mismatch between the basis used for the SCF calculation and that of the requested gCP calculation. For example, the following input with tetramethyltin

! HF def2-SVP GCP(HF/MINIS)
*xyzfile 0 1 tetramethyltin.xyz

should use the parameters of Ge in place of Sn and there should be a mismatch between the basis set in ORCA (def2-SVP) and gCP (MINIS). Sure enough, the output is as follows:

 ** NOTE ** -> element sn will be treated as ge
NOTIFICATION: Different basis set in ORCA and otool_gcp:
ORCA: 142  gCP: 32 
If you are NOT using ECPs, check your basis set inputs!
------------------   -----------------
gCP correction             0.073031339
------------------   -----------------

A mismatch between the basis sets used is allowed since a minor mismatch may only result in a small error. One should still be careful with such results; use your own judgment! This also allows gCP in calculations that use an unparametrized basis set. However, in this case, the number of basis functions and exponents should be very similar!

It should be noted that some elements are not parametrized, depending on the gCP level used. If only a few atoms in a large molecule are treated inaccurately or not at all, the error is expected to be small. To check all parameters used and the individual atomic contributions, specify the print level %output Print[P_gCP] 3 end. For example, the above tetramethyltin input with this print level has the following output:

------------------------------------------------------------------------------
                                                                              
                  g C P  - geometrical counterpoise correction     
                                                                              
------------------------------------------------------------------------------
Method:  hf/minis 
  
 ** NOTE ** -> element sn will be treated as ge
Parameters:   sigma        eta         alpha        beta 
              0.1290      1.1526      1.1549       1.1763 
Nbf:   32  
NAtoms:   17  
  
gCP element parameters: 
  elem   emiss   nbas    elem   emiss   nbas    elem   emiss   nbas
    h    0.04240    1      he   0.02832    1      li   0.25266    2
    be   0.19720    2      b    0.22424    5      c    0.27995    5
    n    0.35791    5      o    0.47901    5      f    0.63852    5
    ne   0.83235    5      na   1.23292    6      mg   1.34339    6
    al   1.44828    9      si   1.61336    9      p    1.76814    9
    s    1.99201    9      cl   2.23311    9      ar   2.49323    9
    k    3.02964   10      ca   3.38998   10      sc   0.00000    0
    ti   0.00000    0      v    0.00000    0      cr   0.00000    0
    mn   0.00000    0      fe   0.00000    0      co   0.00000    0
    ni   0.00000    0      cu   0.00000    0      zn   0.00000    0
    ga   0.00000    0      ge   0.00000    0      as   0.00000    0
    se   0.00000    0      br   0.00000    0      kr   0.00000    0
  
      #     ON  sites    Nvirt           Emiss    BSSE [kcal/mol]
      1      6     15      2.0          0.2799         10.0090
      2     32     16    -16.0          0.0000          0.0000
      3      6     15      2.0          0.2799         10.0084
      4      6     15      2.0          0.2799         10.0095
      5      6     15      2.0          0.2799         10.0083
      6      1     15      0.5          0.0424          0.4827
      7      1     15      0.5          0.0424          0.4828
      8      1     15      0.5          0.0424          0.4828
      9      1     15      0.5          0.0424          0.4827
     10      1     15      0.5          0.0424          0.4827
     11      1     15      0.5          0.0424          0.4827
     12      1     15      0.5          0.0424          0.4828
     13      1     15      0.5          0.0424          0.4827
     14      1     15      0.5          0.0424          0.4828
     15      1     15      0.5          0.0424          0.4827
     16      1     15      0.5          0.0424          0.4827
     17      1     15      0.5          0.0424          0.4827
Egcp:       0.0730313386 a.u.
NOTIFICATION: Different basis set in ORCA and otool_gcp:
ORCA: 142  gCP: 32 
If you are NOT using ECPs, check your basis set inputs!
------------------   -----------------
gCP correction             0.073031339
------------------   -----------------

From this, it can be seen that the Sn atom (atom 2 in the list of atomic contributions) gives no contribution because its Emiss is zero (unparametrized for the given gCP level). This is confirmed by looking at the gCP element parameters section, which lists the emiss of Sn as zero for !GCP(HF/MINIS). Rerunning this example with !GCP(HF/SVP) now gives a contribution for Sn, as seen by the following output. Note that this calculation also has a much smaller basis set mismatch, and so should be the more accurate gCP correction of the two.

------------------------------------------------------------------------------
                                                                              
                  g C P  - geometrical counterpoise correction     
                                                                              
------------------------------------------------------------------------------
Method:  hf/svp 
  
 ** NOTE ** -> element sn will be treated as ge
Parameters:   sigma        eta         alpha        beta 
              0.2054      1.3157      0.8136       1.2572 
Nbf:   148  
NAtoms:   17  
  
gCP element parameters: 
  elem   emiss   nbas    elem   emiss   nbas    elem   emiss   nbas
    h    0.00811    5      he   0.00805    5      li   0.11358    9
    be   0.02837    9      b    0.04937   14      c    0.05538   14
    n    0.07278   14      o    0.10031   14      f    0.13327   14
    ne   0.17360   14      na   0.18114   15      mg   0.12556   18
    al   0.16719   18      si   0.14984   18      p    0.14540   18
    s    0.16431   18      cl   0.18299   18      ar   0.20567   18
    k    0.20096   24      ca   0.29966   24      sc   0.32599   31
    ti   0.30549   31      v    0.29172   31      cr   0.29380   31
    mn   0.29179   31      fe   0.29673   31      co   0.30460   31
    ni   0.24204   31      cu   0.35419   31      zn   0.35072   31
    ga   0.35002   32      ge   0.34578   32      as   0.34953   32
    se   0.36731   32      br   0.38201   32      kr   0.39971   32
  
      #     ON  sites    Nvirt           Emiss    BSSE [kcal/mol]
      1      6     16     11.0          0.0554          2.5093
      2     32     16     16.0          0.3458          6.8274
      3      6     16     11.0          0.0554          2.5092
      4      6     16     11.0          0.0554          2.5094
      5      6     16     11.0          0.0554          2.5092
      6      1     16      4.5          0.0081          0.1703
      7      1     16      4.5          0.0081          0.1703
      8      1     16      4.5          0.0081          0.1703
      9      1     16      4.5          0.0081          0.1703
     10      1     16      4.5          0.0081          0.1703
     11      1     16      4.5          0.0081          0.1703
     12      1     16      4.5          0.0081          0.1703
     13      1     16      4.5          0.0081          0.1703
     14      1     16      4.5          0.0081          0.1703
     15      1     16      4.5          0.0081          0.1703
     16      1     16      4.5          0.0081          0.1703
     17      1     16      4.5          0.0081          0.1703
Egcp:       0.0301322002 a.u.
NOTIFICATION: Different basis set in ORCA and otool_gcp:
ORCA: 142  gCP: 148 
If you are NOT using ECPs, check your basis set inputs!
------------------   -----------------
gCP correction             0.030132200
------------------   -----------------

The gCP input can also be defined in the %method block of the input:

%method
  DoGCP       true/false    # turn gCP on/off
  GCPMETHOD   "method"      # define method string for otool_gcp, e.g. "dft/svp"
  GCP.D3MINIS true/false    # use special DFT-D3 refit for HF/MINIS (default=true)
end

General advice:

• Small basis sets show not only a large BSSE, but general shortcomings. These effects are not always clearly distinguishable.

• If computationally affordable, large basis sets (triple-\(\zeta\) and higher) are always preferable for a given system.

• gCP only makes sense for somewhat large molecules

• gCP should always be applied together with a dispersion correction for DFT: gCP-D3 is well tested, but gCP-NL also works well (see sections Treatment of Dispersion Interactions with DFT-D3 and DFT Calculations with the Non-Local, Density Dependent Dispersion Correction (VV10): DFT-NL)

Note

• !GCP(HF/MINIS) automatically sets the refitted D3 parameter, as proposed in the original gCP paper.

• The gCP method is implemented via an external tool called otool_gcp, which is based on the original Fortran program used in the publication. Thus, the otool_gcp binary can also be called directly via the command line (otool_gcp -h gives an overview of the options).

• It is also possible to read an external parameter file ($HOME/.gcppar) using the !GCP(FILE) keyword. For further information, please look at the manual for the gcp program as provided by Prof. S. Grimme[8].

• During the calculation, some temporary output files are written by ORCA: BASENAME.gcp.in.tmp and BASENAME.gcp.out.tmp will contain the input for otool_gcp and its output, respectively.

• It has been demonstrated that the popular combination of B3LYP with 6-31G(d) can be strongly improved using DFT-D3 and gCP [475]. For convenience, the keyword !B3LYP-gCP-D3/6-31G* has been defined. This is equivalent to !B3LYP 6-31G* D3BJ GCP(DFT/631GD).

7.4.2.11. HF-3c: Hartree-Fock with three corrections

HF-3c is a fast Hartree-Fock based method developed for computation of structures, vibrational frequencies and non-covalent interaction energies in huge molecular systems [834]. The starting point for evaluating the electronic energy is a standard HF calculation with a small Gaussian AO basis set. The used so-called MINIX basis set consists of different sets of basis functions for different groups of atoms as shown in table Table 7.7.

Table 7.7 Composition of the MINIX basis set.

element	basis
H-He, B-Ne	MINIS
Li-Be	MINIS+1(p)
Na-Mg	MINIS+1(p)
Al-Ar	MINIS+1(d)
K-Zn	SV
Ga-Kr	SVP
Rb-Rn	def2-SVP with Stuttgart-Dresden ECPs

Three terms are added to correct the HF energy \(E_{\text{tot} }^{\text{HF/MINIX} }\) in order to include London dispersion interactions, to account for the BSSE and to correct for basis set deficiencies, i.e. overestimated bond lengths. The corrected total energy is therefore calculated as

(7.43)\[ E_{\text{tot} }^{\text{HF-3c} } = E_{\text{tot} }^{\text{HF/MINIX} } + E_{\text{disp} }^{\text{D3(BJ) }} + E_{\text{BSSE} }^{\text{gCP} } + E_{\text{SRB} }. \]

The first correction term \(E_{\text{disp} }^{\text{D3(BJ) }}\) is the atom-pair wise London dispersion energy from the D3 dispersion correction scheme[324] applying Becke-Johnson (BJ) damping [90, 424, 425] (see section Treatment of Dispersion Interactions with DFT-D3). The second term \(E_{\text{BSSE} }^{\text{gCP} }\) denotes the geometrical counterpoise (gCP) correction [476] to treat the BSSE (see section DFT and HF Calculations with the Geometrical Counterpoise Correction: gCP). For the HF-3c method, the three usual D3 parameters \(s_8\), \(a_1\) and \(a_2\) were re-fitted using reference interaction energies of the complexes of the S66 test set [6]. This results in \(s_8=0.8777\), \(a_1=0.4171\) and \(a_2=2.9149\). The parameter \(s_6\) was set to unity as usual to enforce the correct asymptotic limit and the gCP correction was already applied in this fitting step.

The last term \(E_{\text{SRB} }\) is a short-ranged correction to deal with basis set deficiencies which occur when using small or minimal basis sets. It corrects for systematically overestimated covalent bond lengths for electronegative elements and is calculated as a sum over all atom pairs:

\[E_{\text{SRB} } = -s \sum_{A}^{\text{Atoms} } \sum_{B \neq A}^{\text{Atoms} }(Z_A Z_B)^{3/2}\exp(-\gamma (R_{AB}^{0,\text{D3} })^{3/4} R_{AB})\]

Here, \(R_{AB}^{0,\text{D3} }\) are the default cut-off radii as determined ab initio for the D3 scheme [324] and \(Z_A\), \(Z_B\) are the nuclear charges. This correction is applied for all elements up to argon. The empirical fitting parameters \(s=0.03\) and \(\gamma= 0.7\) were determined to produce vanishing HF-3c total atomic forces for B3LYP-D3(BJ)/def2-TZVPP equilibrium structures of 107 small organic molecules. More details can be found in the original publication [834].

The HF-3c method can only be invoked with a simple keyword:

! HF-3c

! HF-3c is a compound keyword and equals ! HF MINIX D3BJ GCP(HF/MINIX) PATOM, hence no basis set etc. has to be specified. The PATOM guess is selected since the grid construction for the default guess can take more time than an actual SCF step. The guess can only be overwritten manually in the %method section.

The default mode for the integral handling is set to Conventional. The storing of the two-electron integrals on disk or in memory if possible leads to large computational savings. In case one want to use the Direct mode, this has to be specified in the %scf input section:

%scf
  SCFmode Direct
end

The output gives the used parameters and the correction itself for D3 and gCP separately. As the SRB correction is also calculated with the otool_gcp, the results are given in the gCP output section. The Total correction to HF/MINIX is the sum of all three corrections (D3, gCP and SRB) and FINAL SINGLE POINT ENERGY is the total HF-3c energy as given in equation (7.43).

-------------------------------------------------------------------------------
                          DFT DISPERSION CORRECTION

                                DFTD3 V2.1  Rev 6
                          USING Becke-Johnson damping
-------------------------------------------------------------------------------
The default Hartree-Fock  is recognized
Active option DFTDOPT                   ...         4

molecular C6(AA) [au] = 1689.256597


            DFT-D V3
 parameters
 using HF/MINIX parameters
 s6 scaling factor         :     1.0000
 a1 scaling factor         :     0.4171
 s8 scaling factor         :     0.8777
 a2 scaling factor         :     2.9149
 ad hoc parameters k1-k3   :    16.0000     1.3333    -4.0000

 Edisp/kcal,au: -32.163184627631  -0.051255291794
 E6   /kcal   : -18.007221978
 E8   /kcal   : -14.155962649
 % E8         :  44.012938437

-------------------------   ----------------
Dispersion correction            -0.051255292
-------------------------   ----------------


------------------------------------------------------------------------------

                  g C P  - geometrical counterpoise correction

------------------------------------------------------------------------------
Method:  hf/minix

Parameters:   sigma        eta         alpha        beta
              0.1290      1.1526      1.1549       1.1763
Egcp:       0.0723150636 a.u.
Ebas:      -0.0636976872 a.u.
------------------   -----------------
gCP+bas correction         0.008617376
------------------   -----------------
----------------------------   ----------------
Total correction to HF/MINIX       -0.042637915
----------------------------   ----------------


-------------------------   --------------------
FINAL SINGLE POINT ENERGY      -163.002895262171
-------------------------   --------------------

For the elements up to Xe, the default initial guess is a Hückel guess. Beyond Xe, the guess mode is changed to HCORE since no Hückel parameters for the respective ECP bases are available and other models are not implemented at the moment. For calculations with only lighter elements and therefore no ECPs, the ECP printouts in the output file can be ignored.

7.4.2.12. PBEh-3c: A PBE hybrid density functional with small AO basis set and two corrections

PBEh-3c is a highly efficient electronic structure approach performing particularly well in the optimization of geometries and for interaction energies of non-covalent complexes.[325] Here, a global hybrid variant of the Perdew-Burke-Ernzerhof (PBE) functional with a relatively large amount of non-local Fock-exchange (42%) is employed with a valence-double-zeta Gaussian AO basis set (def2-mSVP). Basis set superposition errors (BSSE) and London dispersion effects are accounted for by the gCP and D3 schemes, respectively (see above). The basis set is constructed such that:

Table 7.8 Composition of the def2-mSVP basis set.

element	basis
H	def2-SV(P) (\(\alpha\) scaled by 1.2)
He	def2-SVP(-p)
Li-Be,Na-Kr	def2-SV(P)
B,Ne	Ahlrichs’ DZ + Polarization from def2-SVP
C-F	Ahlrichs’ DZ + Polarization from 6-31G*
Rb-Rn	def2-SVP with Stuttgart-Dresden ECPs

For inter- and intramolecular BSSE the gCP expression from Eq. (7.42) is used but with a damping function (similar to the zero-damping in Eq. (7.32)). This damping improves the thermochemistry of the method significantly compared with the non-damped version. London dispersion effects are accounted for by the DFT-D3 (BJ-damping) scheme including the three-body term. Compared to the related HF-3c approach, the PBEh-3c is somewhat more costly, however, it yields much better geometries. These are roughly of MP2-quality (or even better for non-covalent structures) but may be computed at much lower cost. Due to the moderate amount of non-local Fock exchange, the method is less prone to self-interaction errors (as in GGAs) but still applicable in cases when Hartree-Fock fails (strongly correlated systems).

The PBEh-3c method may be invoked with the simple keyword:

! PBEh-3c

Identical to HF-3c, the default initial guess for all elements up to Xe is a Hückel guess. Beyond Xe, the guess mode is changed to HCORE. For calculations with only lighter elements and therefore no ECPs, the ECP printouts in the output file can be ignored.

Recently, a new composite ‘low-cost’ method for accurate thermochemistry, structures, and noncovalent interactions specifically also for transition metal chemistry and other stronger correlated systems was implemented. As it is based on the B97 GGA including D3 with three-body contribution, a short range bond length correction, and a modified, stripped-down triple-\(\zeta\) basis, def2-mTZVP, the computational cost of this method termed B97-3c are between that of HF-3c and PBEh-3c (for large systems roughly two times more expensive than HF-3c). It is invoked with a simple keyword analogously to the latter methods. Some more detailed information can be found in Ref. [121] .

7.4.2.13. \(r^2\)SCAN-3c: A robust “Swiss army knife” composite electronic-structure method

The \(r^2\)SCAN-3c composite method[333] is available as robust “Swiss army knife” electronic structure method for thermochemistry, geometries and non-covalent interactions and has shown in preliminary tests consistent performance for both open and closed shell transition metal complexes. It is based on the \(r^2\)SCAN[282] meta-GGA combined with the D4 dispersion correction[244] and the geometrical counter poise-correction[476]. The modified triple-\(\zeta\) basis set, def2-mTZVPP, is larger and more consistent for the light main-group elements and almost as computationally efficient as the def2-mTZVP basis set of B97-3c. The computational cost of \(r^2\)SCAN-3c is slightly larger than B97-3c. It is invoked with the simple keyword

! r2SCAN-3c

7.4.2.14. \(\omega\)B97X-3c: A composite range-separated hybrid DFT method with a molecule-optimized polarized valence double-\(\zeta\) basis set

The \(\omega\)B97X-3c composite method[540] is based on the \(\omega\)B97X-V functional and combines a tailored and molecule-optimized polarized valence double-\(\zeta\) (vDZP) basis set and a specifically adapted D4 dispersion correction. The vDZP basis set employs large-core ECPs and shows only very small basis set superposition and incompleteness errors compared to conventional double-\(\zeta\) basis sets. In thorough tests on standard benchmarks sets, the \(\omega\)B97X-3c method was shown to be on par with well-performing hybrid DFT methods in a quadruple-\(\zeta\) basis set at a fraction of their computational cost. \(\omega\)B97X-3c is consistently available for all elements up to Rn (Z = 1–86).

It is invoked with the simple keyword:

! wB97X-3c

The vDZP basis set alone is utilized as follows (note that the corresponding large-core ECPs are called automatically):

! vDZP

7.4.3. Semiempirical Methods

The present version of ORCA has inherited the capability of doing semiempirical calculations from the earlier versions. A number of methods based on the “neglect of differential overlap” [692, 779] are currently implemented for energies and analytic gradients (for geometry optimization).

%method
  Method  CNDO
          INDO
          NDDO
  # for Method=CNDO
  Version CNDO_1
          CNDO_2
          CNDO_S
  # for Method=INDO
  Version INDO_1
          INDO_2
          INDO_S
          ZINDO_1
          ZINDO_2
          ZINDO_S
  # for Method=NDDO
  Version ZNDDO_1
          ZNDDO_2
          MNDO
          AM1
          PM3
end

The methods MNDO [206, 207, 848], AM1 [208] and PM3 [822] are available for main group elements only and arise from the work of the Dewar group. They have been optimized to reproduce molecular structure and energetics. The older CNDO/1,2 and INDO/1,2 were developed by the Pople group [61, 175, 176, 177, 693, 696, 697, 749, 750] and were designed to roughly mimic minimal basis ab initio calculations. The methods of the Zerner group (ZINDO/1,2 and ZINDO/S) are closely related to the older methods but have been well parameterized for transition metals too [36, 37, 38, 63, 183, 466, 720, 908, 909, 910, 912]. ZINDO/1 (and less so ZINDO/2) are suitable for geometry optimization. ZINDO/S gives good results for electronically excited states at moderate configuration interaction levels and is also successful for the calculation of electron and spin distributions in large transition metal complexes [37, 183, 466, 908, 909, 910]. The ZNDDO/1,2 methods have been implemented into ORCA as straightforward extensions of the corresponding INDO methods without changing any parameter. However, the methods benefit from the somewhat more accurate representation of the Coulomb interaction within the NDDO approximation [431, 634]. The preliminary experience with these methods is that they are better than the corresponding INDO methods for calculation of transition metal complex structures but on the whole have also similar deficiencies.

The analytic gradients are available for all of these methods and can be used to produce reasonable molecular structures at low computational cost or to get preliminary insight in the behavior of the system under investigation[9].

There is also a mechanism for simplified input. Instead of giving values for Method and Version separately you can also assign the value that would normally belong to Method to Version. The program will recognize that and assign the correct values to both Method and Version.

%method
  # shortcut to Method=NDDO; and Version=AM1;
  Method  AM1
end

• If you want you can also combine semiempirical methods with MP2 (energies only). For example use Method \(=\)AM1; and DoMP2 \(=\)true; It is questionable if this makes the results of semiempirical calculations any better but at least it is possible in ORCA.

You can change the built-in semiempirical parameters in a straightforward fashion. For example:

! ZINDO/S TightSCF DIIS NoMOPrint

%cis    NRoots 20
        MaxDim 3 # Davidson expansion space = MaxDim * NRoots
        end

%ndoparas  P[6,25] 20
           P[6,26] 20
           end

The %ndoparas block is there in order to let you input your favorite personal parameters. The “molecular” parameters are set using “INTFA” (“interaction factors”);

%ndoparas  INTFA[PP_PI] 0.585
  # The interaction factors exist for
  #  ss_sigma
  #  sp_sigma
  #  sd_sigma
  #  pp_sigma
  #  pd_sigma
  #  dd_sigma
  #  pp_pi
  #  pd_pi
  #  dd_pi
  #  dd_delta
  # the parameter entering the Coulomb integrals
  # in INDO/S
  FGAMMA 1.2
  end

All atomic parameters are collected in an array “P”. The first index is the atomic number of the element whose parameters you want to change. The second index identifies which parameter. The list of parameters follows below. Most of them will only be interesting for expert users. The most commonly modified parameters are the Beta’s (number 25 through 28). Note that most programs require a negative number here. In ORCA the resonance integrals are defined in a way that makes the Beta’s positive.

# core integrals (in eV)
US      0
UP      1
UD      2
UF      3
# Basis set parameters (double-zeta for generality)
NSH     4  # number of shells for the element
NZS     5  # number of Slater type orbitals for the s shell
ZS1     6  # first exponent
ZS2     7  # second exponent
CS1     8  # first contraction coefficient
CS2     9  # second contraction coefficient
NZP    10  # number of Slater type orbitals for the p shell
ZP1    11  # ...
ZP2    12
CP1    13
CP2    14
NZD    15  # number of Slater type orbitals for the d shell
ZD1    16  # ...
ZD2    17
CD1    18
CD2    19
NZF    20  # number of Slater type orbitals for the f shell
ZF1    21  # ...
ZF2    22
CF1    23
CF2    24
# Resonance integral parameters (in eV)
BS     25  # s shell beta
BP     26  # p shell beta
BD     27  # d shell beta
BF     28  # f shell beta
# Number of electrons in the g.s.
NEL    29  # total number of electrons (integer)
NS     30  # fractional occupation number of the s shell
NP     31  # fractional occupation number of the p shell
ND     32  # fractional occupation number of the d shell
NF     33  # fractional occupation number of the f shell
# The one center repulsion (gamma) integrals (in eV)
GSS    34
GSP    35
GSD    36
GSF    37
GPP    38
GPD    39
GPF    40
GDD    41
GDF    42
GFF    43
# The Slater Condon parameters (in eV)
F2PP   44
F2PD   45
F2DD   46
F4DD   47
G1SP   48
G1PD   49
G2SD   50
G3PD   51
R1SPPD 52
R2SDPP 53
R2SDDD 54
# The nuclear repulsion parameters for Dewar type models
NR1    55
NR2    56
NR3    57
NR4    58
NR5    59
NR6    60
NR7    61
NR8    62
NR9    63
NR10   64
NR11   65
NR12   66
NR13   67
# The nuclear attraction/repulsion parameter for MNDO/d
RHO    68
# Spin orbit coupling parameters
SOCP   69 # SOC for the p shell
SOCD   70 # SOC for the d shell
SOCF   71 # SOC for the f shell

7.4.3.1. Semi-empirical tight-binding methods: Grimme’s GFN0-xTB, GFN-xTB and GFN2-xTB

ORCA is interfaced to the XTB tool by Grimme and coworkers, allowing the user to request all kinds of calculations using the popular GFN0-xTB, GFN-xTB and GFN2-xTB Hamiltonians.[70, 332, 699] From the technical side, the user has to provide the executable provided by the Grimme group. The xtb program package can be obtained free of charge from https://github.com/grimme-lab/xtb/releases and detailed information on the usage of the xtb standalone program and its utilities can be found at https://xtb-docs.readthedocs.io/en/latest/contents.html. Only the file bin/xtb is used by ORCA . The user should copy this file into the directory where the other ORCA executables are located, and rename it as otool_xtb.

Please use the 6.7.1 version (or any later version) of xtb; older versions are not fully compatible with ORCA anymore or are missing features, for example it may not be possible to invoke the solvation model! Additionally, Windows users should copy libiomp5md.dll from the XTB directory to the ORCA directory.

XTB is invoked by the following keywords:

! XTB0  # for GFN0-xTB. Synonym: GFN0-XTB
! XTB1  # for GFN1-xTB. Synonym: GFN-XTB
! XTB2  # for GFN2-xTB. Synonym: GFN2-XTB
! XTBFF # for GFN-FF. Synonym: GFN-FF

The following methods can be used in conjunction with XTB:

• Single Point Energy

• Energy and Gradient

• Optimization, using all kinds of constraints, relaxed surface scans, etc.

• Nudged-Elastic Band calculations

• Numerical Frequency Calculations

• Intrinsic Reaction Coordinate

• Molecular Dynamics Calculations

• QM/MM calculations

Note

• XTB0 is a non-self-consistent tight-binding method, and as such, its accuracy is generally inferior to XTB1 and XTB2 (and sometimes even XTBFF), despite that it is a few times faster than XTB1 and XTB2. From our experience, we only recommend XTB0 when both XTB1 and XTB2 exhibit qualitative failures for the system of interest.

• Please note that XTB0, XTB1 and XTB2 can also be used for the initial path generation or for the calculation of an initial TS structure on XTB level, both as input for the subsequent NEB calculation on a higher level of theory. For more details, please consult section Nudged Elastic Band Method.

7.4.3.1.1. Solvation

Three implicit solvation models can be requested in XTB calculations: (1) the analytical linearized Poisson-Boltzmann (ALPB) solvation model, (2) the domain decomposition COSMO (ddCOSMO)[136], and (3) the extended conductor-like polarizable continuum model (CPCM-X).[815] These three models can be requested via the following tags in the simple input

! ALPB(solvent) # use ALPB

or

! DDCOSMO(solvent) # use ddCOSMO

or

! CPCMX(solvent) # use CPCMX

where solvent is any of the solvents in Table 7.27.

7.4.3.1.2. Keywords for XTB

A list of additional keywords for XTB is detailed here:

! XTB
! ALPB(water)       # use ALPB for implicit solvation, solvent water, 
! DDCOSMO(water)    # use ddCOSMO for implicit solvation, solvent water,
! CPCMX(water)      # use CPCMX for implicit solvation, solvent water,
                    # they can also be defined in the xtb block
%xtb
  XTBINPUTSTRING  "argument 1" # string passed on to XTB call
  XTBINPUTSTRING2 "keyword 2"  # string passed on to XTB call
  ETemp           0 # electronic temperature, value for --etemp
  DOALPB      false # use implicit solvation, ALPB
  ALPBSOLVENT    "" # ALPB solvent, no default, string for --alpb
  DODDCOSMO   false # use implicit solvation, ddCOSMO
  DDCOSMOSOLVENT "" # ddCOSMO solvent, no default, string for --cosmo
  DOCPCMX     false # use implicit solvation, CPCMX
  CPCMXSOLVENT ""   # CPCMX solvent, no default, string for --cpcmx
  EPSILON       3.5 # Dielectric constant (only for ddCOSMO)       
  ACCURACY        1 # accuracy, value for --acc, default is ORCA's accuracy x 1.e6
  MAXCORE      1024 # memory in MB reserved for XTB calculation, 
                    # default is ORCA's maxcore
  NPROCS          1 # number of processors for running XTB", 
                    # default is ORCA's PAL command  
end

Note

• If jobs are run over several nodes, the number of cores used by the XTB tool might be lower than requested via the pal keyword.

---

[1]

There is some evidence that the version used in the Gaussian program gives slightly better results in molecular applications than the TurboMole variant, but the differences are very small [387].

[2]

https://libxc.gitlab.io

[3]

The basic theory behind the RI method has been known for a long time and since at least the late sixties, methods similar to the RI approximation have been used — mainly in the context of “approximate ab initio methods” such as LEDO, PDDO, and MADO, but also in density functional theory in the mid and late seventies by Baerends, Dunlap, and others [65, 225, 862, 888].

[4]

It probably should be noted that a slightly awkward step in the procedure is the inversion of the auxiliary integral matrix \(\mathrm{\mathbf{V}}\), which can easily become very large. Matrix inversion is an O(N\(^{3})\) process such that for large molecules, this step takes some real time. However, in ORCA, this is only done once during the calculation, whereas other programs that constrain the fit to also exactly reproduce the number of electrons perform a similar process each iteration. Starting from ORCA 2.2.09, the Cholesky decomposition is used in favor of matrix inversion, removing any bottleneck concerning the solution of the linear equation system.

[5]

This algorithm was described by Kendall and Früchtl [440].

[6]

The theory of this approach together with all evaluations and implementation details is described in [383, 624]. References to earlier work can also be found there.

[7]

The theory is not yet published, but will be soon.

[8]

http://www.thch.uni-bonn.de/tc/

[9]

However, do not try to use ZINDO/S (or CNDO/S) for structure optimizations - it does not make sense and will lead to disastrous results because there is no accurate representation of nuclear repulsion in these methods.