# The density functional program

Density-functional theory calculations may be performed using one of the following commands:

calculate functional of a previously computed density.`DFT`

calls the spin-restricted Kohn-Sham program.`RKS`

or`RKS-SCF`

`KS`

and`KS-SCF`

are aliases for`RKS`

.calls the spin-unrestricted Kohn-Sham program`UKS`

or`UKS-SCF`

Each of these commands may be qualified with the key-names of the functional(s) which are to be used, and further options:

*command, key1, key2, key3, $\ldots$*, *options*

If no functional keyname is given, the default is `LDA`

(see below). Following this command may appear directives specifying options for the density-functional modules (see section directives) or the Hartree-Fock program (see section options).

On completion of the functional evaluation, or self-consistent Kohn-Sham calculation, the values of the individual functionals are stored in the Molpro vector variable `DFTFUNS`

; the total is in `DFTFUN`

, and the corresponding individual functional names in `DFTNAME`

.

Energy gradients are available for self-consistent Kohn-Sham calculations.

Density fitting can be used for closed and open-shell spin-restricted KS and is invoked by a prefix `DF-`

(`DF-KS`

or `DF-RKS`

, see section density fitting). For UKS, only Coulomb fitting is possible (`CF-UKS`

). Density fitting is highly recommended (unless explicit symmetry handling is required), as the induced errors are negligible and it offers *massive* speed increases, particularly for pure functionals. For details see R. Polly, H.-J. Werner, F. R. Manby, and Peter J. Knowles, *Fast Hartree-Fock theory using local density fitting approximations*, Mol. Phys. **102**, 2311 (2004). **All publications resulting from DF-HF or DF-KS calculations should cite this work.**

Normally, sensible defaults are used to define the integration grid. The accuracy can be controlled using options as described in section options or directives as described in section directives. More control is provided by the `GRID`

command, as described in section numerical integration grid control (GRID).

## Options

The following options may be specified on the `KS`

or `UKS`

command lines:

Specifies the grid target accuracy (per atom). The default is 1.d-6 unless this has been modified using a global`GRIDTHR`

=*target*`THRESH, GRID`

option.If`COARSEGRID`

(logical)*true*, perform initial iterations with a coarser grid. Default is*false*.In the initial iterations, the grid accuracy is min(`GRIDMAX`

=*gridmax**gridmax, target*coarsefac*) (only if COARSEGRID is set).Factor for initial grid accuracy (see above). The default is 1000.`COARSEFAC`

=*coarsefac*Defines whether grid weight derivatives are included in analytical gradient calculations (default: 0). Disabling these can improve convergence in geometry optimizations.`GRIDGRAD`

=*0 or 1*Factors for each functional. The number of given values must agree with the number of functionals.`DFTFAC`

=`[`

*fac1,fac2,..*`]`

Fraction of exact exchange added to the functional. The default depends on the functional.`EXFAC`

=*factor*Threshold for orbital screening (default depends on energy threshold).`TOLORB`

=*value*

In addition, all options valid for HF (see section options) can be given.

## Directives

The following options may be used to control the operation of the DFT modules. In the Kohn-Sham case, these may come in any order before or after directives for the SCF program as described in Section the SCF program.

### Density source (DENSITY, ODENSITY)

`DENSITY`

,*orbc.filec*,$\dots$ `ODENSITY`

,*orbo.fileo*,$\dots$

For non-self-consistent `DFT`

calculations, specifies the source of the density matrix. The total density is read from *orbc.filec*, with further options specifying density sets in the standard way as described in Section selecting orbitals and density matrices (ORBITAL, DENSITY). `ODENSITY`

can be used to specify the spin density. The defaults are the densities last written by an SCF or MCSCF program.

### Thresholds (DFTTHRESH)

`DFTTHRESH`

,*key1=value1,key2=value2*$\ldots$

Sets various truncation thresholds. *key* can be one of the following.

Overall target accuracy (per atom) of density functional. Defaults to the value of the global thresholdcommand:gthresh`TOTAL`

`GRID`

or the value specified by option`GRID`

. For proper use of this threshold, other thresholds should be left at their default value of zero.Orbital truncation threshold.`ORBITAL`

Density truncation threshold.`DENSITY`

Fock matrix truncation threshold.`FOCK`

### Exact exchange computation (EXCHANGE)

`EXCHANGE`

,*factor*

For Kohn-Sham calculations, compute exchange energy according to Hartree-Fock formalism and add the contribution scaled by *factor* to the fock matrix and the energy functional. Otherwise, the default is *factor=0*, i.e., the exchange is assumed to be contained in the functional, and only the Coulomb interaction is calculated explicitly.

`DFTFACTOR`

,*fac1, fac2, …*

Provide a factor for each functional specified. The functionals will be combined accordingly. By default, all factors are one.

### Double-hybrid functionals (DH, DSDH)

`DH`

, *ax*, *ac*

initiates a double-hybrid calculation (Ref. [1]) where $a_x$ is the fraction of HF exchange and $a_c$ is the fraction of MP2 correlation. A self-consistent KS calculation is performed with functional $a_x E_x^{\text{HF}} + (1-a_x) E_x^{\text{DFT}}[\rho] + (1-a_c) E_c^{\text{DFT}}[\rho]$. One then needs to call the MP2 program to add the MP2 contribution $a_c E_c^{\text{MP2}}$. If $a_c$ is not given, $a_c=a_x^2$ is assumed, according to Ref. [2].

`DSDH`

, *ax*, *ac*

initiates a density-scaled double-hybrid calculation (Ref. [2]) where $a_x$ is the fraction of HF exchange and $a_c$ is the fraction of MP2 correlation. A self-consistent KS calculation is performed with functional $a_x E_x^{\text{HF}} + (1-a_x) E_x^{\text{DFT}}[\rho] + E_c^{\text{DFT}}[\rho] - a_c E_c^{\text{DFT}}[\rho_{1/\sqrt{a_c}}]$, where $\rho_{1/\sqrt{a_c}}$ is the scaled density. In the case of meta-GGA functionals, the kinetic energy density $\tau$ is also scaled, according to Ref. [3]. One then needs to call the MP2 program to add the MP2 contribution $a_c E_c^{\text{MP2}}$. If $a_c$ is not given, $a_c=a_x^2$ is assumed, according to Ref. [2].

Example of input for B2-PLYP calculation (Ref. [1]):

`{ks,b,lyp;dh,0.53,0.27;}`

`mp2;`

Example of input for 1DH-BLYP calculation (Ref. [2]):

`{ks,b,lyp;dh,0.65;}`

`mp2;`

Example of input for DS1DH-TPSS calculation (Ref. [3]):

`{ks,tpss;dsdh,0.725;}`

`mp2;`

Unrestricted calculations (`uks`

, `ump2`

) can also be done.

References:

$[1]$ S. Grimme, J. Chem. Phys. **124**, 034108 (2006).

$[2]$ K. Sharkas, J. Toulouse, A. Savin, J. Chem. Phys. **134**, 064113 (2011).

$[3]$ S. M. O. Souvi, K. Sharkas, J. Toulouse, J. Chem. Phys. **140**, 084107 (2014).

### Rangehybrid methods (RANGEHYBRID)

For coupling of short-range (sr-)DFT with long-range (lr-)ab-initio methods, one first has to specify the coupling parameter $\mu$ in the sr interelectronic interaction $\sum_{i<j} {\rm erf} (\mu r_{ij})/r_{ij}$; this can be done by setting a variable (e.g. `mu=0.5`

). As a next step, long-range ERIs have to be calculated by calling the integral program (e.g. `int;erf,mu;`

).

Then sr-DFT/lr-HF calculations can be performed by calling the RKS program with the additional subcommand `rangehybrid`

. Available short-range functionals are `exerf`

and `ecerf`

for sr-LDA, and `exerfpbe`

and `ecerfpbe`

for sr-PBE; as usual, the functionals have to be specified after the `rks`

command (e.g. `rks,exerf,ecerf;`

). The underlying short-range LDA correlation functional is that of S. Paziani, S. Moroni, P. Gori-Giorgi, G.B. Bachelet, Phys. Rev. B **73**, 155111 (2006).

Finally, sr-DFT/lr-post-HF calculations can be done by adding, within a call of the chosen post-HF program, two subcommands: `srxcdft`

followed by the desired short-range functionals (e.g. `srxcdft,exerf,ecerf;`

), and `dftden`

followed by the record number from which the density for the sr functionals is taken. Implementations are available for `ci`

, `mp2`

, `ccsd`

, `ccsd(t)`

, `rpatddft`

, and the corresponding local MP2 and CC methods w/wo density-fitting.

General two-parameter range-separated double hybrids (RSDH) as described in [1] are also available.

Example of a RSDH calculation with the approximation 3 of Ref. [1] for the short-range density functional:

mu=0.46 lambda=0.58 {int;erflerfc,mu=$mu,srfac=$lambda} {rks,exsrlpbe,ecsqrtlpbe;rangehybrid;orbital,7000.2} {mp2;srxcdft,exsrlpbe,ecsqrtlpbe;dftden,7000.2}

Example of a RSDH calculation with the approximation 4 of Ref. [1] for the short-range density functional:

mu=0.62 lambda=0.60 {int;erflerfc,mu=$mu,srfac=$lambda} {ks,exerfpbe,ecerfpbe;rangehybrid;dsrsdh;orbital,7000.2} mp2

Note that the `erflerfc`

option to the integral program can be used to calculate arbitrary mixtures of short-range and long-range Coulomb repulsion integrals via, e.g.,

mu=0.5 !range-separation parameter sr=0.5 !scaling factor for erfc interaction lr=0.5 !scaling factor for erf interaction {int;erflerfc,mu=$mu,srfac=$sr,lrfac=$lr}

which corresponds to the interaction operator $\mbox{lr}\cdot\frac{\mbox{erf}(\mu r_{12})}{r_{12}} + \mbox{sr}\cdot\frac{\mbox{erfc}(\mu r_{12})}{r_{12}}$. Setting `lr=1`

and `sr=1`

would
therefore mean that Molpro will calculate the full-range Coulomb integrals.

The MP2 calculation can also be replaced by a RPA calculation, like in the RS2H+RPAxSO2 method of Ref. [2].

Example of a RS2H+RPAxSO2 calculation:

mu=0.48 lambda=0.34 {int;erflerfc,mu,lambda} {rks,exsrlpbe,ecsqrtlpbe;rangehybrid;orbital,7000.2} {rpatddft;ecorr,SO2-RCCD;orb,7000.2}

Reference

$[1]$ C. Kalai and J. Toulouse, J. Chem. Phys. **148**, 164105 (2018).

$[2]$ C. Kalai, B. Mussard, and J. Toulouse, J. Chem. Phys. **151**, 074102 (2019).

### Exchange-correlation potential (POTENTIAL)

`POTENTIAL`

,*rec.fil*

For stand-alone DFT calculations, compute exchange-correlation potential pseudo-matrix elements, defined formally as the differential of the sum of all specified functionals with respect to elements of the atomic orbital density matrix. The matrix is written to record *rec* on file *fil*.

### Grid blocking factor (DFTBLOCK)

`DFTBLOCK`

,*nblock*

Respecify the number of spatial integration points treated together as a block in the DFT integration routines (default 128). Increasing *nblock* may enhance efficiency on, e.g., vector architectures, but leads to increased memory usage.

### Dump integrand values(DFTDUMP)

`DFTDUMP`

,*file,status*

Write out values of the integrand at grid points to the file *file*. The first line of *file* contains the number of functional components; there then follows a line for each functional giving the input key of the functional. Subsequent lines give the functional number, cartesian coordinates, integrand value and integration weight with Fortran format `(I2,3F15.10,F23.15)`

.

### Asymptotic correction for xc-potentials(ASYMP)

`ASYMP`

,*shift,$\alpha$,$\beta$,b*

Activates the gradient-regulated asymptotic correction (GRAC) approach for exchange-correlation potentials of Grüning *et al.* (J. Chem. Phys. **114**, 652 (2001)). The user has to supply a *shift* parameter ($\Delta_{\text{xc}}$) for the bulk potential which should approximate the difference between the HOMO energy ($\varepsilon_{\text{HOMO}}$) obtained from the respective standard Kohn-Sham calculation and the (negative) ionisation potential of the monomer ($\text{IP}$): $$\Delta_{\text{xc}}=\varepsilon_{\text{HOMO}}-(-\text{IP})$$ This method accounts for the derivative discontinuity of the exact xc-potential and that is missing in approximate ones. The parameters *$\alpha$* and *$\beta$* determine the interpolation function (see Eq. (2.3) in the above reference) and are set to $\alpha=0.5$ and $\beta=40$ by default, respectively. The parameter *b* is the parameter of the asymptotic xc-potential from van Leeuwen and Baerends (Phys. Rev. A **49**, 2421 (1994), Eqns. (54,55)) and is set to *b*=0.05 by default.

In case of gradient or laplacian functionals the modified GRAC scheme of Bast *et al.* (Chem. Phys. Chem. **9**, 445 (2008)) is used.

If *shift* is set to zero in the input the program will estimate the ionisation energy from the HOMO energy during the SCF (as soon as the HOMO energy is converged to a given threshold) and then sets the bulk shift automatically. This is done by using a linear fit of DFT HOMO energies to ionisation energies calculated with the $\Delta$SCF method for a range of molecules (see also S. Hirata *et al.*, J. Phys. Chem. A, **107**, 10154 (2003)).

## Numerical integration grid control (GRID)

Density functionals are evaluated through numerical quadrature on a grid in three-dimensional space. Although the sensible defaults will usually suffice, the parameters that define the grid can be specified by using the `GRID`

command, which should be presented *before* the `DFT`

or `KS`

commands that will use the grid, but after the geometry definition. Typically, the input file then could look like

{grid,<grid options>; <grid directives>} {ks,...}

All grid parameters that are modified through the options and directives of the `GRID`

program are changed globally, i.e., any subsequent program using a quadrature grid in Molpro will use the grid determined by GRID unless they use any local redefinitions of grid parameters. Alternatively, any available
options for `GRID`

can be used within the `KS`

program:

{ks,...,<grid options>}

Note that then, however, the modifications to the grid parameters are only applied during the KS calculation but they are reset to their default state afterwards. Also note that the directives which are available for `GRID`

can not be used with `KS`

, so these can be exclusively only used for `GRID`

.
The grid parameters are stored internally on a record whose location can be changed by

`GRID,[record=`

*orb.file*`],[status=`

*status*`]`

By default Molpro uses the record *orb.file* = 1800.2 to store the grid parameters that are independent of the geometry. The information on disk consists of two parts: the parameters necessary to define the grid (record on file 2), and a cache of the evaluated grid points and weights (corresponding record on file 1, usually 1800.1). The latter is flagged as ‘dirty’ whenever any parameters are changed, and whenever the geometry changes; if the cache is dirty, then when an attempt is made to use the grid, it will be recalculated, otherwise the cached values are used.

If *status* is `OLD`

, an attempt to restore the grid from a previous calculation is performed; effectively, the old grid provides a template of parameters which can be adjusted using the parameter commands described below. If *status* is `NEW`

, the grid is always created with default parameters. If *status* is `UNKNOWN`

(the default), a new grid is created if record *orb.file* does not exist; otherwise the old grid is used.

The `GRID`

command may be followed by a number of parameter-modifying directives. The currently implemented default parameters are equivalent to the following input commands.

RADIAL,LOG,3,1.0,20,25,25,30 ANGULAR,LEBEDEV,0.0,0.0 LMIN,0,0,0,0 LMAX,53,53,53,53 VORONOI,10 GRIDSAVE GRIDSYM

Options that have a single values directives can be given as an option to `GRID`

:

{GRID,GRIDTHR=...,FACNEIGH=...} or {KS,...; GRID,GRIDTHR=...,FACNEIGH=...}

where `GRIDTHR`

just takes a scalar value for the global grid threshold value, so using the corresponding directive will be still necesary if
different thresholds for the radial and angular grid shall be used, see below.

Molpro can use both an internal implementation to compute the quadrature grid (termed *dftgrid* here) and an external library called *libxcgrid*. For adaptive (pruned) grids, that are computed by default, one can switch between the two cases using

{grid,name=PRUNED}

for using *xcgrid* and

{grid,name=OLD}

to reset to *dftgrid* again. *xcgrid* in addition also contains implementations for some fixed quadrature grids, see below. Most of the options and directives described in the following can be used in conjunction with both programs. However, some can only be exlusively used with either of the two codes. To highlight this in the following, the *dftgrid*-exclusive options will be marked by whereas the *libxcgrid*-exclusive ones by .

### Summary of options and directives for GRID

Almost all grid parameters can be controlled through respective options given on the GRID card:

{grid,option1=value1,option2=value2,...}

A limited number of subcommands (directives) are also available that can be used in the input like:

{grid; directive1,val1a,val1b,...; directive2,val2a,val2b,}

where *val1a*,*val1b*,… are lists of comma-separated values for the first directive, etc.

The following table summarises the options that can be used for GRID. The last two columns displays if the respective option is available with the and programs. A detailed description of the different options will be given in the upcoming sections.

Options for GRID | |||||
---|---|---|---|---|---|

group | name | default value | description | ||

grid type | `NAME` | old | select pruned or fixed standard grids | ✔ | ✔ |

accuracy | `GRIDTHR|THRESHOLD_OVERALL|GRID|THRGRID` | acc=1d-6 | target grid threshold | ✔ | ✔ |

`GRIDTHR_RAD|THRESHOLD_RADIAL|ACCR` | acc/2 | target threshold for radial grid | ✔ | ✔ | |

`GRIDTHR_ANG|THRESHOLD_ANGULAR|ACCA` | acc/10 | target threshold for angular grid | ✔ | ✔ | |

`GRIDACC|MODE` | default | global mode for grid accuracy | ✔ | ✔ | |

`ANGLEVEL` | ✘ | use fixed angular grids in three-zone decomposition of radial grid | ✔ | ||

radial integration grid | `RADIAL_SCHEME` | log | radial integration scheme | ✔ | ✔ |

`GRIDTHR_RAD|THRESHOLD_RADIAL|ACCR` | acc/2 | target threshold for radial grid | ✔ | ✔ | |

`RADIAL_SCALING` | ✘ | scaling factors for atomic radial grids | ✔ | ||

`MIN_NR` | [20,25,25,30] | lower bound for degrees of quadrature for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

`MAX_NR` | [500,500,500,500] | upper bound for degrees of quadrature for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

angular integration grid | `ANGULAR_SCHEME` | lebedev | angular integration scheme | ✔ | ✔ |

`THRESHOLD_ANGULAR|GRIDTHR_ANG|ACCA` | acc/10 | target threshold for angular grid | ✔ | ✔ | |

`ANGLEVEL` | ✘ | use fixed angular grids in three-zone decomposition of radial grid | ✔ | ||

`MIN_L` | [0,0,0,0] | minimum angular order for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

`MAX_L` | [53,53,53,53] | maximum angular order for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

atom partitioning | `VORONOI_SCHEME` | murray | type of step function | ✔ | ✔ |

`SIZEADJ` | treutler | size adjustment method | ✔ | ✔ | |

`BECKE_MMU` | 3 | $m_\mu$ recursion value in Becke scheme | ✔ | ✔ | |

`MURRAY_MMU` | 10 | $m_\mu$ recursion value in Murray scheme | ✔ | ✔ | |

`PRUNING_VORONOI_PAIRS` | all | limit Voronoi pairs taken into account for pruning | ✔ | ||

pruning | `PRUNING_INTEGRAND` | PW91MOL | model functional used for angular pruning | ✔ | |

`FACNEIGH|FAC_NEIGHBOUR|FAC_NEIGHBOR` | 4.0 | scaling factor for the selection of neighbour atoms for the calculation of the test density | ✔ | ✔ | |

`FROM_MIN_ORDER` | true | order of pruning (starting from lowest or highest value of $L$) | ✔ | ||

fixed standard grids | `KK_INDEX` | 4 | index for Krack-Köster type grid | ✔ | |

`KK_MAX_ORDER` | 29 | maximum order for KK type grid when `KK_INDEX` < 3 | ✔ | ||

`KK_INT_ACC` | 4.67 | Krack-Koster parameter for number of radial points for Krack-Köster type grid when `KK_INDEX` < 3 | ✔ | ||

other | `WEIGHT_CUT|WCUT` | 1d-20 | threshold value to discard grid points with small integration weights | ✔ | ✔ |

Directives for GRID | |||||
---|---|---|---|---|---|

group | name | default values | description | ||

radial integration grid | `RADIAL` | log,3,1.0,20,25,25,30 | method and parameters for radial grid | ✔ | ✔ |

angular integration grid | `ANGULAR` | lebedev,0.0,0.0 | method and parameters for angular grid | ✔ | ✔ |

`LMIN` | 0,0,0,0 | minimum angular order for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

`LMAX` | 53,53,53,53 | maximum angular order for H,He, 2nd row, 3rd row, other elements | ✔ | ✔ | |

atom partitioning | `VORONOI` | murray,10,2 | step function method, $m_\mu$ recursion value and type of weights | ✔ | ✔ |

other | `GRIDSAVE` | ✘ | enable grid caching | ✔ | ✔ |

`NOGRIDSAVE` | ✘ | disable grid caching | ✔ | ✔ | |

`GRIDSYM` | ✘ | force use of any point-group symmetry | ✔ | ||

`NOGRIDSYM` | ✘ | disable use of any point-group symmetry | ✔ | ||

`GRIDPRINT` | grid=0 | print grid points at different printing levels | ✔ | ✔ |

### Target quadrature accuracy (GRIDTHR)

`GRIDTHR|THRESHOLD_OVERALL|GRID|THRGRID`

=*acc*

`GRIDTHR_RAD|THRESHOLD_RADIAL|ACCR`

=*accr*

`GRIDTHR_ANG|THRESHOLD_ANGULAR|ACCA`

=*acca*

Specify the target accuracy of integration. If options `GRIDTHR_RAD`

and `GRIDTHR_ANG`

are missing, radial and angular grids are generated adaptively, with the aim of integrating a model functional to the specified accuracy. *acc* is an overall target accuracy, and is the one that should normally be used; radial and angular grid target accuracies are generated algorithmically from it, namely, the scaled values *accr*=*acc*/2 and *acca*=*acc*/10 will then be used. However, they can be adjusted individually by using the options `GRIDTHR_RAD`

and `GRIDTHR_ANG`

respectively.

### Setting a global level for the accuracy of the grid (GRIDACC)

`GRIDACC|MODE`

=[*old*,*default*,*low*,*high*,*veryhigh*]

The accuracy of the numerical quadrature depends on a number of different parameters whose effect may not be very transparent for the user. However, a fine-tuning of the grid parameters in order to control the accuracy of the grid is not necessary when using the global option `GRIDACC`

that can have the four values *default*,*low*,*high* and *veryhigh*. `GRIDACC`

=*default* uses sensible settings for the grid parameters that should suffice for most kinds of calculations and is set by default. In contrast, `GRIDACC`

=*low* may be used when the accuracy of the DFT integration is not so relevant (e.g., for computing starting guesses) and `GRIDACC`

=*high* or `GRIDACC`

=*veryhigh* are recommended for benchmark calculations. The setting of `GRIDACC`

=*old* is for keeping backwards compatibility with older versions of Molpro (< 2021.x) where the default values for the grid parameters differed from the current ones.

### Radial integration grid (RADIAL)

`RADIAL`

,*method,accr,$m_r$,scale,$n_0,n_1,n_2,n_3$*

`RADIAL_SCHEME`

=*method*

`THRESHOLD_RADIAL|GRIDTHR_RAD|ACCR`

=*accr*

`RADIAL_SCALING`

=[*atom1*,*atom2*,…]

`MIN_NR`

,$n_0,n_1,n_2,n_3$

`MAX_NR`

,$n_0,n_1,n_2,n_3$

Specify the details of the radial quadrature scheme. A number of different radial schemes are available, specified by *method* = `EM`

, `BECKE`

, `AHLRICHS`

, `MULTIEXP`

, `DE2`

, `MK2CG`

or `LOG`

, with the latter being the default.

`EM|EULER-MACLAURIN`

is the Euler-Maclaurin scheme defined by C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys. **78**, 997 (1993). $m_r$, for which the default value is 2, is defined in equation (6) of the above as $$r = \alpha {x^{m_r}\over (1-x)^{m_r}}$$ whilst *scale* (default value 1) multiplied by the Bragg-Slater radius of the atom gives the scaling parameter $\alpha$.

`LOG|MURA-KNOWLES`

is the scheme described by M. E. Mura and P. J. Knowles, J. Chem. Phys. **104**, 9848 (1996). It is based on the transformation $$r = - \alpha \log_e (1-x^{m_r})\; ,$$ with $0\le x \le 1$ and simple Gauss quadrature in $x$-space. The recommended value of $m_r$ is 3 for molecular systems, giving rise to the Log3 grid; $m_r$=4 is more efficient for atoms. $\alpha$ is taken to be *scale* times the recommended value for $\alpha$ given by Mura and Knowles, and *scale* defaults to 1.

`BECKE`

is as defined by A. D. Becke, J. Chem. Phys. **88**, 2547 (1988). It is based on the transformation $$r = \alpha {(1+x)\over (1-x)} \; ,
\label{eqn:Becke}$$ using points in $-1\le x \le +1$ and standard Gauss-Chebyshev quadrature of the second kind for the $x$-space quadrature. Becke chose his scaling parameters to be half the Bragg-Slater radius except for hydrogen, for which the whole Bragg-Slater radius was used, and setting *scale* to a value other than 1 allows a different $\alpha$ to be used. $m_r$ is not necessary for this radial scheme.

`AHLRICHS`

is the radial scheme defined by O. Treutler and R. Ahlrichs, J. Chem. Phys. **102**, 346 (1995). It is based on the transformation (their M4 mapping) $$r= {\alpha \over \log_e 2} (1+x)^{0.6} \log_e\left( {2\over 1-x}\right)\; ,$$ with using standard Gauss-Chebyshev quadrature of the second kind for the $x$-space integration. $m_r$ is not necessary for this radial scheme.

`CHEBY`

2nd kind Chebyshev quadrature

`MULTIEXP`

J. Comput. Chem. **24**, 732 (2003)

`MK2GC`

Gauss-Chebyshev quadrature with abscissas $$x_i = \cos\left(i \pi \over n+1 \right)$$

`DE2|DOUBLEXP`

Double Exponential quadrature of Mitani *et al.* Theor. Chem. Acc. **130**, 645 (2011), Theor. Chem. Acc. **131**, 1169 (2012) based on the transformation $$r=\exp(\alpha x_i-e^{-x_i}) $$

$n_0$, $n_1$, $n_2$, $n_3$ are the degrees of quadrature $n_r$ (see equation (3) of Murray et al.), for hydrogen/helium, first row, second row, and other elements respectively.

`accr`

specifies a target accuracy; the number of radial points is chosen according to a model, instead of using an explicit $n_i$. The stricter of $n_i$, `accr`

is used, unless either is zero, in which case it is ignored. Molpro will then set `accr`

=`acc`

/2d0 where `acc`

is the global threshold for the grid accuracy that can be set with the option `GRIDTHR`

.

An atom dependent radial scaling factor can be set with the option `RADIAL_SCALING`

, e.g. `RADIAL_SCALING`

=[1.0,1.1,1.1] could be a possible modification of the radial grids for H$_2$O . Note that factors for all atoms must be given, otherwise the option will be ignored. The option can also be used to modify the scalar *scale* factor for in which case only a single value needs to be given that will then be used to scale the radial grids for all atoms.

The size of the radial grid is evaluated from the `GRIDTHR_RAD`

option:
$$n_{\text{rad}}(A) = -\log_{10} accr - 30 + (n_A - 1) * 15$$
where $accr$ is the threshold and $n_A$ is the period of atom.

The minimum and maximum values of the number of radial points is also restricted by :

`GRID,name=PRUNED,min_nr=[Row1,Row2,Row3,Row4+]`

and

`GRID,name=PRUNED,max_nr=[Row1,Row2,Row3,Row4+]`

for each atom based on its periodic row.

### Angular integration grid (ANGULAR)

`ANGULAR`

,*method,acca,crowd*

`LMIN`

,$l^{\text{min}}_0,l^{\text{min}}_1,l^{\text{min}}_2,l^{\text{min}}_3$

`LMAX`

,$l^{\text{max}}_0,l^{\text{max}}_1,l^{\text{max}}_2,l^{\text{max}}_3$

`MIN_L`

=[Row1,Row2,Row3,Row4+]

`MAX_L`

=[Row1,Row2,Row3,Row4+]

`ANGULAR_SCHEME`

=*lebedev*|*legendre*|*lobatto*

`ANGLEVEL`

=*1*-*8*

`THRESHOLD_ANGULAR|GRIDTHR_ANG|ACCA`

=*acca*

Specify the details of the angular quadrature scheme. The default choice for *method* is `LEBEDEV`

(ie. as in A. D. Becke, J. Chem. Phys. **88**, 2547 (1988)) which provides angular grids of octahedral symmetry. An alternative choice for *method* is `LEGENDRE`

which gives Gauss-Legendre quadrature in $\theta$ and simple quadrature in $\phi$, as defined by C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys. **78**, 997 (1993). Similarly, *method*=`LOBATTO`

defines a spherical quadrature over $\theta$ and $\phi$ space using Lobatto's formula and using the size reduction described by Treutler and Ahlrichs J. Chem. Phys. **102**, 346 (1995). Both, *method*=`LEGENDRE`

and `LOBATTO`

are suitable options if highly accurate angular integrations need to be performed.

Each type of grid specifies a family of which the various members are characterized by a single quantum number $l$; spherical harmonics up to degree $l$ are integrated exactly. $l^{\text{min}}_i$ and $l^{\text{max}}_i, i=0,1,2,3$ specify allowed ranges of $l$ for H–Be, B–Ca, Sc–Ba, and La– respectively. The $l^{\text{min}}_i$ are further moderated at run time so that for any given atom they are not less than $2i+4$ or twice the maximum angular momentum of the basis set on the atom; this constraint can be overridden by giving a negative value in `LMIN`

, and in this case just its absolute value will be used as the lower bound. For the Lebedev grids, if the value of $l$ is not one of the set implemented in Molpro (3, 5, 7, 9, 11, 13, 15, 17, 19, 23, 29, 41, 47, 53), then $l$ is increased to give the next largest angular grid available. In general, different radial points will have different $l$, and in the absence of any moderation described below, will be taken from $l^{\text{max}}_i$.

The maximum allowed angular order is set by the row in the periodic table: `GRID,name=PRUNED,max_l=[Row1,Row2,Row3,Row4+]`

. The default value is [53,53,53,53]. Additionally, the minimum angular order can be set: `GRID,name=PRUNED,min_l=[Row1,Row2,Row3,Row4+]`

. The default value is [3,3,3,3]. This minimum is also modified to be at least $2*l+1$ where $l$ is the maximum angular momentum of the basis set.

*crowd* is a parameter to control the reduction of the degree of quadrature close to the nucleus, where points would otherwise be unnecessarily close together; larger values of crowd mean less reduction thus larger grids. A very large value of this parameter, or, conventionally, setting it to zero, will switch off this feature.

*acca* is a target energy accuracy. It is used to reduce $l$ for a given radial point as far as possible below $l^{\text{max}}_i$ but not lower than $l^{\text{max}}_i$. The implementation uses the error in the angular integral of the kernel of a model functional using a sum of approximate atomic densities. If *acca* is zero, the global threshold is used instead, or else it is ignored.

`ANGLEVEL`

is an option that defines a fixed Lebedev order for three different radial zones defined by the ranges [1…*N*$_{rad}$/3], [*N*$_{rad}$/3+1…*N*$_{rad}$/2] and [*N*$_{rad}$/2+1…*N*$_{rad}$] with *N*$_{rad}$ being the number of radial quadrature points and the order is from close to nucleus to far distant radial points. The idea of this definition is that only low Lebedev orders are needed in the first and second zone while higher orders are required in the outer third zone in order to achieve reasonable accuracies for the integration. It is possible to choose between seven different schemes `ANGLEVEL`

=1 to `ANGLEVEL`

=7 where the grid sizes increases in that order. Fairly accurate results should be accessible with intermediate levels `ANGLEVEL`

=3,4. It is also possible to select `ANGLEVEL`

=8 in which case the $L=59$ (1202-point) Lebedev grid is used at each radial point for highly accurate integrations.

### Pruning integrand and neighbour atoms in grid generation (FACNEIGH)

`PRUNING_INTEGRAND`

=*SAD*,*PW91*,*PW91MOL*

`FROM_MIN_ORDER`

=0,1

`FACNEIGH|FAC_NEIGHBOUR|FAC_NEIGHBOR`

=*facn*

The grid generation for the angular grid depends on the integration of a test density functional around a sphere around each atom at a given radial point.
The procedure is that first a target value for the highest possible Lebedev order $L$ is computed and then that starting from the lowest value of $L$ (up-pruning) the integral is computed again until the deviation to the target value gets smaller than the threshold specified by the value for `GRIDTHR_ANG`

=acca, see also section Angular integration grid. It is also possible to start from the largest possible value of $L$ by setting `FROM_MIN_ORDER`

=0 (down-pruning) and so the optimal value of $L$ will then be found by seeking for the first one where the angular threshold is no longer obeyed. This option is expected to be more accurate but can also lead to larger grid sizes than the default option of `FROM_MIN_ORDER`

=1.

One can choose between the three different integrands `PRUNING_INTEGRAND`

=*SAD*,*PW91*,*PW91MOL* . *SAD* just defines a superposition of Slater atomic densities, *PW91* is a pseudo-density functional based on the PW91 xc functional and *PW91MOL* is a refinement of *PW91* where the procedure to find the optimal value of $L$ is slightly different (e.g., a linear regression is performed in order to determine $L$ with this integrand). Note that *PW91MOL* is set by default and is available both with and . In contrast, integrands *SAD* and *PW91* can only be used with .

The test functional is computed from the superposition of (Slater) atomic densities including the atom itself, plus a set of neighbour atoms. The latter includes all atoms that are within a range of `DIST`

*`FACNEIGH`

apart from the atom with `DIST`

being the (squared!) distance of the closest atom. `FACNEIGH`

is currently set to 4.0 by default and governs the shape of the test density on the sphere. Increasing this value will implicitly lead to more angular points at (particularly) large radial grid points but also increases the overall accuracy of the quadrature. For benchmark calculations it can be therefore recommended to set `FACNEIGH`

larger than the current threshold. Note that the same can also be achieved by reducing the value of the overall `GRIDTHR`

threshold value, yet, it has been found that resetting `FACNEIGH`

can be more effective for intermediate sized grids.

### Atom partitioning of integration grid (VORONOI)

`VORONOI`

,$m_\mu$,`vortyp`

`VORONOI_SCHEME`

=*murray*,*becke*,*stratmann*

`SIZEADJ`

=*becke*,*treutler*,*none*

`BECKE_MMU`

=*mmu*

`MURRAY_MMU`

=*mmu*

`PRUNING_VORONOI_PAIRS`

=*all*,*atom*

Controls Becke-Voronoi partitioning of space. By default, the algorithm of C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys. **78**, 997 (1993) is used, with $m_\mu$ defined by equation (24). The default value is 10.
`vortyp`

defines the type of grid weights to be computed: 1: raw weights, 2: full voronoi, 4: pair voronoi. The default is
`vortyp`

=2.

It is possible to choose between the three different step functions `VORONOI_SCHEME`

=*murray*,*becke* and *stratmann*.
The Murray and Becke scheme employ a recursive formula that increases the steepness of the function. The recursion limit is
given by its *mmu* value. In case of *becke* the default value of $m_\mu$=3 is being used. This can be modified from the default for each case using the options `BECKE_MMU`

or `MURRAY_MMU`

.

For heteronuclear molecules Molpro uses Treutler and Ahlrichs atomic size adjustment for the smoothing function
J. Chem. Phys. **102**, 346 (1995) (Eq.(13)). Alternatively, the scheme by Becke
J. Chem. Phys. **88**, 2547 (1988).
can be used by setting `SIZEADJ`

=*becke*. Possible choices for `SIZEADJ`

are *becke*, *treutler* and *none*, the latter meaning that no size adjustment is being used.

On-the-fly grid pruning can be expensive for large molecules because of the cubic scaling of Voronoi partitioning. Pruning is performed in under mpp parallelisation, with each process pruning each atom. If it is still too expensive, Voronoi partitioning can be limited for the pruning algorithm to quadratic scaling: `GRID,name=PRUNED,pruning_voronoi_pairs=ATOM`

but can result in a less accurate grid.

### Discarding grid points with tiny weights (WEIGHT_CUT)

`WEIGHT_CUT|WCUT`

=*thr*

Grid points very close to the nucleus can have very small grid weights which can drop below machine precsion (~10d-16) when the radial grid is very large. These can be discarded with the option `WEIGHT_CUT`

=*thr*, i.e., grid points with weights smaller than *thr* will then not be used in the numerical integration anymore. While this can make the integration a little bit more efficient, a truncation of the grid with this option can prevent to achieve high accuracies in the numerical integration below, e.g., micro-hartree accuracies for energies. Therefore by default `WEIGHT_CUT`

is set to a very small number that is well below machine precision (`WEIGHT_CUT`

=1d-20) so that effectively all grid points close to the nucleus are taken into account for the integration.

### Grid caching (GRIDSAVE, NOGRIDSAVE)

`NOGRIDSAVE`

disables the disk caching of the grid, i.e, forces the recalculation of the grid each time it is needed.

`GRIDSAVE`

forces the use of a grid cache where possible.

### Grid symmetry (GRIDSYM,NOGRIDSYM)

`NOGRIDSYM`

switches off the use of symmetry in generating the integration grid, whereas

`GRIDSYM`

forces the use of any point-group symmetry.

### Grid printing (GRIDPRINT)

`GRIDPRINT`

,*key=value*,$\dots$

controls printing of the grid, which by default is not done. At present, the only possible value for *key* is `GRID`

, and *value* should be specified as an integer. `GRID=0`

causes the total number of integration points to be evaluated and reported; `GRID=1`

additionally shows the number of points on each atom; `GRID=2`

causes the complete set of grid points and weights to be printed.

## Fixed standard grids

There are three alternative types of pruned grids that can be used; `PRUNED`

, `KK`

and `SG`

. They can be selected by setting the `name`

option:

`GRID,name=PRUNED`

or

`GRID,name=KK`

or

`GRID,name=SGx`

with `SGx`

=`SG0`

,`SG1`

,`SG2`

or `SG3`

.

`KK`

and `SG`

are fixed-pruned grids and are defined in sections KK Grids (GRID,NAME=KK) and SG Grids (GRID,name=SG) respectively. `PRUNED`

grids compute angular order pruning on-the-fly based on the molecular geometry. The `PRUNED`

grid also features the most diverse set of options that can be tuned to balance accuracy and grid size, see sections Radial integration grid, Angular integration grid and Atom partitioning of integration grid.

### Krack-Köster Grids (GRID,name=KK)

The KK grid is a fixed-pruned grid.

The size of the grid should set using the `kk_index`

:

`grid,name=KK,kk_index=value`

which links the angular size using `kk_max_order`

and radial size using `kk_int_acc`

, which are defined below. The `value`

can be set from 3 to 12. Using `kk_index=3`

sets `kk_int_acc=4.34`

and `kk_max_order=23`

. Increasing the index by 1 increases `kk_int_acc`

by 0.33 (effectively 5 radial points) and `kk_max_order`

by 6.

The grid is defined with 5 shells with boundaries at fixed radii. The angular order of each shell is defined relative to a maximum order, $m$, as $m-18, m-12, m-6, m, m-6$. It can be set using:

`grid,name=KK,kk_max_order=value`

Ahlrichs quadrature is used for the radial grid. The number of points used for atom $A$ is defined by the Krack-Köster formula [1]: $$n_{\text{rad}}(A) = -5 (3 \log 10^{-\tt{int\_acc}} - n_A + 8 )$$ where $n_A$ is the period of atom $A$. `int_acc`

is set using:

`grid,name=KK,kk_int_acc=value`

where `value`

would take values from 4.34 upwards.

References:

$[1]$ Krack, M., Köster, A.M., 1998. An adaptive numerical integrator for molecular integrals. J. Chem. Phys. **108**, 3226-3234. https://doi.org/10.1063/1.475719

### SG Grids (GRID,name=SG)

The fixed-pruned Standard Grid (SG) can be used:

`GRID,name=SG0`

`GRID,name=SG1`

`GRID,name=SG2`

or

`GRID,name=SG3`

`SG0`

is a small fixed-pruned grid by Chien and Gill[1]. It is defined for atoms in the first three rows except He, Ne and Ar. Less than 1500 points is produced per atom. The MultiExp radial grid and Becke’s Voronoi partitioning step function are used. The 18 point Lebedev angular grids are replaced with 14 point ones.

`SG2`

and `SG3`

are larger fixed grids by Dasgupta and Herbert [3] that were developed to achieve high accuracies for chemical properties. They are based on the `DE2`

(DoubleExp) radial quadrature scheme, see section Radial integration grid.

References:

$[1]$ Chien, S.-H., Gill, P.M.W., 2006. SG-0: A small standard grid for DFT quadrature on large systems. J. Comput. Chem. **27**, 730-739. https://doi.org/10.1002/jcc.20383

$[2]$ Gill, P.M.W., Johnson, B.G., Pople, J.A., 1993. A standard grid for density functional calculations. Chem. Phys. Let. 209, 506-512. https://doi.org/10.1016/0009-2614(93)80125-9

$[3]$ Dasgupta, S. Herbert, J. M., 2017. Standard Grids for High-Precision Integration of Modern Density Functionals: SG-2 and SG-3. J. Comput. Chem. **38**, 869–882. https://doi.org/10.1002/jcc.24761

## Density Functionals

Widely used functional combinations can be specified via convenient keywords. These can be entered, for example via `{df-rks,pbe}`

(do a single KS calculation with PBE) or `{dfunc,pbe}`

(which will instruct all following KS calculations to use PBE by default). The following table lists the defined alias density functional combinations:

alias | Ref | ||
---|---|---|---|

`B ` | `B88` dftfun:B88 | 1 | 10.1103/PhysRevA.38.3098 |

`B-LYP ` | `B88` dftfun:B88:`LYP` dftfun:LYP | 1:1 | |

`BLYP ` | `B88` dftfun:B88:`LYP` dftfun:LYP | 1:1 | |

`B-P ` | `B88` dftfun:B88:`P86` dftfun:P86 | 1:1 | |

`BP86 ` | `B88` dftfun:B88:`P86` dftfun:P86 | 1:1 | |

`B-VWN ` | `B88` dftfun:B88:`VWN5` dftfun:VWN5 | 1:1 | |

`B3LYP ` | `EXACT` dftfun:EXACT:`B88` dftfun:B88:`DIRAC` dftfun:DIRAC:`LYP` dftfun:LYP:`VWN5` dftfun:VWN5 | 0.2:0.72:0.08:0.81:0.19 | |

`B3LYP3 ` | `EXACT` dftfun:EXACT:`B88` dftfun:B88:`DIRAC` dftfun:DIRAC:`LYP` dftfun:LYP:`VWN3` dftfun:VWN3 | 0.2:0.72:0.08:0.81:0.19 | |

`B3LYP5 ` | `EXACT` dftfun:EXACT:`B88` dftfun:B88:`DIRAC` dftfun:DIRAC:`LYP` dftfun:LYP:`VWN5` dftfun:VWN5 | 0.2:0.72:0.08:0.81:0.19 | |

`B88X ` | `B88` dftfun:B88 | 1 | 10.1103/PhysRevA.38.3098 |

`B97 ` | `EXACT` dftfun:EXACT:`B97DF` dftfun:B97DF | 0.1943:1 | |

`B97R ` | `EXACT` dftfun:EXACT:`B97RDF` dftfun:B97RDF | 0.21:1 | |

`BECKE ` | `B88` dftfun:B88 | 1 | 10.1103/PhysRevA.38.3098 |

`BH-LYP ` | `EXACT` dftfun:EXACT:`B88` dftfun:B88:`LYP` dftfun:LYP | 0.5:0.5:1 | |

`CS ` | `CS1` dftfun:CS1 | 1 | |

`D ` | `DIRAC` dftfun:DIRAC | 1 | |

`HFB ` | `B88` dftfun:B88 | 1 | 10.1103/PhysRevA.38.3098 |

`HFS ` | `DIRAC` dftfun:DIRAC | 1 | |

`LDA ` | `DIRAC` dftfun:DIRAC:`VWN5` dftfun:VWN5 | 1:1 | |

`LSDAC ` | `PW92C` dftfun:PW92C | 1 | 10.1103/PhysRevB.45.13244 |

`LSDC ` | `PW92C` dftfun:PW92C | 1 | 10.1103/PhysRevB.45.13244 |

`LYP88 ` | `LYP` dftfun:LYP | 1 | |

`MM05 ` | `EXACT` dftfun:EXACT:`M05X` dftfun:M05X:`M05C` dftfun:M05C | 0.28:0.72:1 | 10.1063/1.2126975 |

`MM05-2X ` | `EXACT` dftfun:EXACT:`M052XX` dftfun:M052XX:`M052XC` dftfun:M052XC | 0.56:0.44:1 | 10.1021/ct0502763 |

`MM06 ` | `EXACT` dftfun:EXACT:`M06X` dftfun:M06X:`M06C` dftfun:M06C | 0.27:0.73:1 | 10.1007/s00214-007-0310-x |

`MM06-2X ` | `EXACT` dftfun:EXACT:`M062XX` dftfun:M062XX:`M062XC` dftfun:M062XC | 0.54:0.46:1 | 10.1007/s00214-007-0310-x |

`MM06-L ` | `M06LX` dftfun:M06LX:`M06LC` dftfun:M06LC | 1:1 | 10.1063/1.2370993 |

`MM06-HF ` | `EXACT` dftfun:EXACT:`M06HFX` dftfun:M06HFX:`M06HFC` dftfun:M06HFC | 1:1:1 | 10.1021/jp066479k |

`PBE ` | `PBEX` dftfun:PBEX:`PBEC` dftfun:PBEC | 1:1 | 10.1103/PhysRevLett.77.3865 |

`PBE0 ` | `EXACT` dftfun:EXACT:`PBEX` dftfun:PBEX:`PBEC` dftfun:PBEC | 0.25:0.75:1 | 10.1063/1.478522 |

`PBE0MOL ` | `EXACT` dftfun:EXACT:`PBEX` dftfun:PBEX:`PW91C` dftfun:PW91C | 0.25:0.75:1 | |

`PBEREV ` | `PBEXREV` dftfun:PBEXREV:`PBEC` dftfun:PBEC | 1:1 | |

`PW91 ` | `PW91X` dftfun:PW91X:`PW91C` dftfun:PW91C | 1:1 | 10.1103/PhysRevB.46.6671 |

`S ` | `DIRAC` dftfun:DIRAC | 1 | |

`S-VWN ` | `DIRAC` dftfun:DIRAC:`VWN5` dftfun:VWN5 | 1:1 | |

`SLATER ` | `DIRAC` dftfun:DIRAC | 1 | |

`VS99 ` | `VSXC` dftfun:VSXC | 1 | |

`VWN ` | `VWN5` dftfun:VWN5 | 1 | |

`VWN80 ` | `VWN5` dftfun:VWN5 | 1 | |

`M05 ` | `EXACT` dftfun:EXACT:`XC-M05` dftfun:XC-M05 | 0.28:1 | 10.1063/1.2126975 |

`M05-2X ` | `EXACT` dftfun:EXACT:`XC-M05-2X` dftfun:XC-M05-2X | 0.56:1 | 10.1021/ct0502763 |

`M06 ` | `EXACT` dftfun:EXACT:`XC-M06` dftfun:XC-M06 | 0.27:1 | 10.1007/s00214-007-0310-x |

`M06-2X ` | `EXACT` dftfun:EXACT:`XC-M06-2X` dftfun:XC-M06-2X | 0.54:1 | 10.1007/s00214-007-0310-x |

`M06-L ` | `XC-M06-L` dftfun:XC-M06-L | 1 | 10.1063/1.2370993 |

`M06-HF ` | `EXACT` dftfun:EXACT:`XC-M06-HF` dftfun:XC-M06-HF | 1:1 | 10.1021/jp066479k |

`M08-HX ` | `EXACT` dftfun:EXACT:`XC-M08-HX` dftfun:XC-M08-HX | 0.5223:1 | |

`M08-SO ` | `EXACT` dftfun:EXACT:`XC-M08-SO` dftfun:XC-M08-SO | 0.5679:1 | |

`M11-L ` | `XC-M11-L` dftfun:XC-M11-L | 1 | |

`TPSS ` | `TPSSX` dftfun:TPSSX:`TPSSC` dftfun:TPSSC | 1:1 | 10.1103/PhysRevLett.91.146401 |

`TPSSH ` | `EXACT` dftfun:EXACT:`TPSSX` dftfun:TPSSX:`TPSSC` dftfun:TPSSC | 0.1:0.9:1 | 10.1063/1.1626543 |

`M12HFC ` | `EXACT` dftfun:EXACT:`M12C` dftfun:M12C | 1:1 | 10.1063/1.4768228 |

`HJSWPBE ` | `HJSWPBEX` dftfun:HJSWPBEX:`PBEC` dftfun:PBEC | 1:1 | 10.1063/1.2921797 |

`HJSWPBEH ` | `EXACT` dftfun:EXACT:`HJSWPBEX` dftfun:HJSWPBEX:`PBEC` dftfun:PBEC | 0.2:0.8:1 | 10.1063/1.3073302 |

`TCSWPBE ` | `EXERFPBE` dftfun:EXERFPBE:`ECERFPBE` dftfun:ECERFPBE | 1:1 | 10.1063/1.1824896 |

`PBESOL ` | `PBESOLX` dftfun:PBESOLX:`PBESOLC` dftfun:PBESOLC | 1:1 | 10.1103/PhysRevLett.100.136406 |

### Technical definitions of density functional inputs

In the following, $\rho_\alpha$ and $\rho_\beta$ are the $\alpha$ and $\beta$ spin densities; the total spin density is $\rho$;

The gradients of the density enter through $$\begin{aligned} \sigma_{\alpha\alpha} &=& \nabla\rho_\alpha \cdot \nabla\rho_\alpha \; , \sigma_{\beta\beta} = \nabla\rho_\beta \cdot \nabla\rho_\beta\; , \sigma _{\alpha\beta}= \nabla\rho_\alpha \cdot \nabla\rho_\beta \; , \nonumber \\ \sigma &=& \sigma_{\alpha\alpha}+\sigma_{\beta\beta}+2\sigma _{\alpha\beta}. \\ \chi_\alpha&=& \frac{\sqrt{\sigma_{\alpha\alpha}}}{\rho_\alpha^{4/3}}\;, \chi_\beta = \frac{\sqrt{\sigma_{\beta\beta}}}{\rho_\beta^{4/3}}\;.\\ \upsilon_\alpha&=&\nabla^2\rho_\alpha \; , \upsilon_\beta=\nabla^2\rho_\beta \; , \upsilon=\upsilon_\alpha+\upsilon_\beta \;.\end{aligned}$$ Additionally, the kinetic energy density for a set of (Kohn-Sham) orbitals generating the density can be introduced through $$\begin{aligned} \tau_\alpha&=&\sum_i^\alpha \left|{\bf\nabla}\phi_i\right|^2 \; , \tau_\beta=\sum_i^\beta \left|{\bf\nabla}\phi_i\right|^2 \;, \tau=\tau_\alpha+\tau_\beta \;.\end{aligned}$$

All of the available functionals are of the general form $$\begin{aligned} F\left[\rho_s,\rho_{\bar{s}}, \sigma_{ss},\sigma_{\bar{s}\bar{s}},\sigma_{s\bar{s}}, \tau_s,\tau_{\bar{s}}, \upsilon_s,\upsilon_{\bar{s}} \right] &=& \int d^3{\bf r} K\left(\rho_s,\rho_{\bar{s}}, \sigma_{ss},\sigma_{\bar{s}\bar{s}},\sigma_{s\bar{s}}, \tau_s,\tau_{\bar{s}}, \upsilon_s,\upsilon_{\bar{s}} \right)\end{aligned}$$ where $\bar{s}$ is the conjugate spin to $s$.

### Lists of component density functionals

Below is a list of the primary correlation and exchange functionals supported by Molpro. These can be combined into total density functionals. For example, the following input will produce a B3LYP calculation:

`{dfunc,B88,DIRAC,LYP,VWN5,dftfac=[0.72,0.08,0.81,0.19],exfac=0.2}{df-rks}`

Note that many commonly used combinations of primary functionals are predefined as alias keywords (e.g., `{df-rks,b3lyp}`

is an equivalent input. See above).

Xalpha beta gamma. Divergence free semiempirical gradient-corrected exchange energy functional. $\lambda=\gamma$ in ref.`B86`

Xalpha beta gamma with Modified Gradient Correction. B86 with modified gradient correction for large density gradients.`B86MGC`

Xalpha beta gamma Re-optimised. Re-optimised $\beta$ of B86 used in part 3 of Becke’s 1997 paper.`B86R`

Becke 1988 Exchange Functional`B88`

Becke 1988 Correlation Functional. Correlation functional depending on B86MGC exchange functional with empirical atomic parameters, $t$ and $u$. The exchange functional that is used in conjunction with B88C should replace B88MGC here.`B88C`

Becke 1995 Correlation Functional. $\\tau$ dependent Dynamical correlation functional.`B95`

Density functional part of B97. This functional needs to be mixed with 0.1943*exact exchange.`B97DF`

Density functional part of B97 Re-parameterized by Hamprecht et al. Re-parameterization of the B97 functional in a self-consistent procedure by Hamprecht et al. This functional needs to be mixed with 0.21*exact exchange.`B97RDF`

Becke-Roussel Exchange Functional. A. D. Becke and M. R. Roussel,Phys. Rev. A`BR`

**39**, 3761 (1989) $$K=\frac{1}{2}\sum_s \rho_s U_s ,$$ where $$U_s=-(1-e^{-x}-xe^{-x}/2)/b ,$$ $$b=\frac{x^3e^{-x}}{8\pi\rho_s}$$ and $x$ is defined by the nonlinear equation $$\frac{xe^{-2x/3}}{x-2}=\frac{2\pi^{2/3}\rho_s^{5/3}}{3Q_s} ,$$ where $$Q_s=(\upsilon_s-2\gamma D_s)/6 ,$$ $$D_s=\tau_s-\frac{\sigma_{ss}}{4\rho_s}$$ and $$\gamma=1.$$Becke-Roussel Exchange Functional — Uniform Electron Gas Limit. A. D. Becke and M. R. Roussel,Phys. Rev. A`BRUEG`

**39**, 3761 (1989) As for`BR`

but with ${\gamma=0.8}$.Becke-Wigner Exchange-Correlation Functional. Hybrid exchange-correlation functional comprimising Becke’s 1998 exchange and Wigner’s spin-polarised correlation functionals.`BW`

Colle-Salvetti correlation functional. R. Colle and O. Salvetti, Theor. Chim. Acta 37, 329 (1974); C. Lee, W. Yang and R. G. Parr, Phys. Rev. B`CS1`

**37**, 785(1988)`CS1`

is formally identical to`CS2`

, except for a reformulation in which the terms involving $\upsilon$ are eliminated by integration by parts. This makes the functional more economical to evaluate. In the limit of exact quadrature,`CS1`

and`CS2`

are identical, but small numerical differences appear with finite integration grids.Colle-Salvetti correlation functional. R. Colle and O. Salvetti, Theor. Chim. Acta 37, 329 (1974); C. Lee, W. Yang and R. G. Parr, Phys. Rev. B`CS2`

**37**, 785(1988)`CS2`

is defined through $$\begin{aligned} K &=& -a \left({ \rho+2b\rho^{-5/3} \left[ \rho_\alpha t_{\alpha} + \rho_\beta t_{\beta} -\rho t_W \right] e^{-c\rho^{-1/3}} \over 1+d \rho^{-1/3} }\right) \end{aligned}$$ where $$\begin{aligned} t_{\alpha} &=&\frac{\tau_\alpha}{2}-\frac{\upsilon_\alpha}{8} \\ t_{\beta} &=&\frac{\tau_\beta}{2}-\frac{\upsilon_\beta}{8} \\ t_{W} &=& {1\over 8} {\sigma \over \rho} - {1\over 2} \upsilon \end{aligned}$$ and the constants are $a=0.04918, b=0.132, c=0.2533, d=0.349$.Slater-Dirac Exchange Energy. Automatically generated Slater-Dirac exchange.`DIRAC`

Short-range LDA correlation functional. Local-density approximation of correlation energy`ECERF`

for short-range interelectronic interaction ${\rm erf}(\mu r_{21})/r_{12}$,

S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys. Rev. B **73**, 155111 (2006). $$\begin{aligned}
& &
\epsilon_c^{\rm SR}(r_s,\zeta,\mu) = \epsilon_c^{\rm PW92}(r_s,\zeta)-
\nonumber \\
& & \frac{[\phi_2(\zeta)]^3Q\left(\frac{\mu\sqrt{r_s}}{\phi_2(\zeta)}\right)+a_1 \mu^3+a_2 \mu^4+
a_3\mu^5+a_4\mu^6+a_5\mu^8}{(1+b_0^2\mu^2)^4},\end{aligned}$$ where $$Q(x)=\frac{2\ln(2)-2}{\pi^2}\ln\left(\frac{1+a\,x+b\,x^2+c\,x^3}{1+a\,x+
d\,x^2}\right),$$ with $a=5.84605$, $c=3.91744$, $d=3.44851$, and $b=d-3\pi\alpha/(4\ln(2)-4)$. The parameters $a_i(r_s,\zeta)$ are given by $$\begin{aligned}
a_1 & = & 4 \,b_0^6 \,C_3+b_0^8 \,C_5, \nonumber \\
a_2 & = & 4 \,b_0^6 \,C_2+b_0^8\, C_4+6\, b_0^4 \epsilon_c^{\rm PW92}, \nonumber \\
a_3 & = & b_0^8 \,C_3, \nonumber \\
a_4 & = & b_0^8 \,C_2+4 \,b_0^6\, \epsilon_c^{\rm PW92} \nonumber, \\
a_5 & = & b_0^8\,\epsilon_c^{\rm PW92}, \nonumber\end{aligned}$$ with $$\begin{aligned}
C_2 & = & -\frac{3(1\!-\!\zeta^2)\,g_c(0,r_s,\zeta\!=\!0)}{8\,r_s^3} \nonumber \\
C_3 & = & - (1\!-\!\zeta^2)\frac{g(0,r_s,\zeta\!=\!0)}{\sqrt{2\pi}\, r_s^3}
\nonumber \\
C_4 & = & -\frac{9\, c_4(r_s,\zeta)}{64 r_s^3} \nonumber \\
C_5 & = & -\frac{9\, c_5(r_s,\zeta)}{40\sqrt{2 \pi} r_s^3}\nonumber \\
c_4(r_s,\zeta)
& = & \left(\frac{1\!+\!\zeta}{2}\right)^2g''\left(0,
r_s\left(\frac{2}{1\!+\!\zeta}\right)^{1/3}\!\!\!\!\!\!\!\!, \,\,\,\,\,\zeta\!=\!1\right)+ \left(\frac{1\!-\!\zeta}{2}\right)^2 \times \nonumber \\ & & g''\left(0,
r_s\left(\frac{2}{1\!-\!\zeta}\right)^{1/3}\!\!\!\!\!\!\!\!, \,\,\,\,\,\zeta\!=\!1\right)
+ (1\!-\!\zeta^2)D_2(r_s)-\frac{\phi_8(\zeta)}{5\,\alpha^2\,r_s^2}
\nonumber \\
c_5(r_s,\zeta)
& = & \left(\frac{1\!+\!\zeta}{2}\right)^2g''\left(0,
r_s\left(\frac{2}{1\!+\!\zeta}\right)^{1/3}\!\!\!\!\!\!\!\!, \,\,\,\,\,\zeta\!=\!1\right)+
\left(\frac{1\!-\!\zeta}{2}\right)^2 \times \nonumber \\ & & g''\left(0,
r_s\left(\frac{2}{1\!-\!\zeta}\right)^{1/3}\!\!\!\!\!\!\!\!, \,\,\,\,\,\zeta\!=\!1\right)+ (1\!-\!\zeta^2)D_3(r_s),\end{aligned}$$ and $$\begin{aligned}
\phantom{\bigl[} b_0(r_s) = 0.784949\,r_s \\
\phantom{\Biggl[} g''(0,r_s,\zeta\!=\!1) = \frac{2^{5/3}}{5\,\alpha^2 \,r_s^2} \,
\frac{1-0.02267 r_s}{\left(1+0.4319 r_s+0.04 r_s^2\right)} \\
\phantom{\Biggl[}D_2(r_s) = \frac{e^{- 0.547 r_s}}{r_s^2}\left(-0.388 r_s+0.676 r_s^2\right) \\
\phantom{\Biggl[}D_3(r_s) = \frac{e^{-0.31 r_s}}{r_s^3}\left(-4.95 r_s+ r_s^2\right).\end{aligned}$$ Finally, $\epsilon_c^{\rm PW92}(r_s,\zeta)$ is the Perdew-Wang parametrization of the correlation energy of the standard uniform electron gas [J.P. Perdew and Y. Wang, Phys. Rev. B **45**, 13244 (1992)], and $$g(0,r_s,\zeta\!=\!0)=\frac{1}{2}(1-Br_s+Cr_s^2+Dr_s^3+Er_s^4)\,{\rm e}^{-dr_s},$$ is the on-top pair-distribution function of the standard jellium model [P. Gori-Giorgi and J.P. Perdew, Phys. Rev. B **64**, 155102 (2001)], where $B=-0.0207$, $C=0.08193$, $D=-0.01277$, $E=0.001859$, $d=0.7524$. The correlation part of the on-top pair-distribution function is $g_c(0,r_s,\zeta\!=\!0)=g(0,r_s,\zeta\!=\!0)-\frac{1}{2}$.

Range-Separated Correlation Functional. Toulouse-Colonna-Savin range-separated correlation functional based on PBE, see J. Toulouse et al., J. Chem. Phys.`ECERFPBE`

**122**, 014110 (2005).

Exact Exchange Functional. Hartree-Fock exact exchange functional can be used to construct hybrid exchange-correlation functional.`EXACT`

Short-range LDA correlation functional. Local-density approximation of exchange energy`EXERF`

for short-range interelectronic interaction ${\rm erf}(\mu r_{12})/r_{12}$,

A. Savin, in Recent Developments and Applications of Modern Density Functional Theory, edited by J.M. Seminario (Elsevier, Amsterdam, 1996). $$\begin{aligned}
\epsilon_x^{\rm SR}(r_s,\zeta,\mu) & = & \frac{3}{4\pi}\frac{\phi_4(\zeta)}{\alpha\,r_s}- \frac{1}{2}(1\!+\!\zeta)^{4/3}
f_x\left(r_s,\mu(1\!+\!\zeta)^{-1/3}\right)+ \nonumber \\
& & \frac{1}{2}(1\!-\!\zeta)^{4/3}
f_x\left(r_s,\mu(1\!-\!\zeta)^{-1/3}\right)\end{aligned}$$ with $$\phi_n(\zeta)=\frac{1}{2}\left[ (1\!+\!\zeta)^{n/3}+(1\!-\!\zeta)^{n/3} \right],$$ $$\begin{aligned}
f_x(r_s,\mu) & = & -\frac{\mu}{\pi}\biggl[(2y-4y^3)\,e^{-1/4y^2}-
3y+4y^3+ \sqrt{\pi}\,{\rm erf}\left(\frac{1}{2y}\right)\biggr], \nonumber \\
\qquad y & = & \frac{\mu\,\alpha\,r_s}{2},\end{aligned}$$ and $\alpha=(4/9\pi)^{1/3}$.

Range-Separated Exchange Functional. Toulouse-Colonna-Savin range-separated exchange functional based on PBE, see J. Toulouse et al., J. Chem. Phys.`EXERFPBE`

**122**, 014110 (2005).

Gill’s 1996 Gradient Corrected Exchange Functional`G96`

Handy least squares fitted functional`HCTH120`

Handy least squares fitted functional`HCTH147`

Handy least squares fitted functional`HCTH93`

Meta GGA Correlation Functional. Henderson-Janesko-Scuseria range-separated exchange functional based on a model of an exchange hole derived by a constraint-satisfaction technique, see T. M. Henderson et al., J. Chem. Phys.`HJSWPBEX`

**128**, 194105 (2008).

Local tau Approximation. LSDA exchange functional with density represented as a function of $\tau$.`LTA`

Lee, Yang and Parr Correlation Functional. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B`LYP`

**37**, 785(1988); B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett.**157**, 200 (1989).

doi:10.1103/PhysRevB.37.785,10.1016/0009-2614(89)87234-3

M05-2X Meta-GGA Correlation Functional`M052XC`

M05-2X Meta-GGA Exchange Functional`M052XX`

M05 Meta-GGA Correlation Functional`M05C`

M05 Meta-GGA Exchange Functional`M05X`

M06-2X Meta-GGA Correlation Functional`M062XC`

M06-2X Meta-GGA Exchange Functional`M062XX`

M06 Meta-GGA Correlation Functional`M06C`

M06-HF Meta-GGA Correlation Functional`M06HFC`

M06-HF Meta-GGA Exchange Functional`M06HFX`

M06-L Meta-GGA Correlation Functional`M06LC`

M06-L Meta-GGA Exchange Functional`M06LX`

M06 Meta-GGA Exchange Functional`M06X`

Meta GGA Correlation Functional. Meta-GGA correlation functional based on first principles, see M. Modrzejewski et al., J. Chem. Phys.`M12C`

**137**, 204121 (2012).

Exchange Functional for Accurate Virtual Orbital Energies`MK00`

Exchange Functional for Accurate Virtual Orbital Energies. MK00 with gradient correction of the form of B88X but with different empirical parameter.`MK00B`

.. Gradient correction to VWN.`P86`

PBE Correlation Functional`PBEC`

doi:10.1103/PhysRevLett.77.3865

PBEsol Correlation Functional`PBESOLC`

doi:10.1103/PhysRevLett.100.136406

PBEsol Exchange Functional`PBESOLX`

doi:10.1103/PhysRevLett.100.136406

PBE Exchange Functional`PBEX`

doi:10.1103/PhysRevLett.77.3865

Revised PBE Exchange Functional. Changes the value of the constant R from the original PBEX functional`PBEXREV`

doi:10.1103/PhysRevLett.80.890

.. GGA Exchange Functional.`PW86`

Perdew-Wang 1991 GGA Correlation Functional`PW91C`

Perdew-Wang 1991 GGA Exchange Functional`PW91X`

Perdew-Wang 1992 GGA Correlation Functional. Electron-gas correlation energy.`PW92C`

Test for number of electrons`STEST`

Thomas-Fermi Kinetic Energy. Automatically generated Thomas-Fermi Kinetic Energy.`TFKE`

Tozer and Handy 1998. Density and gradient dependent first row exchange-correlation functional.`TH1`

.. Density and gradient dependent first row exchange-correlation functional.`TH2`

.. Density and gradient dependent first and second row exchange-correlation functional.`TH3`

.. Density an gradient dependent first and second row exchange-correlation functional.`TH4`

.. Density and gradient dependent first row exchange-correlation functional for closed shell systems. Total energies are improved by adding $DN$, where $N$ is the number of electrons and $D=0.1863$.`THGFC`

doi:10.1016/S0009-2614(97)00586-1

.. Density and gradient dependent first row exchange-correlation functional. FCFO = FC + open shell fitting.`THGFCFO`

doi:10.1016/S0009-2614(97)00586-1

.. Density and gradient dependent first row exchange-correlation functional.`THGFCO`

doi:10.1016/S0009-2614(97)00586-1

.. Density dependent first row exchange-correlation functional for closed shell systems.`THGFL`

doi:10.1016/S0009-2614(97)00586-1

TPSS Correlation Functional. J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett.`TPSSC`

**91**, 146401 (2003).

doi:10.1103/PhysRevLett.91.146401

TPSS Exchange Functional. J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett.`TPSSX`

**91**, 146401 (2003).

doi:10.1103/PhysRevLett.91.146401

.`VSXC`

von Weizsäcker kinetic energy. Automatically generated von Weizsäcker kinetic energy.`VW`

Vosko-Wilk-Nusair (1980) III local correlation energy. VWN 1980(III) functional`VWN3`

Vosko-Wilk-Nusair (1980) V local correlation energy. VWN 1980(V) functional. The fitting parameters for $\Delta\varepsilon_{c}(r_{s},\zeta)_{V}$ appear in the caption of table 7 in the reference.`VWN5`

M05-2X Meta-GGA Exchange-Correlation Functional. Here it means M05-2X exchange-correlation part which excludes HF exact exchange term. Y. Zhao, N. E. Schultz, and D. G. Truhlar, J. Chem. Theory Comput.`XC-M05-2X`

**2**, 364 (2006).

M05 Meta-GGA Exchange-Correlation Functional. Here it means M05 exchange-correlation part which excludes HF exact exchange term. Y. Zhao, N. E. Schultz, and D. G. Truhlar, J. Chem. Phys.`XC-M05`

**123**, 161103 (2005).

M06-2X Meta-GGA Exchange-Correlation Functional. Here it means M06-2X exchange-correlation part which excludes HF exact exchange term. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008).`XC-M06-2X`

M06-HF Meta-GGA Exchange-Correlation Functional. Here it means M06-HF exchange-correlation part which excludes HF exact exchange term. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A 110, 13126 (2006).`XC-M06-HF`

M06-L Meta-GGA Exchange-Correlation Functional. Y. Zhao and D. G. Truhlar, J. Chem. Phys.`XC-M06-L`

**125**, 194101 (2006).

M06 Meta-GGA Exchange-Correlation Functional. Here it means M06 exchange-correlation part which excludes HF exact exchange term. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008).`XC-M06`

M08-HX Meta-GGA Exchange-Correlation Functional. Here it means M08-HX exchange-correlation part which excludes HF exact exchange term. Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput.`XC-M08-HX`

**4**, 1849 (2008).

M08-SO Meta-GGA Exchange-Correlation Functional. Here it means M08-SO exchange-correlation part which excludes HF exact exchange term. Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput.`XC-M08-SO`

**4**, 1849 (2008).

M11-L Exchange-Correlation Functional. R. Peverati and D. G. Truhlar, Journal of Physical Chemistry Letters 3, 117 (2012).`XC-M11-L`

SOGGA Exchange-Correlation Functional. Y. Zhao and D. G. Truhlar, J. Chem. Phys.`XC-SOGGA`

**128**, 184109 (2008).

SOGGA11-X Exchange-Correlation Functional. Here it means SOGGA11-X exchange-correlation part which excludes HF exact exchange term. R. Peverati and D. G. Truhlar, J. Chem. Phys.`XC-SOGGA11-X`

**135**, 191102 (2011).

SOGGA11 Exchange-Correlation Functional. R. Peverati, Y. Zhao and D. G. Truhlar, J. Phys. Chem. Lett. 2 (16), 1991 (2011).`XC-SOGGA11`

### LibXC functionals

Functionals from the LibXC library https://www.tddft.org/programs/libxc can be used with the dft/ks programs by entering the functional names as given in the LibXC documentation. For example

ks,XC_GGA_X_B88,XC_GGA_C_LYP

defines an input for the B-LYP functional using the LibXC functional implementation instead of the internal Molpro ones. Correspondingly,

ks,XC_HYB_GGA_XC_B3LYP

defines a hybrid B3LYP calculation that would correspond to

ks,b88,dirac,lyp,vwn3; dftfac,0.72,0.08,0.81,0.19; exchange,0.2

using the corresponding individual functionals from Molpro. Note the LibXC's B3LYP does *not* correspond to Molpro's B3LYP, since the latter uses the VWN5 version of the VWN local correlation
functional instead of the VWN3 version.

Calculations with range-separated functionals require some more input in order to get integrals for the corresponding long-range erf interaction. For example

omega=0.3 !range-separation parameter srx=0.157706 !short-range exchange {int; ERFLERFC,mu=$omega,srfac=$srx} ks,XC_HYB_GGA_XC_WB97X

will have to be entered in the input file for performing calculations with the wB97X functional by Chai and Head-Gordon https://doi.org/10.1063/1.2834918.

The explicit input for the omega and short-range exchange scaling factor (`srx`

) can be avoided if the calculation is run in density-fitting mode, like:

df-ks,XC_HYB_GGA_XC_WB97X

For getting the corresponding wB97X-D functional https://doi.org/10.1039/B810189B, add the D2-type dispersion correction

df-ks,XC_HYB_GGA_XC_WB97X; disp,d2

see section Empirical damped dispersion correction for more details.

Another class of functionals are vdwDF's (van-der-Waals density functionals) which include also a nonlocal correlation functional part, see for instance the review https://doi.org/10.1088/0034-4885/78/6/066501. Molpro currently supports calculations with the VV10 (Vydrov & Van Voorhis 2010) functional https://doi.org/10.1063/1.3521275. The local part for this functional can be taken from LibXC:

ks,XC_GGA_XC_VV10

Regarding input descriptions for the nonlocal functional contribution, browse to section van-der-Waals DF's.

### Implementing new functionals

New functionals are implemented based upon the automatic code generation (ACG) program (doi:10.1016/S0010-4655(01)00148-5 ). In order to work the program requires the maple mathematics program and an XSLT parser, defined by the variable XSLT in CONFIG.

The format of the input file is an XML file containing all of the information about the new functional. All density functional XML files are placed in the directory `lib/df`

and are automatically activated on the next instance of the `make`

command in the Molpro base directory.

The root element of the XML document is `content`

. At the next level the element, `functional`

is expected, 1 per file.

The `functional`

element has an `id`

attribute which is used as the keyword for the functional in Molpro, and optional `doi`

attribute for specifying a reference. The allowed elements are defined in the following table:

Elements allowed for defining functionals | |
---|---|

title | Text to appear as a heading for the functional documentation |

tex | Text to document the functional |

The final element is `maple`

for which multiple cases are allowed. A typical maple expression such as

A:=1.2:

is written as

<maple lhs="A">1.2</maple>.

To input a Maple procedure such as

add_together:=proc(a,b) a+b end:

one should write

<maple lhs="add_together" proc="a,b">a+b</maple>.

As an example the Perdew-Wang 1991 GGA exchange functional is given below:

<?xml version="1.0" encoding="ISO-8859-1"?> <content> <functional id="PW91X"> <title>Perdew-Wang 1991 GGA Exchange Functional</title> <maple lhs="g">1/2*E(2*rho(s))</maple> <maple lhs="G">1/2*E(2*rho(s))</maple> <maple lhs="E" proc="n"> -3/(4*Pi)*(3*Pi^2)^(1/3)*n^(4/3)*F(S)</maple> <maple lhs="S">chi(s)/(2*(6*Pi^2)^(1/3))</maple> <maple lhs="F" proc="S"> (1+0.19645*S*arcsinh(7.7956*S) + (0.2743-0.1508*exp(-100*S^2))*S^2)/ (1+0.19645*S*arcsinh(7.7956*S)+0.004*S^4) </maple> </functional> </content>

### Implementing new hybrid-functionals

Hybrid functionals are defined in the file `lib/dalias.registry`

. Entries can be added to this file as required, and then run `make`

to update the other registry files.

## Exact exchange Kohn-Sham methods

A number of exact exchange Kohn-Sham methods are implemented in Molpro. Here, in contrast to standard hydrid functionals, the functional derivative of the exact exchange energy is taken with respect to the electron density:
\begin{equation}
v_{\text{x}}({\bf r})=\frac{\delta E_{\text{x}}[\rho]}{\delta\rho({\bf r})} \label{eq:vx}
\end{equation}
which yields a *local* exchange potential opposed to the Hartree-Fock method in which the potential is nonlocal. Due to this, while total energies from exact exchange KS methods are close to HF energies, orbitals and orbital energies can strongly differ from each other. The advantage of exact exchange KS methods is that the orbitals and orbital energies can directly be used as input data for time-dependent DFT methods employing local KS kernels for calculating response properties or excitation energies, see also section Time-dependent density functional theory. Moreover, orbitals from exact exchange KS calculations are also suitable for use in DFT-SAPT calculations (section symmetry-adapted intermolecular perturbation theory) because the local exact exchange potential $v_{\text{x}}({\bf r})$ possesses the correct asymptotic behaviour and thus yields an improved description of the electronic density in the asymptotic range.

Since the exact exchange energy is only an implicit functional of the density (through the direct dependency on the occupied molecular orbitals) the functional derivative in Eq. \eqref{eq:vx} can not be obtained by a direct differentiation. Two different methods are implemented in Molpro by which the local exchange potential can be computed, the LHF (Local Hartree-Fock) method by Della Sala and Görling [1] and the Optimised Effective Potential (OEP) method employing the auxiliary basis set implementation of Ref. [5]. The two different approaches are described in the following subsections.

Note that (almost) all methods described in this section can **only be run without using symmetry**. Furthermore, for most of the approaches also no open-shell implementation is available.

**Bibliography**:

$[1]$ F. Della Sala and A. Görling, J. Chem. Phys. **115**, 5718 (2001)

$[2]$ A. Heßelmann and F.R. Manby, J. Chem. Phys. **123**, 164116 (2005)

$[3]$ A. Görling, A. Heßelmann, M. Jones and M. Levy, J. Chem. Phys. **128**, 104104 (2008)

$[4]$ A. Heßelmann and A. Görling, Chem. Phys. Lett. **455**, 110 (2008)

$[5]$ A. Heßelmann, A. Götz, F. Della Sala and A. Görling, J. Chem. Phys. **127**, 054102 (2007)

$[6]$ A. Heßelmann, J. Chem. Theory Comput. **14**, 1943 (2018)

### Local Hartree-Fock (LHF) method

The LHF method is an *approximate* exact exchange KS method which derives the local exchange potential by starting from the assumption that the HF and the exchange-only KS determinants are identical, see Ref. [1]. An advantage of the method is that the approach always yields smooth and physically correct potentials (cf. the next section and Refs. [3,4]). A disadvantage of the LHF method is that the potential has to be computed on a numerical grid which can be quite costly, particularly due to the fact that the Slater potential contribution to $v_{\text{x}}$ requires a computation of the electrostatic potential of the product of two AO basis functions on the grid [1]. The use of the density-fitting approach of Ref. [2], however, significantly reduces this cost, see below.

A standard LHF calculation with Molpro can be performed with the Kohn-Sham DFT program using

ks,lhf

During the SCF iterations, the program prints out the HOMO eigenvalue and the Slater exchange energy which should be identical to the exact exchange energy. It is always good to check in the output that the latter holds true.

Instead of performing the LHF calculation self-consistently, it is also possible to do a one-step LHF calculation using orbitals and orbitalk energies from a different HF/KS calculation as input. This can be done by calling the LHF program directly, e.g.

{hf; save,2100.2} {lhf,orb=2100.2}

The following options can be used with the LHF program:

orbital record containing previous orbitals (mandatory)`ORB`

use density fitting to compute the exchange matrix`LHFFIT`

record in which to write the local exchange potential`POT`

recompute the LHF orbitals and orbital energies from an already computed exchange potential`RESTART`

if set to 2, compute the Slater potential with the method described in appendix A in Ref. [1] (default: 1)`SLA`

batch size for numerical grid`BLOCK`

if set to 1, the Slater potential and the response potential will be written to the files ’slater.dat’ and ’response.dat’ (names can be changed by setting the variables SLATER and RESPONSE in the Molpro input file before calling the LHF program). The ’response.dat’ file will also contain the total exchange potential. The default is to not punch out the potentials.`PUNCH`

scale the response potential by $\rho/(\rho+\mathrm{damp})$ in order to force the potential to decay more rapidly in the asymptotic range (typical value $\mathrm{damp}=1d-5$). Not used by default, but can be helpful in cases where the convergence of the KS calculation is slow.`DAMP`

To significantly speed up the computation of the Slater potential contribution, the density-fitting method of Ref. [2] can be used. It is automatically switched to if a density-fitted KS calculation is performed:

df-ks,lhf

Note, however, that now an additional auxiliary basis set with the name DFLHF needs to be added in the basis input section which, e.g.:

... set,dflhf default,avtz/jkfit ...

Notice that the use of the underlying JKFit basis sets that correspond to the given orbital basis set are most suitable, because the fitting is performed for products of two occupied molecular orbitals. Moreover, note that **no contracted auxiliary basis sets** can be used with the DF-LHF program.

There also exists a one-step density-fitting analogue to the LHF program:

{hf; save,2100.2} {dflhf,orb=2100.2}

In addition to the options for the LHF program given above, the DFLHF program has the further options

set to $\ne 0$ for using a mixed Poisson/Gaussian auxiliary basis set to further speed up the calculation of the Slater potential, see Ref. [2]. The Poisson auxiliary basis set needs to be supplied in the basis input section.`POISSON`

switch 3-index integral routine to be used (can take the values 1,2 or 3 (default is 2))`INT`

It is also possible to reuse the exchange potential computed by the LHF/DFLHF programs in a subsequent KS calculation. In this case, the potential is kept fixed during the SCF and only the Coulomb potential is variationally optimised. This can be done by, e.g. (’tlhf’ denoting ’Transformation Local Hartree-Fock’)

{hf; save,2100.2} {lhf,orb=2100.2} {ks,tlhf}

The LHF functional can also be combined with other functionals as in common hybrid-DFT methods. For example, the hybrid PBE0 functional (inluding 25 percent exact exchange) would be implemented by

{ks,lhf,pbex,pbec; dftfac,0.25,0.75,1.0}

### Optimised Effective Potential (OEP) method

The functional derivative of the exchange energy can also be calculated by using the optimised effective potential method (OEP) which was originally derived by Talman and Shadwick ( Phys. Rev. A **14**, 36 (1976)). This method sometines is also denoted as EXX (Exact Exchange) method if used in conjunction with exchange-only KS methods.

Unfortunately, OEP methods can be numerically unstable if finite basis sets are used (which is, of course, generally the case), see Refs. [3] and [4]. The reason for this originates from the fact that the OEP equation is a response-type equation where the underlying (static) KS response function is described in terms of products of occupied and unoccupied orbitals. It turns out that it generally is difficult to adequately represent the (potentially linearly dependent) occ-virt space by an auxiliary function space [3,4]. To be able to obtain physically sound OEP (exchange) potentials, it is necessary to either employ regularisation techniques to eliminate the small eigenvalues of the response matrix (see, e.g., T. Heaton-Burgess, F. A. Bulat, W. Yang, Phys. Rev. Lett. **98**, 256401 (2007)), or to use auxiliary basis sets that are *balanced* to the orbital basis set. As a rule of thumb, the attribute ’balanced’ here means that the potential does not change anymore when the orbital basis space is increased while the auxiliary basis set is kept fixed [4,5].

Balanced orbital-auxiliary basis sets have been derived for a limited number of elements in Ref. [5]. They can be specified in the basis input section by

set,orbital default,avtz-orb-oep set,oep default,avtz-aux-oep

(see the oepbas.libmol file in the library directory for all available OEP basis sets).

Auxiliary basis set OEP/EXX calculations based on the approach of Ref. [5] can be performed with the EXX program, which, using the most essential options, can be called as

{exx,orb=<orb.rec>,homo=1,charge=1}

where the supplement of previous orbitals in a given record is mandatory. Furthermore, the options `homo=1`

and `charge=1`

activate the HOMO and the charge constraints to the potential which adds to the numerical stability of the OEP method, see Ref. [5].

Note that the OEP orbital basis sets are uncontracted and therefore calculations for larger molecules can become quite expensive. Meaningful OEP calculations can, however, also be performed with standard contracted basis sets if some regularisation techniques are used in the solution to the OEP equation. For this, it is recommended to use

{exx,...,oep=7,gamma=3d-4}

where $\gamma$ is a regularisation factor used in a Tikhonov regularisation approach to the solution of the OEP equation, see Ref. [6].

Further options to the EXX program are:

maximum number of iterations`MAXIT`

threshold for energy convergence`SCFTHR`

choose between 7 (!) different OEP routines (an explanation of all these would blow up this section significantly, default is OEP=1)`OEP`

debug print level`PRINT`

use Fermi-Amaldi potential (instead of OEP potential)`FERAMA`

if set to $\ne$ 0, use density fitting for two-electron integrals`DFIT`

set to $\ne$ 0 for integral-direct calculation (not very efficient at the moment)`DIRECT`

the reference potential can be calculated with the ELP method using OEP=4 or OEP=5. With the setting of ITVREF to a value $n$, the refernce potential can be held fixed from iteration $n+1$ onwards (to speed up the calculation).`ITVREF`

The underlying OEP module, that is called inside the EXX program, has the further options:

can be set to ’BASIS’ or ’BASIS2’ according to the two schemes described in Ref. [5] to determine the coefficients for the reference potential`VREF`

can be set to zero to omit calculation of the rest potential`VREST`

norm (should be ’J’)`AUX`

set (x,y,z) coordinates for a starting position used for plotting the potential`ATOM1`

set (x,y,z) coordinates for a terminal position used for plotting the potential`ATOM2`

do not use!`SYMCEN1/2`

can be used to switch on several tests during the generation of the OEP potential`TEST`

enter a name of the file to which the potential is written to (along a line that can be specified with ATOM1/2, see above)`POT`

set to $\ne$ 0 for activating the charge condition`CHARGE`

set to $\ne$ 0 for activating the HOMO condition`HOMO`

can additionally be used to fix the orbital energy of an orbital whose energy lies below the HOMO level (is not available for all OEP routines)`HOMO1`

set to $\ne$ 0 to print out the eigenvalues of the response matrix (useful to test stability of the OEP method)`RESPEIG`

set to $\ne$ 0 for printing the right-hand side in the output`RHS`

set to $\ne$ 0 for printing the response matrix`XPRINT`

matrix inversion methods to solve the OEP equation. The different options are: GESV (default), SVD, SVD2, TIKH`SOLVE`

threshold value for singular value decomposition (SVD) or Tikhonov regularisation (TIKH)`THR`

an integer number can be entered to filter out the response function eigenvectors with the lowest eigenvalues in an SVD matrix inversion`FILTER`

here, SVD is used to perform a singular value decomposition of the (aux$|$occ$\times$virt) 3-index integrals`SVD`

the corresponding SVD threshold for this`SVDTHR`

set to $\ne$ 0 to use contraction of auxiliary basis through the SVD of (aux$|$occ$\times$virt) integrals`CONTR`

solve OEP within the AO basis set making use of the Laplace transformation`LAPLACE`

dump the nonlocal and local exchange potentials to a record`VXDUMP`

dump record for OEP which contains the orbitals, the density and the eigenvalues`OEPDUMP`

record in which the coefficients of the potential can be written to`SAVE`

record from which the coefficients can be read`READ`

debug print level`OEPPRINT`

specify an orbital dump record from a previous calculation used to determine the vector for the HOMO condition (useful if the HOMO swaps with a lower lying level with a different method)`HOMOORB`

can be used to specify the number of an orbital whose eigenvalue should be fixed (see also HOMO1 option)`FIXEIG`

set the norm for the OEP equation (not available for all OEP routines)`EQNORM`

switch the mode in ELP calculations`ELPMODE`

can be used to deallocate some arrays that are used during the OEP calculation`EXIT`

threshold for testing the energy convergence. When the energy change is below this threshold, the exchange potential will be kept fixed and will not be calculated anymore.`CONV`

specify a record number for a record in which to write the matrix $\langle i|v_x^{\text{NL}}-v_x|j\rangle$ containing the difference of the nonlocal and the local exact exchange potentials. This is required in TDDFT calculations employing the exact-exchange KS kernel.`VXDIFF`

an eigenvalue threshold used in the decomposition of the overlap matrix, used to reduce the auxiliary function space (not available in all OEP routines)`THRS`

Tikhonov regularisation factor (only used for OEP=7 case). A value of $\gamma=3\cdot 10^{-4}$ should conventionally be used.`GAMMA`

Like in case of the LHF method, there also exists a one-step OEP program which takes orbitals from a preceeding SCF calculations, solves the OEP equations in this basis, and writes out the new orbitals and eigenvalues to a new record. A typical example that additionally also writes out the matrix elements $\langle i|v_x^{\text{NL}}-v_x|j\rangle$ to a record is given by:

{hf; save,2100.2} {oep,orb=2100.2,vxdiff=5000.2,gamma=3d-4,homo=1}

It is possible to do hybrid EXX-GGA calculations by using the `EXXHYXB`

program. In addition to the options given above for the `EXX`

program it can read the directives for standard KS calculations (see section directives). Further options for the `EXXHYXB`

program are:

factor for exact exchange`XFAX`

set to $\ne$ 0 in order to add the GGA contribution to the xc matrix to the rhs of the OEP equation`ADD`

if set to a value $\ne$ 0 the asymptotic correction method for the xc potential described in Ref. [6] will be used`SHIFT`

As an example, for using the local PBE0 functional employing the asymptotic correction method of Ref. [6] (LPBE0AC xc potential) the following input may be used:

{ks,lda; save,2100.2} {exxhyb,orb=2100.2,oep=7,homo=1,charge=1,add=1,xfac=0.25,shift=1 func,pbex,pbec; dftfac,0.75,1.0}

Spin-unrestriced calculations with the EXX or the hybrid-EXX methods can be done by using the `UEXX`

and the `UEXXHYXB`

programs. Here the initial orbitals to be supplied have to stem from a preceeding unrestricted HF or KS calculation.

## van-der-Waals DF's

Long-range correlation interactions can be described by so-called van-der-Waal's density functionals (vdwDF's) that are designed to correct the qualitatively wrong asymptotic behaviour of the xc energy between two noncovalently bonded systems. vdwDF's are clearly different to standard density functionals in that they are defined in terms of a nonlocal integral kernel $\phi({\bf r}_1,{\bf r}_2)$: $E_c^{vdwDF}=\int{\bf r}_1\int{\bf r_2}\rho({\bf r}_1) \phi({\bf r}_1,{\bf r}_2) \rho({\bf r}_2) $ where $\phi$ itself is a (nonlocal) function of the density $\rho$ and its derivative(s) $\nabla\rho$. Due to this, the integral can normally not be calculated analytically and requires a double-space integration using numerical quadrature techniques.

vdwDF correlation energies can be calculated both non-selfconsistently, calling the `vdwDF`

program, and selfconsistently
using the `ks`

program. An non-selfconsistent calculation for the VV10 functional by Vydrov & Van Voorhis https://doi.org/10.1063/1.3521275 would
look like

{ks,XC_GGA_XC_VV10; save,2100.2} eloc=energy {vdwDF; vv10,den=2100.2,toten=1} evv10=energy

where variable `eloc`

stores the total energy for the standard Kohn-Sham calculation with the local part to the functional (`XC_GGA_XC_VV10`

) and
`evv10`

stores the sum of `eloc`

plus the energy contribution for the nonlocal correlation part (`vv10`

). Note that the VV10 functional is currently the only available
option to vdwDF. It can take the following suboptions:

Dump record containing the density matrix (must be supplied)`DEN`

If`TOTEN`

`TOTEN>0`

the total energy will be printed (and stored in the intrinsic`energy`

variable) reading in the previous SCF energyThis can be used to alter the threshold for the grid accuracy. It is particularly useful for larger systems where the double-space integration becomes more expensive and where the integration (probably) need not be highly accurate. Default:`GTHR`

`GTHR=<global grid threshold set in dft/ks>`

.Block size for grid batches. Default:`NBLOCK`

`NBLOCK=0`

which actually means that no grid batching is performed at all (requiring more memory consumption).The`B`

`B`

parameter from the VV10 functional that is functional dependent. The program can detect a few number of functionals automatically and sets the value of`B`

to the optimised values determined by Hujo and Grimme https://doi.org/10.1039/C1CP20591A. The default value is`B=5.9`

that should be used in conjunction with the local`XC_GGA_XC_VV10`

functional part.

To use the complete VV10 functional in a self-consistent KS calculation, use:

ks,XC_GGA_XC_VV10,vv10nl

i.e., append the nonlocal functional keyword `vv10nl`

to the functional list. Further functional specific input can be added via, e.g.

ks,XC_GGA_XC_VV10,vv10nl; vdwdf,print=1,gthr=1d-8,b=5.9

see the corresponding options to vdwDF above. The only new option in SCF calculations is a print option that can be used to swich on `print=1`

and
off `print=0`

the printing of the VV10NL correlation energies during the SCF. By default this is activated `print>0`

.

## Empirical damped dispersion correction

Empirical damped dispersion corrections can be calculated in addition to Kohn-Sham calculations. This is particularly important in cases where long-range correlation effects become dominant.

The dispersion correction can be added to the DFT energy by using

ks,<func>; disp,<x>

or

ks,<func>,disp=<x>

The total energy will then be calculated as $$E_{\text{DFT-D2}}=E_{\text{DFT}}+E_{\text{disp}}(\mbox{D2})$$ if $x=d2$ (see Ref. [2]), $$E_{\text{DFT-D3}}=E_{\text{DFT}}+E_{\text{disp}}(\mbox{D3})$$ if $x=d3$ (see Ref. [3]), and $$E_{\text{DFT-D4}}=E_{\text{DFT}}+E_{\text{disp}}(\mbox{D4})$$ if $x=d4$ (see Ref. [4]).

Currently the default dispersion correction added to the DFT energy is the D4 dispersion correction developed by Grimme *et al.*, see Ref. [4].

The `disp,d3`

keyword can have the following additional options:

Functional name (default:`FUNC`

`FUNC='pbe'`

).Can have values 2 and 3 according to parametrisations from Refs. [2] and [3] (default:`VERSION`

`VERSION=3`

)Performs a detailed analysis of pair contributions.`ANAL`

Cartesian gradients are computed.`GRAD`

Calculate three-body dispersion contribution if`ABC`

`ABC`

$\ne$0.Do not calculate three-body dispersion contribution if`NOABC`

`NOABC`

$\ne$0.Calculate gradient numerically.`NUM`

Same as`OLD`

`VERSION=2`

.DFT-D3 original zero-damping.`ZERO`

DFT-D3 with Becke-Johnson finite-damping.`BJ`

Revised DFT-D3 with Becke-Johnson damping.`BJM`

Neglect threshold in Bohr for CN (default=40).`CNTHR`

Neglect threshold in Bohr for $E_{\text{disp}}$, (default=95).`CUTOFF`

Use special parameters for calculations with triple-zeta basis sets. Preliminary results in the SI of Ref. [3] indicate that results are slightly worse than with the default parameters and QZVP type basis sets. This option should be carfully tested for future use in very large computations.`TZ`

(see also 'https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dft-d3 for further documentation).

The DFT-D4 program (`disp,d4`

) accepts the options:

Molecular charge.`CHRG`

Scaling factor for leading order term.`S6`

Scaling factor for $C_8$ dependent term.`S8`

Scaling factor for $C_9$ dependent term.`S9`

Parameter in Becke-Johnson damping function.`A1`

Parameter in Becke-Johnson damping function.`A2`

Name of density functional (useful if the name conventions for functionals in the D4 program are different from Molpro's).`DFA`

Treatment of many-body dispersion.`LMBD`

`LMBD=0`

: skip calculation of many-body dispersion interactions,`LMBD=1`

: full RPA-like MBD,`LMBD=2`

: Axilrod-Teller-Muto three-body term,`LMBD=3`

: D3-like approximated ATM termsame as`NOMB`

`LMBD=0`

Calculate gradient if`GRAD`

`LMBD=0`

$\ne 0$.Calculate hessian (currently not used in Molpro).`HESS`

Set print level (always highest is set automatically with dispcorr4 command).`PRINT`

Charge scale height.`G_A`

Charge scale steepness.`G_C`

Kind of charge model.`REFQ`

(see also 'https://github.com/dftd4/dftd4 for further documentation).

Alternatively, the D2, D3 or D4 dispersion corrections can also be calculated separately from the DFT calculation using the following template:

ks,<func> eks=energy <dispcorr3,dispcorr4> eks_plus_disp=eks+edisp

using `dispcorr3`

for DFT-D3 or `dispcorr4`

for DFT-D4.

The older DFT-D1 [1] or DFT-D2 [2] dispersion corrections by Grimme can be calculated by using the input

ks,<func> eks=energy dispcorr eks_plus_disp=eks+edisp

with the following options to `dispcorr`

:

Adjusts the parametrisation used:`MODE`

`MODE=1`

uses parameters from Ref. [1] and`MODE=2`

uses parameters from Ref. [2] (default:`MODE=1`

)Overall scaling parameter $s_6$ (see Eq. \eqref{eq:edisp}) and Refs. [2,3] for optimised values).`SCALE,S6`

Damping function parameter (see Eq. \eqref{eq:fdamp}). Smaller values lead to larger corrections for intermediate distances.`ALPHA`

Choose damping function 'GRIMME' or 'CHAI' (default:'GRIMME'). With`DAMP`

`DAMP='CHAI`

' the damping function by Chai and Head-Gordon can be activated that is used, e.g., in the wB97XD functional, see Ref.[5]Parameter`A`

`a`

in the damping function of Chai et al., see Ref. [5]Calculate gradient if`GRAD`

`GRAD>0`

In the DFT-D1 and DFT-D2 method the dispersion energy is calculated as \begin{equation}\label{eq:edisp} E_{\text{disp}}=-s_6\sum_{i,j<i}^{N_{\text{atoms}}}f_{\text{damp}}(R_{ij}) \frac{C_6^{ij}}{R_{ij}^6}~. \end{equation} where $N_{\text{atoms}}$ is the total number of atoms, $R_{ij}$ is the interatomic distance of atoms $i$ and $j$, $s_6$ is a global scaling parameter depending on the choice of the functional used and the $C_6^{ij}$ values are calculated from atomic dispersion coefficients $C_6^i$ and $C_6^j$ in the following way: $$\begin{aligned} C_6^{ij}&=&2\frac{C_6^i C_6^j}{C_6^i + C_6^j} \qquad \text{(Ref.\ [1])}\\ C_6^{ij}&=&\sqrt{C_6^i C_6^j} \qquad \text{(Ref.\ [2])}\end{aligned}$$ The function $f_{\text{damp}}$ damps the dispersion correction for shorter interatomic distances and is given by: \begin{equation}\label{eq:fdamp} f_{\text{damp}}(R_{ij})=\frac{1}{1+\exp{[-\alpha(R_{ij}/ (R_{\text{vdW}}^i+R_{\text{vdW}}^j)-1)]}} \end{equation} with $R_{\text{vdW}}^i$ being the van-der-Waals radius for atom $i$ and $\alpha$ is a parameter that is usually set to 23 (Ref. [1]) or 20 (Ref. [2]).

Major changes in the D3 and D4 dispersion corrections compared to D2 are:

- The hybridisation states of the atoms are taken into account in the evaluation of the dispersion coefficients [3].
- The dispersion coefficients are dependent on local atomic charges obtained from a self-consistent tight-binding calculation [4]. Many-body interactions can be taken into account additionally.

Gradient contributions from the D3 and D4 dispersion correction are automatically computed in DFT geometry optimisations.

References:

$[1]$ S. Grimme, J. Comp. Chem. **25**, 1463 (2004).

$[2]$ S. Grimme, J. Comp. Chem. **27**, 1787 (2006).

$[3]$ S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys. **132**, 154104 (2010)

$[4]$ E. Caldeweyher, C. Bannwarth, and S. Grimme J. Chem. Phys. **147**, 034112 (2017)

$[5]$ J.-D. Chai, M. Head-Gordon Phys. Chem. Chem. Phys. **10**, 6615 (2008)

## Nonlocal DFT (NLDFT)

The deficiency of standard DFT functionals to describe long-range correlation energies (a.k.a. dispersion interactions) can be corrected through a double-Hirshfeld partitioning of the correlation functional energy density, using different weight functions $w_A$ and $\widetilde{w}_A$ for the local ($A=B$) and nonlocal ($A\ne B$) terms ($A,B$: atom indices):

\begin{align} E_{\text{c}}^{\text{NLDFT}}[\rho]=\sum_{A}^{N}\int d{\bf r} \, w_{A}^{2}({\bf r})\, z(\rho,\nabla\rho,\dots)+2\sum_{A}^{N} \sum_{B<A}^{N}\int d{\bf r}\,\widetilde{w}_{A}({\bf r}) \widetilde{w}_{B}({\bf r})\, z(\rho,\nabla\rho,\dots)\label{eq:EcNLDFT} \end{align}

The Hirshfeld weights $w_A$ of the local part of the correlation functional (first term on the rhs in Eq. \eqref{eq:EcNLDFT} are determined by using analytical expressions for the atomic densities ($\rho_A$) using Slaters rules. In contrast to this, the weights for the nonlocal part of the correlation energy (second term on the rhs in Eq. \eqref{eq:EcNLDFT} are determined by modified atomic densities using

\begin{align} \widetilde{\rho}_{A}({\bf r}) & = & \mathcal{N}\left(\rho_{A}({\bf r})+ f(r,p_{1},p_{2})\frac{p_{3}}{r^{6}}\right)\nonumber \\ & = & \mathcal{N}\left(\rho_{A}({\bf r})+\Big(\frac{1}{2} \Big[\mathrm{erf}(p_{1}(r+p_{2}))+\mathrm{erf}(p_{1}(r-p_{2}))\Big]\Big)^{6} \frac{p_{3}}{r^{6}}\right)\label{eq:RhoTilde} \end{align}

with $\mathcal{N}$: normalisation constant, $r$: distance from nucleus and $p_1,p_2,p_3$: parameters (see Ref. [1] for details). Due to the $1/r^{6}$ long-range tail of the modified densities $\widetilde{\rho}_{A}$, the atom-atom contribution to $E_{\text{c}}^{\text{NLDFT}}$ in Eq. \eqref{eq:EcNLDFT} adapts to the correct long-range behaviour of the correlation energy so that the NLDFT functional, unlike standard DFT functionals, can describe dispersion interactions, see Ref. [1]. Note that compared to the explicit damped dispersion corrections described in section empirical damped dispersion correction the NLDFT functional of Eq. \eqref{eq:EcNLDFT} also describes changes in the exchange-correlation potential due to the modification of the correlation energy density.

Since the densities $\widetilde{\rho}_{A}$ defined in Eq. \eqref{eq:RhoTilde} contain atomic parameters, calculations with the NLDFT functionals should be carried out only for systems containing atoms which were already parametrised, namely H, C, N, O, F, Si, P, S and Cl, see table 1 in Ref. [1]. Moreover, the NLDFT functional should (and can) only be used in conjunction with the PBE correlation functional. For the exchange functional part, the revised PBE exchange functional (with $\kappa=0.9$) from Ref. [1] should be used. NLDFT calculations can then be done by adding the keyword `nldft`

to the input for the Kohn-Sham calculation:

ks,pbec,revpbex; nldft

Note in this case that it is important to name the `pbec`

functional first in the list of functionals. This is necessary, because the xc potential contributions from various functionals are accumulated during the construction of the xc matrix while the NLDFT double-Hirshfeld decomposition is to be done for the correlation energy functional only, see Eq. \eqref{eq:EcNLDFT}.

At the end of a NLDFT calculation the atom-atom decomposition of the correlation energy as well as the sum of the local and nonlocal correlation energy terms are printed. Furthermore, the sum of atom-atom contributions for which $R_{AB}<3.0$ and $R_{AB}<6.0$ bohr are displayed.

All publications resulting from the use of this method should acknowledge Ref. [1] which describes the functional and analyses its performance for a wide range of different benchmark data bases. In addition, Ref. [2], which shows the performance of the NLDFT functional for describing the binding energies of large supramolecular complexes, may be cited. The performance of the NLDFT functional for geometry optimisations is investigated in Ref. [3].

References:

$[1]$ A. Heßelmann, J. Chem. Theory Comput. **9**, 273 (2013)

$[2]$ A. Heßelmann and T. Korona, J. Chem. Phys. **141**, 094107 (2014)

$[3]$ A. Heßelmann, J. Chem. Phys. **149**, 044103 (2018)

## Random-phase approximation

The random-phase approximation program (`rpatddft`

) can be used to calculate RPA correlation energies after a SCF calculation. Additionnally, it can be used to calculate dynamic dipole polarizabilities, C$_6$ dispersion coefficients, and excitation energies. The program currently works without point-group symmetry.

List of the main keywords:

`ECORR, <list of methods>`

Calculation of RPA correlation energies [1] (see options below)

`PROPERTIES, <list of methods>`

Calculation of dynamic dipole polarizabilities and C$_6$ dispersion coefficients [10] (see options below)

`EXCIT, <list of methods>`

Calculation of excitation energies [11] (see options below)

Record for input orbitals (required).`ORB,<orbrec>`

as well as contextual options (see later for an explanation on this):

`INTAC,NLAMBDA=<n>[,LAMBDA=<lambda>,WEIGHT=<weight>]`

Number of quadrature points for the Gauss-Legendre numerical integration along the adiabatic connection for RPA calculations (default is 7). If `LAMBDA`

and `WEIGHT`

are given, assumes a one-point quadrature with given abscissa and weight.

`INTFREQ,METHFREQ=<methfrq>,NFREQ=<nfreq>`

Options for the numerical integration over the frequency variable of RPA calculations. `METHFREQ`

governs the type of quadrature used (0(default) is Gauss-Chebyshev, 1 and 2 are Gauss-Legendre, 3 is Clenshaw-Curtis) and `NFREQ`

governs the number of quadrature points (default is 16).

Use the solid-state variant when performing DIEL calculations.`DIELMODE,mode`

and global options, shared by all commands inside the `rpatddft`

block:

Check matrices stability conditions in RPA calculations. When used without an`STAB`

`ECORR`

,`EXCIT`

or`PROPERTIES`

keyword, check the Hessian and RPA matrices eigenvalues and do nothing more.`DFTKERNEL,<funcx>[,<funcc>]`

Specify the exchange and correlation kernel for `EXCIT`

(if only one argument is given, it is understood to be the exchange-correlation kernel).

Computes C6 coefficients from last two saved polarizabilities.`C6`

Tamm-Dancoff approximation for`TDA`

`EXCIT`

and`PROPERTIES`

.The MP2 energy is calculated in certain situations where it is available almost for free, provided that some matrices are allocated. This behavior can be switched off by this`NOMP2`

`NOMP2`

keyword.For spin-unrestricted calculations, use a formalism where matrices are of $\alpha\alpha+\alpha\beta+\beta\alpha+\beta\beta$ dimensions (the default is to use a formalism with a nospinflip/splinflip block structure)`NOSPINBLOCK`

Exclude spin-flip dimensions of unrestricted RPA calculations that use the`NOSPINFLIP`

`NOSPINBLOCK`

formalism (not suitable for all RPA variants).Write files with eigenvalues, virtual orbital energies, dipole moments, dipole velocities, dipole accelerations and amplitudes from a TDA calculation`WRITEFILE`

#to comment#`VIAXPY`

Specify the two-electron integral transformation routine: 0 (still the default for spin-unrestricted and gradient calculations) is the ‘old’ one, 1 is the ‘old’ one that has been cleaned up, and 2 (default otherwise) is a much more efficient transformation using Molpro’s`INTEGRAL,<nbr>`

`transform`

routine.Explicitly specify the number of occupied orbitals (useful for fake pseudopotential calculations).`OCC,<nocc>`

Specify core orbitals (default: last specified core orbitals or, if none, atomic inner shells)`CORE,<core>`

Level of print expected from the output (from 0(default) to 3).`PRINT,<nbr>`

Level of print of integrals (AO,MO,Orbtials,…), from 0 to 4.`PRINT_INT,<nbr>`

If greater or equal to 1, will print out information on time spent in routines.`PRINT_TIME,<nbr>`

Calculation of RPA correlation energies `ECORR, <list of methods>`

If no method is given, a `SO2-RCCD`

calculation will be done (see below).

There are two main RPA *variants* [1]: dRPA (direct RPA, without the inclusion of an Hartree-Fock exchange kernel in the response function) and RPAx (with the Hartree-Fock exchange kernel included in the response function).

There are four main *formulations* in which the RPA equations can be derived. The adiabatic-connection fluctuation-dissipation theorem (ACFDT) equation involves integrations both over frequency and coupling constant: an analytical integration along the frequency variable followed by a quadrature on the coupling constant yields the adiabatic connection formulation (AC) [1], while an analytic integration on the coupling constant followed by a numerical integration over the frequency yields the dielectric formulation (DIEL) [2]. Two other formulations are obtained when the two integrations are carried out analytically: the plasmon formula (PLASMON) [1] and the ring coupled cluster doubles formulation (RCCD) [3]. These four formulations are not in general equivalent.

Most variants+formulations can readily by used in a spin-unrestricted context [6]. This is implemented in the code and does not need any further input from the user: the RPA program recognizes the spin-unrestricted character of a SCF calculation that was done beforehand and acts accordingly.

Gradients of most of the RCCD-formulation RPA energies are available, both without range-separation with RHF orbitals and with range-separation with RSH orbitals [9]. The calculations are triggered by the presence of the keyword `FORCE`

or `OPTG`

*after* the energy-related section (see examples at the end of the section).

The user can test the RPA program using `make rpatddfttest`

, which proposes a variety of tests for RPA correlation energy calculations.

The keywords for the methods are constructed on the model:

`<variant>-<formulation>-<alternative>`

For the AC formulation, the available methods are:

dRPA calculation (see Refs. [1] and [4]).`DRPAI-AC`

dRPA calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [4]).`DRPAII-AC`

RPAx calculation (see Refs. [1] and [5]).`RPAXI-AC`

RPAx calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [5]).`RPAXII-AC`

dRPA calculation using an alternative derivation (see Ref. [7]).`DRPAI-AC-ALT`

RPAx calculation using an alternative derivation (see Ref. [7]).`RPAXI-AC-ALT`

dRPA calculation without integration along the adiabatic connection (using the “kinetic” and “potential” contributions, sometimes called the “alternative plasmon formula”, see Ref. [1]).`DRPAI-AC-NOINT`

RPAx calculation without integration along the adiabatic connection (using the “kinetic” and “potential” contributions, sometimes called the “alternative plasmon formula”, see Ref. [1]).`RPAXII-AC-NOINT`

For the DIEL formulation, the available methods are:

dRPA calculation (see Ref. [2]).`DRPAI-DIEL`

dRPA-IIa approximation (see Ref. [2]).`DRPAIIA-DIEL`

RPAx-Ia approximation (see Ref. [2]).`RPAXIA-DIEL`

For the RCCD formulation, the available methods are:

dRPA-I calculation (see Ref. [3]).`DRPAI-RCCD`

RPAx-II calculation (Szabo-Ostlund variant 1 is calculated too, see Ref. [3]).`RPAXII-RCCD`

Szabo-Ostlund variant 2 (see Ref. [3]).`SO2-RCCD`

dRPA-I+SOSEX correction (see Ref. [3]).`SOSEX-RCCD`

RPAX2 approximation (see Ref. [8]).`RPAX2-RCCD`

For the PLASMON formulation, the available methods are:

dRPA-I calculation (see Ref. [1]).`DRPAI-PLASMON`

RPAx-II calculation (see Ref. [1]).`RPAXII-PLASMON`

Note that to all these keywords are associated energy variables defined as :

`ECORR_VARIANT_FORMULATION_ALTERNATIVE`

(see the examples below).

Example of a dRPA-I calculation using the PBE functional:

{rks,pbe,orbital,2101.2} {rhf,start=2101.2,maxit=0} {rpatddft; orb,2101.2; ecorr,DRPAI-AC } e=ECORR_DRPAI_AC

Example of a range-separated RPAx-I calculation using the short-range PBE exchange-correlation functional and the range-separated parameter mu=0.5:

{int;erf,0.5} {rks,exerfpbe,ecerfpbe;rangehybrid;orbital,2101.2} {rpatddft; orb,2101.2; ecorr,RPAXI-AC }

Example of several RPA calculations in the same run:

{rhf,orbital,2101.2} {rpatddft; orb,2101.2; ecorr,DRPAI-AC,RPAXII-RCCD,DRPAI-DIEL } e1=ECORR_DRPAI_AC e2=ECORR_RPAXII_RCCD e3=ECORR_DRPAI_DIEL

(this way, the calculations are done with the same transformed integrals, *i.e.* without redoing the integral transformation).

Example of a dRPA-I gradient calculation:

{rks,pbe,orbital,2101.2} {rpatddft; orb,2101.2; ecorr,DRPAI-RCCD } force

Example of a geometry optimization at the LDA+dRPA-I level:

{int;erf,0.5} {rks,exerf,ecerf;rangehybrid;orbital,2101.2} {rpatddft; orb,2101.2; ecorr,DRPAI-RCCD } optg

Calculation of properties, excitation energies and oscillator strengths

`EXCIT, METHOD=<method>`

The `EXCIT`

calculations output shows the excitation energies in ua, eV and nm, the oscillator strengths in length and velocity gauge, as well as the major excitations involved in each mode. The methods available are:

Direct random-phase approximation (or time-dependent Hartree).`DRPA`

Time-dependent Hartree-Fock.`TDHF`

Time-dependent density-functional theory.`TDDFT`

Range-separated time-dependent density-functional theory [11].`RS-TDDFT`

The exchange density functionals (FUNCX) available are:

(short-range LDA exchange density functional for the erf interaction [12]).`LDAXERF`

The correlation density functionals (FUNCC) available are:

(Perdew-Wang-92 LDA correlation density functional)`LDAC`

(short-range LDA correlation density functional for the erf interaction [13]).`LDACERF`

Example of a range-separated time-dependent density-functional theory calculation using the short-range LDA exchange-correlation functional and the range-separated parameter mu=0.5:

mu=0.5 {int;erf,mu;save} {rks,exerf,ecerf;rangehybrid;orbital,2101.2} {int} {setmu,mu} {rpatddft; orb,2101.2; excit,method=rs-tddft; dftkernel,funcx=ldaxerf,funcc=ldacerf }

Example of a TDHF-TDA calculation with writing of several files for interfacing with a real-time propagation code (see Ref. [14]):

{rpatddft; integral,1; writefile,transmom.dat,energies.dat,virtual.dat,transmom_v.dat,amplitudes.dat,transmom_a.dat; orb,2330.2; excit,method=tdhf; NOSPINBLOCK; tda}

The files are: `transmom.dat`

: transition moments in position form; `energies.dat`

: total energies of the states; `virtual.dat`

: virtual positive-energy orbitals; `transmom_v.dat`

: transition moments in velocity form; `amplitudes.dat`

: coefficients of excited-state wave functions over single-excited determinants; `transmom_a.dat`

: transition moments in acceleration form

References

$[1]$ J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. **7**, 3116 (2011).

$[2]$ G. Jansen, B. Mussard, D. Rocca, J. G. Ángyán (in prep).

$[3]$ J. Toulouse, W. Zhu, A. Savin, G. Jansen, and J. G. Ángyán, J. Chem. Phys. **135**, 084119 (2011).

$[4]$ F. Furche, Phys. Rev. B **64**, 195120 (2001).

$[5]$ J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. **102**, 096404 (2009).

$[6]$ B. Mussard, P. Reinhardt, J. G. Ángyán, and J. Toulouse, J. Chem. Phys. (submitted).

$[7]$ Heßelmann, A., Görling, A., Phys. Rev. Lett. **106**, 093001 (2011).

$[8]$ Heßelmann, A., Phys. Rev. A **85**, 012517 (2012).

$[9]$ B. Mussard, P. G. Szalay, J. G. Ángyán, J. Chem. Theory Comput. **10**, 1968 (2014).

$[10]$ J. Toulouse, E. Rebolini, T. Gould, J. F. Dobson, P. Seal, J. G. Ángyán, J. Chem. Phys. **138**, 194106 (2013).

$[11]$ E. Rebolini, A. Savin, J. Toulouse, Mol. Phys. **111**, 1219 (2013).

$[12]$ J. Toulouse, A. Savin, and H.-J. Flad, Int. J. Quantum Chem. **100**, 1047 (2004).

$[13]$ S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys. Rev. B **73**, 155111 (2006).

$[14]$ E. Coccia, B. Mussard, M. Lebeye, J. Caillat, R. Taieb, J. Toulouse, and E. Luppi, Int. J. Quant. Chem. **116**, 1120 (2016).

## Fractional orbital occupation calculations

Spin-unrestricted UHF or UKS calculations, as well as RPA calculations, can be done with fractional orbital occupation numbers. In particular, this allows one to perform calculations with fractional electron numbers, useful for example for analysing approximations (see Ref. [1]).

Occupation numbers can be given separately for alpha and beta orbitals, respectively, by using the following commands *before* calling the `uhf`

or `uks`

program:

`FRACOCCA,<orbital_index>,<occupation_number>`

`FRACOCCB,<orbital_index>,<occupation_number>`

where `orbital_index`

is the index of the orbital in the order used for constructing the density matrix, and `occupation_number`

is the occupation number of this orbital between 0 and 1. These commands can be repeated in case of several fractionally occupied orbitals. When specifying the charge and the wave function, all fractionally occupied orbitals must be considered as occupied orbitals.

For DFT calculations, the implementation has only been done for the option `matrix=0`

in the `uks`

command.

Example of a PBE calculation on the fractional C cation with 5.3 electrons:

symmetry,nosym geom={C} fracocca,4,0.3 {uks,pbe,matrix=0;wf,6,0,2}

Example of a RSH calculation on the fractional CO anion with 14.8 electrons:

symmetry,nosym angstrom geom={ 2 C 0 0 0 O 0 0 1.128} set,charge=-1 fracocca,8,0.8 {int;erf,0.5} {uks,ecerfpbe,exerfpbe,matrix=0;rangehybrid;wf,15,0,1}

Example of a HF calculation on the H atom with 0.5 alpha electron and 0.5 beta electron:

symmetry,nosym geom={H} set,charge=-1 fracocca,1,.5 fracoccb,1,.5 {uhf;wf,2,0,0}

Subsequent RPA correlation calculations (see here) will be automatically done with fractional orbital occupation numbers.

References

$[1]$ B. Mussard, J. Toulouse, Mol. Phys., **115**, 161 (2017).

## Examples

The following shows the use of both non-self-consistent and self-consistent DFT.

- examples/cndft.inp
geometry={c;n,c,r} r=1.1 angstrom dfunc,b-lyp rhf;method(1)=program dft;edf(1)=dftfun uhf;method(2)=program dft;edf(2)=dftfun uks;method(3)=program,edf(3)=dftfun dft;method(4)=program,edf(4)=dftfun table,dftname,dftfuns table,method,edf