XCFun DFT library

XCFun is a library of DFT exchange-correlation (XC) functionals. It is based on automatic differentiation and can therefore generate arbitrary order derivatives of these functionals. This documentation refers to the development version of the XCFun library, make sure that you are up to date.

Installation and linking

Compile the library using one of the included makefiles (for example Makefile.gcc). This will generate a library file lib/libxcfun.a which can be statically linked to your application. Include files listing all available functionals will be generated during the compilation of XCFun. These files are include/xcfun_autogen.h (for the C interface), and fortran/xcfun_autogen.f90 for the Fortran module. C or C++ programs that uses XCFun should include the xcfun.h header file, while Fortran programs should use the xcfun module defined in fortran/xcfun_module.f90.


To use the library you should first create one xc_functional object for each functional and each thread you want to use. Each functional object can only be used by one thread at a time. After creating these objects using xc_new_functional() you should set them up by defining the variables to differentiate with respect to, as well as the parameters of the functional you want to use. After having done so you can use xc_eval() to evaluate the exchange correlation energy and its derivatives.

Example C program that evaluates BLYP to order 2 using alpha/beta variables:

#include "xcfun.h"

int main(int argc, char *argv[])
  int derivative_order = 2;
  int nr_points = 1;
  double density[5] = {1,2,3,5,4}; /* na nb ga.ga gb.gb ga.gb */
  double output[21]; /* We have 21 output numbers for derivatives up to order 2 */
  xc_functional fun = xc_new_functional();
  xc_set_mode(fun, XC_VARS_AB);
  xc_set_param(fun, XC_LYPC, 1.0);
  xc_set_param(fun, XC_BECKEX, 1.0);
  xc_eval(fun, derivative_order, nr_points, density, output);

Please see the complete API documentation below.

Input, output and units

The library uses atomic units for all input and output variables.

The XC energy density and derivatives can be evaluated using local spin-up (α) and spin-down (β) quantities. In the most general case these are

  • nα The spin-up electron number density.
  • nβ The spin-down density.
  • σαα = ∇nα·∇nα The square magnitude of the spin-up density gradient.
  • σαβ = ∇nα·∇nβ The dot product between the spin-up and spin-down gradient vectors.
  • σββ = ∇nβ·∇nβ The square magnitude of the spin-down density gradient.
  • τα = 1/2 Σi |∇ψ|2 The spin-up Kohn-Sham kinetic energy density.
  • τβ The spin-down Kohn-Sham kinetic energy density.
Alternatively you can use total density (n = nα + nβ) and spin density (s = nα - nβ) variables. These also have corresponding gradient and kinetic energy components. See xc_set_mode() below for more information.

The output is given in graded reverse lexicographical order. For example a spin-polarized second order GGA functional will give 21 output elements, starting with the XC energy density. Symbolically we may write this as a list starting with the energy E, followed by five gradient elements Eα Eβ Eσαα Eσαβ Eσββ and 15 second derivatives Eαα Eαβ Eασαα ... Eββ Eβσαα ... Eσββσββ. See the function xc_derivative_index() for information about addressing these elements.


The library is written in C++, but can also be directly used in a C project or Fortran project. The C interface is described in include/xcfun.h, while the Fortran interface is described in the module file fortran/xcfun_module.f90. This documentation tries to describe both the C and Fortran API at the same time, even though small differences may exist between the two interfaces.

Setup and Testing

const char *xcfun_splash()

Return a multi-line string describing the library. Please print this string so that your users find the right citation for the library. The Fortran version is xcfun_splash(text), and the message is put in the string text.


Return a double precision version number of the library.


Run all internal tests and return the number of failed tests.

Creating and destroying functional objects


Create a new functional object. The C version returns an object of type xc_functional. The Fortran version returns an integer. The creation of this object may be rather slow; create an object once for each calculation, not once for each grid point.


Release the memory associated with functional (previously allocated by xc_new_functional()).

Defining a functional

xc_set_mode(functional, mode)

Set functional to operate in mode, one of
  • XC_VARS_A Full spin polarization using only alpha quantities, in the order nα, σαα, τα
  • XC_VARS_N Completely unpolarized mode using only total density quantities. n, σnn, τn
  • XC_VARS_AB Using alpha and beta variables: nα, nβ, σαα, σαβ, σββ, τα, τβ
  • XC_VARS_NS Using total density/spin polarization variables: n, s, σnn, σns, σss, τn, τs
The mode and the type of functional (LDA, GGA, metaGGA) determines the number of variables used by the functional.

xc_set_param(functional, param, value)

Set the parameter param in functional to some double precision value. Param is one of the named constants listed in include/xcfun_autogen.h (or fortran/xcfun_autogen.f90). The parameters can be the weight of different functionals, or some other parameter in the functional definition. Example of parameters are XC_VWN5C, the amount of VWN5 LDA correlation, and XC_RANGESEP_MU, the range separation constant μ used in range separated hybrid functionals. You can iterate over all parameters by using the XC_NR_PARAMS constant. C parameters start at 0 and Fortran parameters start at 1.

The functional weights default to 0, while other parameters default to some sensible value. The current value of a parameter can be checked with the double precision function

xc_get_param(functional, param)

which returns the current value of param. Note that different functionals may have different parameter values set.


Returns the name of param. This name is the same as the symbolic constant used in the programming interface, so the first three characters (XC_) can be stripped off before printing the string to the user. The Fortran version is xc_name(param, name).


Return a single line (without newline) describing param. This may be a null string if the parameter was not described in the library. The Fortran version puts the description in the second argument to this subroutine.


Return a multi-line string (if present) ending with a newline. This string describes the functional or parameter in more detail, and should include references for citation and further information. The Fortran version puts the description in the second argument to this subroutine.

Evaluating a functional

Given a functional that has been set up using the functions above we may evaluate this functional up to the highest order defined in src/config.h. Increasing the maximum order increases the compile time and size of the library, so it is important to keep this limit rather low. The current limits of the library may be queried with the function


which returns the highest derivative order supported by the previously defined functional.


Returns one of the integers
  • XC_LDA
  • XC_GGA
depending on the parameters set in the functional.


Returns the length of the input array (i.e. the number of density variables. This depends on the type of functional and on the evaluation mode. Typical values are 2 (LDA with incomplete spin polarization) or 5 (GGA with incomplete spin polarization).

xc_output_length(functional, order)

Return the number of output coefficients computed by xc_eval() if functional is evaluate up to order. Note that all derivatives up to order are calculated, not only those of the particular order.

xc_regularize_density(functional, density)

Modify density in-place so that functional can be evaluated without infinities or floating point errors. This mainly involves making sure that the density is physical and nonzero. In cases where the derivatives should be infinite they will instead be very large, but otherwise the effect on the output should be minimal.

C xc_eval(functional, order, density,result)

C Vector version xc_eval_vec(functional, order, nr_points, density, density_pitch, result, result_pitch)

Fortran xc_eval(functional, order, nr_points, densities, results)

Evaluate the previously defined functional at nr_points different points, with derivatives up to order. The input array density contains nr_points*N density values, where N is the value returned by xc_input_length(). Similarly for the output array result. The input and output arrays have the values for each point packed together. The output values are stored in graded lexicographical order. For an LDA functional with alpha/beta variables we would therefore have the order E Ea Eb Eaa Eab Ebb Eaaa Eaab Eabb Ebbb ...

The C vector version accepts two 'pitch' arguments, which is the distance in memory between the first elements of each point in the input and output arrays.

The Fortran always works on an array of points, and automatically determines the pitch values. Note that the input numbers for each point must be consecutive in memory (stride must be one).

xc_potential(functional, density, energy, potential)

Compute the XC energy and potential (potential[0] is alpha and potential[1] is beta) at a single point. For LDA this is simply vαxc = dE/dnα and for GGA the expression is vαxc = dE/dnα - ∇·dE/d∇nα. The latter case thus needs the evaluation of the second derivatives of the functional. In the case of GGA the values of the density laplacians should be given at the end of the density array, giving at total of seven numbers in the XC_VARS_AB case (which is the only one supported at the moment).

This function has to be called one point at a time.

MetaGGA's are not supported at the moment, although this is in principle possible to implement for functionals that depend on the density laplacian (but not on the kinetic energy density).

xc_derivative_index(functional, derivative)

Given an integer array of "exponents" this function returns the index into the output array of xc_eval where the corresponding derivative is located. Note that this position depends on the functional type and mode. Alternatively you can, for low orders, use the predefined constants of the form XC_D01000 (which corresponds to the first derivative with respect to the second variable, for setups with five variables). To use these predefined constants you must therefore know the type of your functional beforehand. For LDA in a two variable mode the corresponding constant would be XC_D01, but this is a different index from XC_D01000.

Compile Time Options

The file src/config.h contains a number of compile time options. These are
  • XC_MAX_ORDER Define this to the highest order of derivatives that the library should support. Too high values make the library huge in size.
  • XC_NO_REGULARIZATION Define this to turn off regularization of the input densities in xc_eval(). If you use this option and feed the library unphysical densities it can and will stop or crash. Use mainly to investigate numerical behaviour near singular points.
  • XC_NO_SCHWARZ_REGULARIZATION Do not strictly enforce the Schwarz inequality for the gradient norms. This inequality can be violated without encountering any problems for most functionals, and it's suggested that this option is kept active (i.e. the regularization is turned off).
  • XCFUN_IN_DEVELOPMENT Enable code that is still "in development". This code is probably marked "in development" for a reason, and may contain bugs.
  • XCFUN_REF_PW92C Use the truncated valued from the reference implementation for some exact constants in the PW92C correlation functional. This can help getting an exact agreement with other implementations. Functionals that use PW92C (i.e. PBEC) are also affected by this options.