The density functional program

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

  • DFT calculate functional of a previously computed density.
  • RKS or RKS-SCF calls the spin-restricted Kohn-Sham program. KS and KS-SCF are aliases for RKS.
  • UKS or UKS-SCF calls the spin-unrestricted Kohn-Sham program

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).

The following options may be specified on the KS or UKS command lines:

  • GRID=target Specifies the grid target accuracy (per atom). The default is 1.d-6 unless this has been modified using a global THRESH, GRID option.
  • COARSE (logical) If true, perform initial iterations with a coarser grid. Default is false.
  • GRIDMAX=gridmax In the initial iterations, the grid accuracy is min(gridmax, target*coarsefac) (only if COARSE is set).
  • COARSEFAC=coarsefac Factor for initial grid accuracy (see above). The default is 1000.
  • GRIDGRAD=0 or 1 Defines whether grid weight derivatives are included in analytical gradient calculations (default: 0). Disabling these can improve convergence in geometry optimizations.
  • DFTFAC=[fac1,fac2,..] Factors for each functional. The number of given values must agree with the number of functionals.
  • EXFAC=factor Fraction of exact exchange added to the functional. The default depends on the functional.
  • TOLORB=value Threshold for orbital screening (default depends on energy threshold).

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

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,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.


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

  • TOTAL Overall target accuracy (per atom) of density functional. Defaults to the value of the global thresholdcommand:gthresh 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 Orbital truncation threshold.
  • DENSITY Density truncation threshold.
  • FOCK Fock matrix truncation threshold.


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.

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]):

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

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

Unrestricted calculations (uks, ump2) can also be done.

$[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).

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:


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


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

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.

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


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.


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.


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).


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)).

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 top-level command, which should be presented before the the DFT or KS commands that will use the grid. Alternatively, GRID and its subcommands can be presented as directives within the KS program.


The integration grid is stored on record orb.file (default 1800.2). The information on disk consists of two parts: the parameters necessary to define the grid, and a cache of the evaluated grid points and weights. 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 subcommands. The currently implemented default parameters are equivalent to the following input commands.


For alternative grid settings, see section alternative numerical integration grid control (GRID).


Specify the target accuracy of integration. Radial and angular grids are generated adaptively, with the aim of integrating the Slater-Dirac 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. However, they can be adjusted individually by specifying accr and acca respectively.


Specify the details of the radial quadrature scheme. Four different radial schemes are available, specified by method = EM, BECKE, AHLRICHS or LOG, with the latter being the default.

EM 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 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.

$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 as given by the THR command 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.


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. The 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).

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$.

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 the Slater-Dirac exchange functional using a sum of approximate atomic densities. If acca is zero, the global threshold is used instead, or else it is ignored.

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. This test functional is computed from the superposition of 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 distance of the closest atom. FACNEIGH is currently set to 2.5 by default and governs the shape of the test density on the sphere. Increasing this value will implicitly lead to more angular points at (particular) 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.


Controls Becke-Voronoi partitioning of space. 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.


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


forces the use of a grid cache where possible.


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


forces the use of any point-group symmetry.


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.

There are three alternative types of pruned grids that can be used; PRUNED, NEESE and SG. They can be selected by setting the name option:








NEESE and SG are fixed-pruned grids and are defined in sections neese Grids (GRID,NAME=NEESE) 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 PRUNED Grid Thresholds, PRUNED Grid Radial Scheme and PRUNED Grid Voronoi Partitioning.

The on-the-fly pruned grid size should normally be controlled using threshold_overall=value. This is sets the value of threshold_radial and threshold_angular, which are described below. The default is:


or the two thresholds can be set independently:


The size of the radial grid is evaluated from the threshold_radial option: $$n_{\text{rad}}(A) = -\log_{10} t - 30 + (n_A - 1) * 15$$ where $t$ is the threshold and $n_A$ is the period of atom. The default value of this threshold is 1e-6. The minimum and maximum values of the number of radial points is also restricted by:




for each atom based on its periodic row.

On-the-fly grid pruning generates sets the angular grid order for each radial points based on the molecular geometry. A target integrand is evaluated using the maximum allowed Lebedev angular order, with the shells angular order increased until it agrees to less than the threshold_angular option. This is set:


The default threshold is 1e-6.

The maximum allowed angular order is set by the row in the periodic table:


The default value is [53,53,53,53]. Additionally, the minimum angular order can be set:


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.

There are two target integrands; Slater atomic density (SAD) and a pseudo-density functional (PW91). The Slater atomic density is default and can be explicitly specified:


it is a purely spherical density at each atom and each shells contribution is geometry dependent. The angular nature of the shell is generated from the surrounding atoms and the Voronoi partitioning.

Alternatively a pseudo-density functional can be used:


whereby the density is the Slater atomic density.

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:


but can result in a less accurate grid.

The radial quadrature used in the on-the-fly pruning can be set to the default:










Each scheme typically has a scaling value. These can be adapted using:


where a scaling for every atom must be specified or they will be ignored.

The molecular grid is a sum of atomic grids. To avoid overcounting in the overlapping regions, Voronoi partitioning is typically used. Fuzzy cell boundaries are employed using either Murray (default), Becke or Stratmann’s step functions.






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. This can be modified from the default values for each case:




The NEESE grid is a fixed-pruned grid.

The size of the grid should set using the neese_index:


which links the angular size using neese_max_order and radial size using neese_int_acc, which are defined below. The value can be set from 3 to 12. Using neese_index=3 sets neese_int_acc=4.34 and neese_max_order=23. Increasing the index by 1 increases neese_int_acc by 0.33 (effectively 5 radial points) and neese_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:


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:


where value would take values from 4.34 upwards.

$[1]$ Krack, M., Köster, A.M., 1998. An adaptive numerical integrator for molecular integrals. J. Chem. Phys. 108, 3226-3234.

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




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.

$[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.
$[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.

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 B88dftfun:B88 1 10.1103/PhysRevA.38.3098
B-LYP B88dftfun:B88:LYPdftfun:LYP 1:1
BLYP B88dftfun:B88:LYPdftfun:LYP 1:1
B-P B88dftfun:B88:P86dftfun:P86 1:1
BP86 B88dftfun:B88:P86dftfun:P86 1:1
B-VWN B88dftfun:B88:VWN5dftfun:VWN5 1:1
B3LYP EXACTdftfun:EXACT:B88dftfun:B88:DIRACdftfun:DIRAC:LYPdftfun:LYP:VWN5dftfun:VWN5 0.2:0.72:0.08:0.81:0.19
B3LYP3 EXACTdftfun:EXACT:B88dftfun:B88:DIRACdftfun:DIRAC:LYPdftfun:LYP:VWN3dftfun:VWN3 0.2:0.72:0.08:0.81:0.19
B3LYP5 EXACTdftfun:EXACT:B88dftfun:B88:DIRACdftfun:DIRAC:LYPdftfun:LYP:VWN5dftfun:VWN5 0.2:0.72:0.08:0.81:0.19
B88X B88dftfun:B88 1 10.1103/PhysRevA.38.3098
B97 EXACTdftfun:EXACT:B97DFdftfun:B97DF 0.1943:1
B97R EXACTdftfun:EXACT:B97RDFdftfun:B97RDF 0.21:1
BECKE B88dftfun:B88 1 10.1103/PhysRevA.38.3098
BH-LYP EXACTdftfun:EXACT:B88dftfun:B88:LYPdftfun:LYP 0.5:0.5:1
CS CS1dftfun:CS1 1
HFB B88dftfun:B88 1 10.1103/PhysRevA.38.3098
LDA DIRACdftfun:DIRAC:VWN5dftfun:VWN5 1:1
LSDAC PW92Cdftfun:PW92C 1 10.1103/PhysRevB.45.13244
LSDC PW92Cdftfun:PW92C 1 10.1103/PhysRevB.45.13244
LYP88 LYPdftfun:LYP 1
MM05 EXACTdftfun:EXACT:M05Xdftfun:M05X:M05Cdftfun:M05C 0.28:0.72:1 10.1063/1.2126975
MM05-2X EXACTdftfun:EXACT:M052XXdftfun:M052XX:M052XCdftfun:M052XC 0.56:0.44:1 10.1021/ct0502763
MM06 EXACTdftfun:EXACT:M06Xdftfun:M06X:M06Cdftfun:M06C 0.27:0.73:1 10.1007/s00214-007-0310-x
MM06-2X EXACTdftfun:EXACT:M062XXdftfun:M062XX:M062XCdftfun:M062XC 0.54:0.46:1 10.1007/s00214-007-0310-x
MM06-L M06LXdftfun:M06LX:M06LCdftfun:M06LC 1:1 10.1063/1.2370993
MM06-HF EXACTdftfun:EXACT:M06HFXdftfun:M06HFX:M06HFCdftfun:M06HFC 1:1:1 10.1021/jp066479k
PBE PBEXdftfun:PBEX:PBECdftfun:PBEC 1:1 10.1103/PhysRevLett.77.3865
PBE0 EXACTdftfun:EXACT:PBEXdftfun:PBEX:PBECdftfun:PBEC 0.25:0.75:1 10.1063/1.478522
PBE0MOL EXACTdftfun:EXACT:PBEXdftfun:PBEX:PW91Cdftfun:PW91C 0.25:0.75:1
PW91 PW91Xdftfun:PW91X:PW91Cdftfun:PW91C 1:1 10.1103/PhysRevB.46.6671
S-VWN DIRACdftfun:DIRAC:VWN5dftfun:VWN5 1:1
VS99 VSXCdftfun:VSXC 1
VWN VWN5dftfun:VWN5 1
VWN80 VWN5dftfun:VWN5 1
M05 EXACTdftfun:EXACT:XC-M05dftfun:XC-M05 0.28:1 10.1063/1.2126975
M05-2X EXACTdftfun:EXACT:XC-M05-2Xdftfun:XC-M05-2X 0.56:1 10.1021/ct0502763
M06 EXACTdftfun:EXACT:XC-M06dftfun:XC-M06 0.27:1 10.1007/s00214-007-0310-x
M06-2X EXACTdftfun:EXACT:XC-M06-2Xdftfun:XC-M06-2X 0.54:1 10.1007/s00214-007-0310-x
M06-L XC-M06-Ldftfun:XC-M06-L 1 10.1063/1.2370993
M06-HF EXACTdftfun:EXACT:XC-M06-HFdftfun:XC-M06-HF 1:1 10.1021/jp066479k
M08-HX EXACTdftfun:EXACT:XC-M08-HXdftfun:XC-M08-HX 0.5223:1
M08-SO EXACTdftfun:EXACT:XC-M08-SOdftfun:XC-M08-SO 0.5679:1
M11-L XC-M11-Ldftfun:XC-M11-L 1
TPSS TPSSXdftfun:TPSSX:TPSSCdftfun:TPSSC 1:1 10.1103/PhysRevLett.91.146401
TPSSH EXACTdftfun:EXACT:TPSSXdftfun:TPSSX:TPSSCdftfun:TPSSC 0.1:0.9:1 10.1063/1.1626543
M12HFC EXACTdftfun:EXACT:M12Cdftfun:M12C 1:1 10.1063/1.4768228
HJSWPBE HJSWPBEXdftfun:HJSWPBEX:PBECdftfun:PBEC 1:1 10.1063/1.2921797
HJSWPBEH EXACTdftfun:EXACT:HJSWPBEXdftfun:HJSWPBEX:PBECdftfun:PBEC 0.2:0.8:1 10.1063/1.3073302
PBESOL PBESOLXdftfun:PBESOLX:PBESOLCdftfun:PBESOLC 1:1 10.1103/PhysRevLett.100.136406

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$.

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:


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).

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


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


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


  • B88 Becke 1988 Exchange Functional


  • B88C 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.


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


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


  • B97RDF 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.


  • BR Becke-Roussel Exchange Functional. A. D. Becke and M. R. Roussel,Phys. Rev. A 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.$$
  • BRUEG Becke-Roussel Exchange Functional — Uniform Electron Gas Limit. A. D. Becke and M. R. Roussel,Phys. Rev. A 39, 3761 (1989) As for BR but with ${\gamma=0.8}$.
  • BW Becke-Wigner Exchange-Correlation Functional. Hybrid exchange-correlation functional comprimising Becke’s 1998 exchange and Wigner’s spin-polarised correlation functionals.


  • CS1 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 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.
  • CS2 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 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$.
  • DIRAC Slater-Dirac Exchange Energy. Automatically generated Slater-Dirac exchange.


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

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}$.

  • ECERFPBE Range-Separated Correlation Functional. Toulouse-Colonna-Savin range-separated correlation functional based on PBE, see J. Toulouse et al., J. Chem. Phys. 122, 014110 (2005).


  • EXACT Exact Exchange Functional. Hartree-Fock exact exchange functional can be used to construct hybrid exchange-correlation functional.
  • EXERF Short-range LDA correlation functional. Local-density approximation of exchange energy

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}$.

  • EXERFPBE Range-Separated Exchange Functional. Toulouse-Colonna-Savin range-separated exchange functional based on PBE, see J. Toulouse et al., J. Chem. Phys. 122, 014110 (2005).


  • G96 Gill’s 1996 Gradient Corrected Exchange Functional


  • HCTH120 Handy least squares fitted functional


  • HCTH147 Handy least squares fitted functional


  • HCTH93 Handy least squares fitted functional


  • HJSWPBEX 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. 128, 194105 (2008).


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


  • LYP Lee, Yang and Parr Correlation Functional. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37, 785(1988); B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett. 157, 200 (1989).


  • M052XC M05-2X Meta-GGA Correlation Functional


  • M052XX M05-2X Meta-GGA Exchange Functional


  • M05C M05 Meta-GGA Correlation Functional


  • M05X M05 Meta-GGA Exchange Functional


  • M062XC M06-2X Meta-GGA Correlation Functional


  • M062XX M06-2X Meta-GGA Exchange Functional


  • M06C M06 Meta-GGA Correlation Functional


  • M06HFC M06-HF Meta-GGA Correlation Functional


  • M06HFX M06-HF Meta-GGA Exchange Functional


  • M06LC M06-L Meta-GGA Correlation Functional


  • M06LX M06-L Meta-GGA Exchange Functional


  • M06X M06 Meta-GGA Exchange Functional


  • M12C Meta GGA Correlation Functional. Meta-GGA correlation functional based on first principles, see M. Modrzejewski et al., J. Chem. Phys. 137, 204121 (2012).


  • MK00 Exchange Functional for Accurate Virtual Orbital Energies


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


  • P86 .. Gradient correction to VWN.


  • PBEC PBE Correlation Functional


  • PBESOLC PBEsol Correlation Functional


  • PBESOLX PBEsol Exchange Functional


  • PBEX PBE Exchange Functional


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


  • PW86 .. GGA Exchange Functional.


  • PW91C Perdew-Wang 1991 GGA Correlation Functional


  • PW91X Perdew-Wang 1991 GGA Exchange Functional


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


  • STEST Test for number of electrons
  • TFKE Thomas-Fermi Kinetic Energy. Automatically generated Thomas-Fermi Kinetic Energy.


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


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


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


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


  • THGFC .. 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$.


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


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


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


  • TPSSC TPSS Correlation Functional. J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).


  • TPSSX TPSS Exchange Functional. J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).


  • VSXC .


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


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


  • VWN5 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.


  • XC-M05-2X 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. 2, 364 (2006).


  • XC-M05 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. 123, 161103 (2005).


  • XC-M06-2X 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-HF 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-L M06-L Meta-GGA Exchange-Correlation Functional. Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).


  • XC-M06 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-M08-HX 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. 4, 1849 (2008).


  • XC-M08-SO 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. 4, 1849 (2008).


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


  • XC-SOGGA SOGGA Exchange-Correlation Functional. Y. Zhao and D. G. Truhlar, J. Chem. Phys. 128, 184109 (2008).


  • XC-SOGGA11-X 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. 135, 191102 (2011).


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


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


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


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}

will have to be entered in the input file for performing calculations with the wB97X functional by Chai and Head-Gordon

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:


For getting the corresponding wB97X-D functional, add the D2-type dispersion correction

df-ks,XC_HYB_GGA_XC_WB97X; disp,1

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 Molpro currently supports calculations with the VV10 (Vydrov & Van Voorhis 2010) functional The local part for this functional can be taken from LibXC:


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

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


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"?>
 <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)/

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.

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.

$[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)

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


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}

The following options can be used with the LHF program:

  • ORB orbital record containing previous orbitals (mandatory)
  • LHFFIT use density fitting to compute the exchange matrix
  • POT record in which to write the local exchange potential
  • RESTART recompute the LHF orbitals and orbital energies from an already computed exchange potential
  • SLA if set to 2, compute the Slater potential with the method described in appendix A in Ref. [1] (default: 1)
  • BLOCK batch size for numerical grid
  • PUNCH 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.
  • DAMP 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.

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:


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.:


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}

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

  • POISSON 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.
  • INT switch 3-index integral routine to be used (can take the values 1,2 or 3 (default is 2))

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}

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}

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


(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


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


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:

  • MAXIT maximum number of iterations
  • SCFTHR threshold for energy convergence
  • OEP choose between 7 (!) different OEP routines (an explanation of all these would blow up this section significantly, default is OEP=1)
  • PRINT debug print level
  • FERAMA use Fermi-Amaldi potential (instead of OEP potential)
  • ELP use Effective Local Potential method by Staroverov et al. (J. Chem. Phys. 125, 081104 (2006))
  • DFIT if set to $\ne$ 0, use density fitting for two-electron integrals
  • DIRECT set to $\ne$ 0 for integral-direct calculation (not very efficient at the moment)
  • ITVREF 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).

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

  • VREF can be set to ’BASIS’ or ’BASIS2’ according to the two schemes described in Ref. [5] to determine the coefficients for the reference potential
  • VREST can be set to zero to omit calculation of the rest potential
  • AUX norm (should be ’J’)
  • ATOM1 set (x,y,z) coordinates for a starting position used for plotting the potential
  • ATOM2 set (x,y,z) coordinates for a terminal position used for plotting the potential
  • SYMCEN1/2 do not use!
  • TEST can be used to switch on several tests during the generation of the OEP potential
  • POT 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)
  • CHARGE set to $\ne$ 0 for activating the charge condition
  • HOMO set to $\ne$ 0 for activating the HOMO condition
  • HOMO1 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)
  • RESPEIG set to $\ne$ 0 to print out the eigenvalues of the response matrix (useful to test stability of the OEP method)
  • RHS set to $\ne$ 0 for printing the right-hand side in the output
  • XPRINT set to $\ne$ 0 for printing the response matrix
  • SOLVE matrix inversion methods to solve the OEP equation. The different options are: GESV (default), SVD, SVD2, TIKH
  • THR threshold value for singular value decomposition (SVD) or Tikhonov regularisation (TIKH)
  • FILTER an integer number can be entered to filter out the response function eigenvectors with the lowest eigenvalues in an SVD matrix inversion
  • SVD here, SVD is used to perform a singular value decomposition of the (aux$|$occ$\times$virt) 3-index integrals
  • SVDTHR the corresponding SVD threshold for this
  • CONTR set to $\ne$ 0 to use contraction of auxiliary basis through the SVD of (aux$|$occ$\times$virt) integrals
  • LAPLACE solve OEP within the AO basis set making use of the Laplace transformation
  • VXDUMP dump the nonlocal and local exchange potentials to a record
  • OEPDUMP dump record for OEP which contains the orbitals, the density and the eigenvalues
  • SAVE record in which the coefficients of the potential can be written to
  • READ record from which the coefficients can be read
  • OEPPRINT debug print level
  • HOMOORB 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)
  • FIXEIG can be used to specify the number of an orbital whose eigenvalue should be fixed (see also HOMO1 option)
  • EQNORM set the norm for the OEP equation (not available for all OEP routines)
  • ELPMODE switch the mode in ELP calculations
  • EXIT can be used to deallocate some arrays that are used during the OEP calculation
  • CONV 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.
  • VXDIFF 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.
  • THRS an eigenvalue threshold used in the decomposition of the overlap matrix, used to reduce the auxiliary function space (not available in all OEP routines)
  • GAMMA Tikhonov regularisation factor (only used for OEP=7 case). A value of $\gamma=3\cdot 10^{-4}$ should conventionally be used.

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}

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:

  • XFAX factor for exact exchange
  • ADD set to $\ne$ 0 in order to add the GGA contribution to the xc matrix to the rhs of the OEP equation
  • SHIFT if set to a value $\ne$ 0 the asymptotic correction method for the xc potential described in Ref. [6] will be used

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}
    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.

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 would look like

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

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:

  • DEN Dump record containing the density matrix (must be supplied)
  • TOTEN If TOTEN>0 the total energy will be printed (and stored in the intrinsic energy variable) reading in the previous SCF energy
  • GTHR This 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=<global grid threshold set in dft/ks>.
  • NBLOCK Block size for grid batches. Default: NBLOCK=0 which actually means that no grid batching is performed at all (requiring more memory consumption).
  • B The 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 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:


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 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>

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

Currently the default dispersion correction added to the DFT energy is the D3 dispersion correction developed by Grimme et al., see Ref. [3]. The disp,2 keyword can have the following additional options:

  • FUNC Functional name (default: FUNC='pbe').
  • VERSION Can have values 2 and 3 according to parametrisations from Refs. [2] and [3] (default: VERSION=3)
  • ANAL Performs a detailed analysis of pair contributions.
  • GRAD Cartesian gradients are computed.
  • ABC Calculate three-body dispersion contribution if ABC$\ne$0.
  • NOABC Do not calculate three-body dispersion contribution if NOABC$\ne$0.
  • NUM Calculate gradient numerically.
  • OLD Same as VERSION=2.
  • ZERO DFT-D3 original zero-damping.
  • BJ DFT-D3 with Becke-Johnson finite-damping.
  • BJM Revised DFT-D3 with Becke-Johnson damping.
  • CNTHR Neglect threshold in Bohr for CN (default=40).
  • CUTOFF Neglect threshold in Bohr for $E_{\text{disp}}$, (default=95).
  • TZ 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.

(see also ' for further documentation).

The DFT-D4 program (disp,3) accepts the options:

  • CHRG Molecular charge.
  • S6 Scaling factor for leading order term.
  • S8 Scaling factor for $C_8$ dependent term.
  • S9 Scaling factor for $C_9$ dependent term.
  • A1 Parameter in Becke-Johnson damping function.
  • A2 Parameter in Becke-Johnson damping function.
  • DFA Name of density functional (useful if the name conventions for functionals in the D4 program are different from Molpro's).
  • LMBD Treatment of many-body dispersion. 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 term
  • NOMB same as LMBD=0
  • GRAD Calculate gradient if LMBD=0$\ne 0$.
  • HESS Calculate hessian (currently not used in Molpro).
  • PRINT Set print level (always highest is set automatically with dispcorr4 command).
  • G_A Charge scale height.
  • G_C Charge scale steepness.
  • REFQ Kind of charge model.

(see also ' for further documentation).

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


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


with the following options to dispcorr:

  • MODE Adjusts the parametrisation used: MODE=1 uses parameters from Ref. [1] and MODE=2 uses parameters from Ref. [2] (default: MODE=1)
  • SCALE,S6 Overall scaling parameter $s_6$ (see Eq. \eqref{eq:edisp}) and Refs. [2,3] for optimised values).
  • ALPHA Damping function parameter (see Eq. \eqref{eq:fdamp}). Smaller values lead to larger corrections for intermediate distances.
  • DAMP Choose damping function 'GRIMME' or 'CHAI' (default:'GRIMME'). With 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]
  • A Parameter a in the damping function of Chai et al., see Ref. [5]
  • GRAD Calculate gradient if 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.

$[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)

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].


$[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)

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)

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

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).

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

and global options, shared by all commands inside the rpatddft block:

  • STAB Check matrices stability conditions in RPA calculations. When used without an 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).

  • C6 Computes C6 coefficients from last two saved polarizabilities.
  • TDA Tamm-Dancoff approximation for EXCIT and PROPERTIES.
  • NOMP2 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 keyword.
  • NOSPINBLOCK 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)
  • NOSPINFLIP Exclude spin-flip dimensions of unrestricted RPA calculations that use the NOSPINBLOCK formalism (not suitable for all RPA variants).
  • WRITEFILE Write files with eigenvalues, virtual orbital energies, dipole moments, dipole velocities, dipole accelerations and amplitudes from a TDA calculation
  • VIAXPY #to comment#
  • INTEGRAL,<nbr> 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 transform routine.
  • OCC,<nocc> Explicitly specify the number of occupied orbitals (useful for fake pseudopotential calculations).
  • CORE,<core> Specify core orbitals (default: last specified core orbitals or, if none, atomic inner shells)
  • PRINT,<nbr> Level of print expected from the output (from 0(default) to 3).
  • PRINT_INT,<nbr> Level of print of integrals (AO,MO,Orbtials,…), from 0 to 4.
  • PRINT_TIME,<nbr> If greater or equal to 1, will print out information on time spent in routines.

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:


For the AC formulation, the available methods are:

  • DRPAI-AC dRPA calculation (see Refs. [1] and [4]).
  • DRPAII-AC dRPA calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [4]).
  • RPAXI-AC RPAx calculation (see Refs. [1] and [5]).
  • RPAXII-AC RPAx calculation, using antisymmetrized two-electron integrals (see Refs. [1] and [5]).
  • DRPAI-AC-ALT dRPA calculation using an alternative derivation (see Ref. [7]).
  • RPAXI-AC-ALT RPAx calculation using an alternative derivation (see Ref. [7]).
  • DRPAI-AC-NOINT dRPA 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 RPAx calculation without integration along the adiabatic connection (using the “kinetic” and “potential” contributions, sometimes called the “alternative plasmon formula”, see Ref. [1]).

For the DIEL formulation, the available methods are:

  • DRPAI-DIEL dRPA calculation (see Ref. [2]).
  • DRPAIIA-DIEL dRPA-IIa approximation (see Ref. [2]).
  • RPAXIA-DIEL RPAx-Ia approximation (see Ref. [2]).

For the RCCD formulation, the available methods are:

  • DRPAI-RCCD dRPA-I calculation (see Ref. [3]).
  • RPAXII-RCCD RPAx-II calculation (Szabo-Ostlund variant 1 is calculated too, see Ref. [3]).
  • SO2-RCCD Szabo-Ostlund variant 2 (see Ref. [3]).
  • SOSEX-RCCD dRPA-I+SOSEX correction (see Ref. [3]).
  • RPAX2-RCCD RPAX2 approximation (see Ref. [8]).

For the PLASMON formulation, the available methods are:

  • DRPAI-PLASMON dRPA-I calculation (see Ref. [1]).
  • RPAXII-PLASMON RPAx-II calculation (see Ref. [1]).

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


(see the examples below).

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


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


Example of several RPA calculations in the same run:


(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:


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


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:

  • DRPA Direct random-phase approximation (or time-dependent Hartree).
  • TDHF Time-dependent Hartree-Fock.
  • TDDFT Time-dependent density-functional theory.
  • RS-TDDFT Range-separated time-dependent density-functional theory [11].

The exchange density functionals (FUNCX) available are:

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

The correlation density functionals (FUNCC) available are:

  • LDAC (Perdew-Wang-92 LDA correlation density functional)
  • LDACERF (short-range LDA correlation density functional for the erf interaction [13]).

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:


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


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

$[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).

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:



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:


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


C 0 0 0
O 0 0 1.128}

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


Subsequent RPA calculations will be automatically done with fractional orbital occupation numbers.

$[1]$ B. Mussard, J. Toulouse, Mol. Phys., to appear (2016).

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

r=1.1 angstrom