Libxc 5.1.x

Compatibility notice

This manual is for Libxc 5.1.x and newer.

See the Change Log for a complete list of changes made to the library API over the different Libxc versions.

An example

Probably the best way to explain the usage of Libxc is through an example. The following small program calculates the xc energy for a given functional for several values of the density; the available C bindings can be found in header file xc.h. When compiling this example it is necessary to link with libxc.a or libxc.so, depending if you compiled the static or shared versions of the library.

#include <stdio.h>

#include <xc.h>

int main()
{
  xc_func_type func;
  double rho[5] = {0.1, 0.2, 0.3, 0.4, 0.5};
  double sigma[5] = {0.2, 0.3, 0.4, 0.5, 0.6};
  double exc[5];
  int i, vmajor, vminor, vmicro, func_id = 1;
  /* Get the libxc version */
  xc_version(&vmajor, &vminor, &vmicro);
  printf("Libxc version: %d.%d.%d\n", vmajor, vminor, vmicro);
  /* Initialize the functional */
  if(xc_func_init(&func, func_id, XC_UNPOLARIZED) != 0){
    fprintf(stderr, "Functional '%d' not found\n", func_id);
    return 1;
  }
  /* Evaluate the energy density, depending on the family */
  switch(func.info->family)
    {
    case XC_FAMILY_LDA:
      xc_lda_exc(&func, 5, rho, exc);
      break;
    case XC_FAMILY_GGA:
    case XC_FAMILY_HYB_GGA:
      xc_gga_exc(&func, 5, rho, sigma, exc);
      break;
    }
  /* Print out density and energy density per particle */
  for(i=0; i<5; i+=1){
    printf("%lf %lf\n", rho[i], exc[i]);
  }
  /* Deallocate memory */
  xc_func_end(&func);
}

The functionals are divided in families (LDA, GGA, etc.). Given a functional identifier xc.func_id, the functional is initialized by xc_func_init, and evaluated by xc_XXX_exc, which returns the energy per particle exc. Finally the function xc_func_end cleans up.

The output of the above example should look something like this:

Libxc version: 5.1.0
0.100000 -0.342809
0.200000 -0.431912
0.300000 -0.494416
0.400000 -0.544175
0.500000 -0.586194

Fortran 90 bindings (following the Fortran 2003 standard) are also included in Libxc. These can be found in the file [libxc_master.F90] (https://gitlab.com/libxc/libxc/blob/master/src/libxc_master.F90). In general, calling Libxc from Fortran is as simple as from C. Note that you also need to link with libxcf03.a or libxcf03.so when using the Fortran 2003 interface. Here is the previous example in Fortran:

program lxctest
  use xc_f03_lib_m
  implicit none
  TYPE(xc_f03_func_t) :: xc_func
  TYPE(xc_f03_func_info_t) :: xc_info
  integer(8), parameter :: npoints = 5
  real(8) :: rho(npoints) = (/0.1, 0.2, 0.3, 0.4, 0.5/)
  real(8) :: sigma(npoints) = (/0.2, 0.3, 0.4, 0.5, 0.6/)
  real(8) :: exc(npoints)
  integer :: i, vmajor, vminor, vmicro, func_id = 1
  ! Print out the version
  call xc_f03_version(vmajor, vminor, vmicro)
  write(*,'("Libxc version: ",I1,".",I1,".",I1)') vmajor, vminor, vmicro
  ! Initialize the functional
  call xc_f03_func_init(xc_func, func_id, XC_UNPOLARIZED)
  ! Get the information
  xc_info = xc_f03_func_get_info(xc_func)
  ! Evaluate energy depending on the family
  select case (xc_f03_func_info_get_family(xc_info))
    case(XC_FAMILY_LDA)
      call xc_f03_lda_exc(xc_func, npoints, rho(1), exc(1))
    case(XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
      call xc_f03_gga_exc(xc_func, npoints, rho(1), sigma(1), exc(1))
  end select
  ! Print out density and energy density per particle
  do i = 1, npoints
    write(*,"(F8.6,1X,F9.6)") rho(i), exc(i)
  end do
  ! Deallocate memory
  call xc_f03_func_end(xc_func)
 end program lxctest

The info structure

For each functional there is a structure that holds several pieces of information. The relevant part of this structure for the end user is defined as:

#define XC_FLAGS_HAVE_EXC         (1 <<  0) /*     1 */
#define XC_FLAGS_HAVE_VXC         (1 <<  1) /*     2 */
#define XC_FLAGS_HAVE_FXC         (1 <<  2) /*     4 */
#define XC_FLAGS_HAVE_KXC         (1 <<  3) /*     8 */
#define XC_FLAGS_HAVE_LXC         (1 <<  4) /*    16 */
#define XC_FLAGS_1D               (1 <<  5) /*    32 */
#define XC_FLAGS_2D               (1 <<  6) /*    64 */
#define XC_FLAGS_3D               (1 <<  7) /*   128 */
#define XC_FLAGS_HYB_CAM          (1 <<  8) /*   256 */
#define XC_FLAGS_HYB_CAMY         (1 <<  9) /*   512 */
#define XC_FLAGS_VV10             (1 << 10) /*  1024 */
#define XC_FLAGS_HYB_LC           (1 << 11) /*  2048 */
#define XC_FLAGS_HYB_LCY          (1 << 12) /*  4096 */
#define XC_FLAGS_STABLE           (1 << 13) /*  8192 */
#define XC_FLAGS_DEVELOPMENT      (1 << 14) /* 16384 */
#define XC_FLAGS_NEEDS_LAPLACIAN  (1 << 15) /* 32768 */

typedef struct{
  int   number;   /* identifier number */
  int   kind;     /* XC_EXCHANGE, XC_CORRELATION, XC_EXCHANGE_CORRELATION, XC_KINETIC */

  char *name;     /* name of the functional, e.g. "PBE" */
  int   family;   /* type of the functional, e.g. XC_FAMILY_GGA */
  func_reference_type *refs[XC_MAX_REFERENCES];  /* index of the references */

  int   flags;    /* see above for a list of possible flags */

  double dens_threshold;

  /* this allows to have external parameters in the functional */
  int n_ext_params;
  const func_params_type *ext_params;
  void (*set_ext_params)(struct xc_func_type *p, const double *ext_params);

  ...
} xc_func_info_type;

For example, for the Slater exchange functional, this structure is defined as

const xc_func_info_type xc_func_info_lda_x = {
  XC_LDA_X,
  XC_EXCHANGE,
  "Slater exchange",
  XC_FAMILY_LDA,
  {&xc_ref_Dirac1930_376, &xc_ref_Bloch1929_545, NULL, NULL, NULL},
  XC_FLAGS_3D | XC_FLAGS_HAVE_EXC | XC_FLAGS_HAVE_VXC | XC_FLAGS_HAVE_FXC | XC_FLAGS_HAVE_KXC,
  1e-24,
  0, NULL, NULL,
  ...
};

The user of the library can access the information in the following way:

#include <stdio.h>
#include <xc.h>

int main()
{
  int ii;
  xc_func_type func;
  /* Initialize the functional */
  xc_func_init(&func, XC_GGA_X_B88, XC_UNPOLARIZED);
  /* Print out information */
  printf("The functional '%s' is ", func.info->name);
  switch (func.info->kind) {
  case (XC_EXCHANGE):
    printf("an exchange functional");
    break;
  case (XC_CORRELATION):
    printf("a correlation functional");
    break;
  case (XC_EXCHANGE_CORRELATION):
    printf("an exchange-correlation functional");
    break;
  case (XC_KINETIC):
    printf("a kinetic energy functional");
    break;
  default:
    printf("of unknown kind");
    break;
  }
  printf(", it belongs to the '");
  /* Print out family */
  switch (func.info->family) {
  case (XC_FAMILY_LDA):
    printf("LDA");
    break;
  case (XC_FAMILY_GGA):
    printf("GGA");
    break;
  case (XC_FAMILY_HYB_GGA):
    printf("Hybrid GGA");
    break;
  case (XC_FAMILY_MGGA):
    printf("MGGA");
    break;
  case (XC_FAMILY_HYB_MGGA):
    printf("Hybrid MGGA");
    break;
  default:
    printf("unknown");
    break;
  }
  printf("' family and is defined in the reference(s):\n");
  /* Print out references */
  for(ii = 0; func.info->refs[ii] != NULL; ii++){
    printf("[%d] %s\n", ii+1, func.info->refs[ii]->ref);
  }
  /* Deallocate memory */
  xc_func_end(&func);
}

The output of this example should look like this:

The functional 'Becke 88' is an exchange functional, it belongs to the 'GGA' family and is defined in the reference(s):
[1] A. D. Becke, Phys. Rev. A 38, 3098 (1988)

Several routines are available to access the information about the functional from Fortran. Here is the previous example using the Fortran 90 interface:

program xcinfo
  use xc_f03_lib_m
  implicit none
  type(xc_f03_func_t) :: xc_func
  type(xc_f03_func_info_t) :: xc_info
  type(xc_f03_func_reference_t) :: xc_ref
  integer :: i, number
  character(len=120) :: kind, family
  ! Initialize the functional
  call xc_f03_func_init(xc_func, XC_GGA_X_B88, XC_UNPOLARIZED)
  xc_info = xc_f03_func_get_info(xc_func)
  ! Get the type of the functional
  select case(xc_f03_func_info_get_kind(xc_info))
  case (XC_EXCHANGE)
    write(kind, '(a)') 'an exchange functional'
  case (XC_CORRELATION)
    write(kind, '(a)') 'a correlation functional'
  case (XC_EXCHANGE_CORRELATION)
    write(kind, '(a)') 'an exchange-correlation functional'
  case (XC_KINETIC)
    write(kind, '(a)') 'a kinetic energy functional'
  case default
    write(kind, '(a)') 'of unknown kind'
  end select
  ! Get the family
  select case (xc_f03_func_info_get_family(xc_info))
  case (XC_FAMILY_LDA);
    write(family,'(a)') "LDA"
  case (XC_FAMILY_GGA);
    write(family,'(a)') "GGA"
  case (XC_FAMILY_HYB_GGA);
    write(family,'(a)') "Hybrid GGA"
  case (XC_FAMILY_MGGA);
    write(family,'(a)') "MGGA"
  case (XC_FAMILY_HYB_MGGA);
    write(family,'(a)') "Hybrid MGGA"
  case default;
    write(family,'(a)') "unknown"
  end select
  ! Print out information
  write(*,'("The functional ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
    trim(xc_f03_func_info_get_name(xc_info)), trim(kind), trim(family)
  ! Print out references
  do i = 0, XC_MAX_REFERENCES-1
    ! get_references updates its index argument, setting it negative when
    ! there are no more references; use a copy since it is intent(inout)
    number = i
    xc_ref = xc_f03_func_info_get_references(xc_info, number)
    if (number < 0) exit
    write(*, '(a,i1,2a)') '[', i+1, '] ', trim(xc_f03_func_reference_get_ref(xc_ref))
  end do
  ! Free the memory
  call xc_f03_func_end(xc_func)
end program xcinfo

Initialization

A functional is initialized by the routine

int xc_func_init(xc_func_type *p, int functional, int nspin);

input:
  functional: which functional do we want?
  nspin: either XC_UNPOLARIZED or XC_POLARIZED
output:
  p: structure that holds our functional
returns: 0 (OK) or -1 (ERROR)

External parameters

There are some functionals that require extra information, like the number of particles or the temperature. The number of parameters and a short description of their meaning is included in the info structure

typedef struct{
  double value;
  char *description;
} func_params_type;

typedef struct{
  ...

  /* this allows to have external parameters in the functional */
  int n_ext_params;
  const func_params_type *ext_params;
  void (*set_ext_params)(struct xc_func_type *p, const double *ext_params);

  ...
} xc_func_info_type;

Note that several functions are provided in the Fortran interface to access this information.

These parameters are set by default to some value and changing them is done using the xc_func_set_ext_params function

xc_func_set_ext_params(xc_func_type *p, double *ext_params);

input:
  p: structure obtained from calling xc_func_init
  ext_params: array of external parameters

The ext_params array must contain a value for every parameter, in the order in which they are listed in the info structure. Individual parameters can also be set and queried by name, and the current values read back, with

void   xc_func_set_ext_params_name(xc_func_type *p, const char *name, double par);
double xc_func_get_ext_params_name(const xc_func_type *p, const char *name);
void   xc_func_get_ext_params(const xc_func_type *p, double *ext_params);

Evaluation

The routines that should be used depend on the functional family. Note: The routines expect a nonnegative density.

LDA

This is done by the routines

void xc_lda_exc(const xc_func_type *p, size_t np, const double *rho,
     double *exc);
void xc_lda_vxc(const xc_func_type *p, size_t np, const double *rho,
     double *vrho);
void xc_lda_exc_vxc(const xc_func_type *p, size_t np, const double *rho,
     double *exc, double *vrho);
void xc_lda_exc_vxc_fxc(const xc_func_type *p, size_t np, const double *rho,
     double *exc, double *vrho, double *v2rho2);
void xc_lda_fxc(const xc_func_type *p, size_t np, const double *rho,
     double *v2rho2);
void xc_lda_kxc(const xc_func_type *p, size_t np, const double *rho,
     double *v3rho3);
void xc_lda_lxc(const xc_func_type *p, size_t np, const double *rho,
     double *v4rho4);
void xc_lda(const xc_func_type *p, size_t np, const double *rho, double *exc,
     double *vrho, double *v2rho2, double *v3rho3, double *v4rho4);

input:
  p: structure obtained from calling xc_func_init
  np: number of points
  rho[]: the density
output:
  exc[]: the energy density per unit particle
  vrho[]: first derivative of the energy per unit volume
  v2rho2[]: second derivative of the energy per unit volume
  v3rho3[]: third derivative of the energy per unit volume
  v4rho4[]: fourth derivative of the energy per unit volume

The components of the density rho are

$$ \mathrm{rho}[0] = \rho_\uparrow \qquad \mathrm{rho}[1] = \rho_\downarrow $$

The total energy is given by

$$ E = \int \epsilon({\bf r}) {\rm d}^3 r = \int \rho({\bf r}) e({\bf r}) {\rm d}^3 r $$

where $e$ is the the energy density per unit particle, and the derivatives of the energy density $\epsilon$ are defined as

1st order $$ \mathrm{vrho}_\alpha \quad;\quad \frac{d \epsilon}{d \rho_\alpha} \quad;\quad [(0), (1)] $$

2nd order $$ \mathrm{v2rho2}_{\alpha\beta} \quad;\quad \frac{d^2 \epsilon}{d \rho_\alpha d \rho_\beta} \quad;\quad [(0,0), (0,1), (1,1)] $$

3rd order $$ \mathrm{v3rho3}_{\alpha\beta\gamma} \quad;\quad \frac{d^3 \epsilon}{d \rho_\alpha d \rho_\beta d \rho_\gamma} \quad;\quad [(0,0,0), (0,0,1), (0,1,1), (1,1,1)] $$

4th order $$ \mathrm{v4rho4}_{\alpha\beta\gamma\kappa} \quad;\quad \frac{d^4 \epsilon}{d \rho_\alpha d \rho_\beta d \rho_\gamma d \rho_\kappa} \quad;\quad [(0,0,0,0), (0,0,0,1), (0,0,1,1), (0,1,1,1), (1,1,1,1)] $$

If the functional was initialized with nspin=XC_UNPOLARIZED, the spin indices should be dropped from the previous expressions, resulting in arrays of size np. Otherwise, the parameters have dimensions rho[2*np], vrho[2*np], v2rho2[3*np], v3rho3[4*np], and v4rho4[5*np]. The energy per unit particle array always has dimensions exc[np].

GGA

This is done by the routines

void xc_gga_exc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *exc)
void xc_gga_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *vrho, double *vsigma)
void xc_gga_exc_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *exc, double *vrho, double *vsigma)
void xc_gga_exc_vxc_fxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *exc, double *vrho, double *vsigma, double *v2rho2, double *v2rhosigma, double *v2sigma2)
void xc_gga_fxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *v2rho2, double *v2rhosigma, double *v2sigma2)
void xc_gga_kxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3)
void xc_gga_lxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *v4rho4, double *v4rho3sigma, double *v4rho2sigma2, double *v4rhosigma3, double *v4sigma4)
void xc_gga(const xc_func_type *p, size_t np, const double *rho, const double *sigma,
     double *exc, double *vrho, double *vsigma,
     double *v2rho2, double *v2rhosigma, double *v2sigma2,
     double *v3rho3, double *v3rho2sigma, double *v3rhosigma2, double *v3sigma3,
     double *v4rho4, double *v4rho3sigma, double *v4rho2sigma2, double *v4rhosigma3, double *v4sigma4)

input:
  p: structure obtained from calling xc_func_init
  np: number of points
  rho[]: the density
  sigma[]: contracted gradients of the density
output:
  exc[]: the energy per unit particle
  vrho[]: first partial derivative of the energy per unit volume in terms of the density
  vsigma[]: first partial derivative of the energy per unit volume in terms of sigma
  v2rho2[]: second partial derivative of the energy per unit volume in terms of the density
  v2rhosigma[]: second partial derivative of the energy per unit volume in terms of the density and sigma
  v2sigma2[]: second partial derivative of the energy per unit volume in terms of sigma
  etc.

The components of the contracted gradient sigma are

$$ \sigma[0] = \nabla \rho_\uparrow \cdot \nabla \rho_\uparrow \qquad \sigma[1] = \nabla \rho_\uparrow \cdot \nabla \rho_\downarrow \qquad \sigma[2] = \nabla \rho_\downarrow \cdot \nabla \rho_\downarrow \qquad $$

The extra derivatives are defined as

1st order $$ \mathrm{vsigma}_{\alpha} \quad;\quad \frac{\partial \epsilon}{\partial \sigma_\alpha} \quad;\quad [(0), (1), (2)] $$

2nd order $$ \mathrm{v2rhosigma}_{\alpha\beta} \quad;\quad \frac{\partial^2 \epsilon}{\partial \rho_\alpha \partial \sigma_\beta} \quad;\quad [(0,0), (0,1), (0,2), (1,0), (1,1), (1,2)] $$ $$ \mathrm{v2sigma2}_{\alpha\beta} \quad;\quad \frac{d^2 \epsilon}{\partial \sigma_\alpha \partial \sigma_\beta} \quad;\quad [(0,0), (0,1), (0,2), (1,1), (1,2), (2,2)] $$

3rd order $$ \mathrm{v3rho2sigma}_{\alpha\beta\gamma} \quad;\quad \frac{d^3 \epsilon}{\partial \rho_\alpha \partial \rho_\beta \partial \sigma_\gamma} \quad;\quad [(0,0,0), (0,0,1), (0,0,2), (0,1,0), (0,1,1), (0,1,2), (1,1,0), (1,1,1), (1,1,2)] $$ $$ \mathrm{v3rhosigma2}_{\alpha\beta\gamma} \quad;\quad \frac{d^3 \epsilon}{\partial \rho_\alpha \partial \sigma_\beta \partial \sigma_\gamma} \quad;\quad [(0,0,0), (0,0,1), (0,0,2), (0,1,1), (0,1,2), (0,2,2) (1,0,0), (1,0,1), (1,0,2), (1,1,1), (1,1,2), (1,2,2)] $$ $$ \mathrm{v3sigma3}_{\alpha\beta\gamma} \quad;\quad \frac{d^3 \epsilon}{\partial \sigma_\alpha \partial \sigma_\beta \partial \sigma_\gamma} \quad;\quad [(0,0,0), (0,0,1), (0,0,2), (0,1,1), (0,1,2), (0,2,2), (1,1,1), (1,1,2), (1,2,2), (2,2,2)] $$

4th order $$ \mathrm{v4rho3sigma}_{\alpha\beta\gamma\kappa} \quad;\quad \frac{d^4 \epsilon}{\partial \rho_\alpha \partial \rho_\beta \partial \rho_\gamma \partial \sigma_\kappa} \quad;\quad [ (0,0,0,0), (0,0,0,1), (0,0,0,2), (0,0,1,0), (0,0,1,1), (0,0,1,2), (0,1,1,0), (0,1,1,1), (0,1,1,2), (1,1,1,0), (1,1,1,1), (1,1,1,2) ] $$ $$ \mathrm{v4rho2sigma2}_{\alpha\beta\gamma\kappa} \quad;\quad \frac{d^4 \epsilon}{\partial \rho_\alpha \partial \rho_\beta \partial \sigma_\gamma \partial \sigma_\kappa} \quad;\quad [ (0,0,0,0), (0,0,0,1), (0,0,0,2), (0,0,1,1), (0,0,1,2), (0,0,2,2), (0,1,0,0), (0,1,0,1), (0,1,0,2), (0,1,1,1), (0,1,1,2), (0,1,2,2), (1,1,0,0), (1,1,0,1), (1,1,0,2), (1,1,1,1), (1,1,1,2), (1,1,2,2) ] $$ $$ \mathrm{v4rhosigma3}_{\alpha\beta\gamma\kappa} \quad;\quad \frac{d^4 \epsilon}{\partial \rho_\alpha \partial \sigma_\beta \partial \sigma_\gamma \partial \sigma_\kappa} \quad;\quad [ (0,0,0,0), (0,0,0,1), (0,0,0,2), (0,0,1,1), (0,0,1,2), (0,0,2,2), (0,1,1,1), (0,1,1,2), (0,1,2,2), (0,2,2,2), (1,0,0,0), (1,0,0,1), (1,0,0,2), (1,0,1,1), (1,0,1,2), (1,0,2,2), (1,1,1,1), (1,1,1,2), (1,1,2,2), (1,2,2,2) ] $$ $$ \mathrm{v4sigma4}_{\alpha\beta\gamma\kappa} \quad;\quad \frac{d^4 \epsilon}{\partial \sigma_\alpha \partial \sigma_\beta \partial \sigma_\gamma \partial \sigma_\kappa} \quad;\quad [ (0,0,0,0), (0,0,0,1), (0,0,0,2), (0,0,1,1), (0,0,1,2), (0,0,2,2), (0,1,1,1), (0,1,1,2), (0,1,2,2), (0,2,2,2), (1,1,1,1), (1,1,1,2), (1,1,2,2), (1,2,2,2), (2,2,2,2) ] $$

MetaGGA

This is done by the routines

void xc_mgga_exc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *exc)
void xc_mgga_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *vrho, double *vsigma, double *vlapl, double *vtau)
void xc_mgga_exc_vxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double * exc, double *vrho, double *vsigma, double *vlapl, double *vtau)
void xc_mgga_fxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *v2rho2, double *v2rhosigma, double *v2rholapl, double *v2rhotau, double *v2sigma2,
     double *v2sigmalapl, double *v2sigmatau, double *v2lapl2, double *v2lapltau,  double *v2tau2)
void xc_mgga_kxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *v3rho3, double *v3rho2sigma, double *v3rho2lapl, double *v3rho2tau, double *v3rhosigma2,
     double *v3rhosigmalapl, double *v3rhosigmatau, double *v3rholapl2, double *v3rholapltau,
     double *v3rhotau2, double *v3sigma3, double *v3sigma2lapl, double *v3sigma2tau,
     double *v3sigmalapl2, double *v3sigmalapltau, double *v3sigmatau2, double *v3lapl3,
     double *v3lapl2tau, double *v3lapltau2, double *v3tau3)
void xc_mgga_lxc(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *v4rho4, double *v4rho3sigma, double *v4rho3lapl, double *v4rho3tau, double *v4rho2sigma2,
     double *v4rho2sigmalapl, double *v4rho2sigmatau, double *v4rho2lapl2, double *v4rho2lapltau,
     double *v4rho2tau2, double *v4rhosigma3, double *v4rhosigma2lapl, double *v4rhosigma2tau,
     double *v4rhosigmalapl2, double *v4rhosigmalapltau, double *v4rhosigmatau2,
     double *v4rholapl3, double *v4rholapl2tau, double *v4rholapltau2, double *v4rhotau3,
     double *v4sigma4, double *v4sigma3lapl, double *v4sigma3tau, double *v4sigma2lapl2,
     double *v4sigma2lapltau, double *v4sigma2tau2, double *v4sigmalapl3, double *v4sigmalapl2tau,
     double *v4sigmalapltau2, double *v4sigmatau3, double *v4lapl4, double *v4lapl3tau,
     double *v4lapl2tau2, double *v4lapltau3, double *v4tau4)
void xc_mgga(const xc_func_type *p, size_t np, const double *rho, const double *sigma, const double *lapl, const double *tau,
     double *exc, double *vrho, double *vsigma, double *vlapl, double *vtau,
     double *v2rho2, double *v2rhosigma, double *v2rholapl, double *v2rhotau, double *v2sigma2,
     double *v2sigmalapl, double *v2sigmatau, double *v2lapl2, double *v2lapltau, double *v2tau2,
     double *v3rho3, double *v3rho2sigma, double *v3rho2lapl, double *v3rho2tau, double *v3rhosigma2,
     double *v3rhosigmalapl, double *v3rhosigmatau, double *v3rholapl2, double *v3rholapltau,
     double *v3rhotau2, double *v3sigma3, double *v3sigma2lapl, double *v3sigma2tau,
     double *v3sigmalapl2, double *v3sigmalapltau, double *v3sigmatau2, double *v3lapl3,
     double *v3lapl2tau, double *v3lapltau2, double *v3tau3,
     double *v4rho4, double *v4rho3sigma, double *v4rho3lapl, double *v4rho3tau, double *v4rho2sigma2,
     double *v4rho2sigmalapl, double *v4rho2sigmatau, double *v4rho2lapl2, double *v4rho2lapltau,
     double *v4rho2tau2, double *v4rhosigma3, double *v4rhosigma2lapl, double *v4rhosigma2tau,
     double *v4rhosigmalapl2, double *v4rhosigmalapltau, double *v4rhosigmatau2,
     double *v4rholapl3, double *v4rholapl2tau, double *v4rholapltau2, double *v4rhotau3,
     double *v4sigma4, double *v4sigma3lapl, double *v4sigma3tau, double *v4sigma2lapl2,
     double *v4sigma2lapltau, double *v4sigma2tau2, double *v4sigmalapl3, double *v4sigmalapl2tau,
     double *v4sigmalapltau2, double *v4sigmatau3, double *v4lapl4, double *v4lapl3tau,
     double *v4lapl2tau2, double *v4lapltau3, double *v4tau4)

input:
  p: structure obtained from calling xc_func_init
  np: number of points
  rho[]: the density
  sigma[]: contracted gradients of the density
  lapl[]: the laplacian of the density
  tau[]: the kinetic energy density
output:
  exc[]: energy density per unit particle
  vrho[]: first partial derivative of the energy per unit volume in terms of the density
  vsigma[]: first partial derivative of the energy per unit volume in terms of sigma
  vlapl[]: first partial derivative of the energy per unit volume in terms of the laplacian of the density
  vtau[]: first partial derivative of the energy per unit volume in terms of the kinetic energy density
  etc.

Note that the kinetic energy density is defined with the 1/2 factor (until Libxc version 2.0.0, tau was defined as twice the kinetic energy density):

$$ \mathrm{lapl}[0] = \nabla^2 \rho_\uparrow \qquad \mathrm{lapl}[1] = \nabla^2 \rho_\downarrow $$ $$ {\tau}_{\alpha} = \frac{1}{2}\sum_i | \nabla \psi_{\alpha i} |^2 \qquad \mathrm{tau}[0] = {\tau}_\uparrow \qquad \mathrm{tau}[1] = {\tau}_\downarrow $$

Many meta-GGA functionals do not depend on the density Laplacian. A functional uses lapl only if the XC_FLAGS_NEEDS_LAPLACIAN flag is set in func.info->flags; otherwise the lapl input (and the corresponding vlapl and higher Laplacian derivatives) is ignored and may be passed as NULL.

The extra derivatives are defined as

1st order $$ \mathrm{vlapl}_{\alpha} \quad;\quad \frac{\partial \epsilon}{\partial \nabla^2 \rho_\alpha} \quad;\quad [(0), (1)] $$ $$ \mathrm{vtau}_{\alpha} \quad;\quad \frac{\partial \epsilon}{\partial \tau_\alpha} \quad;\quad [(0), (1)] $$

By now you should have understood how the higher derivatives (2nd, 3rd, and 4th) work. Just recall that the order of variables is $[\rho,\sigma,\nabla^2 \rho, \tau]$.

Array dimensions

The number of components stored for a given output array is fixed by the following rule. For each variable that appears in the array name, count the number of distinct symmetric multi-indices for its spin components: a variable with $s$ spin components differentiated $k$ times contributes

$$ \binom{s+k-1}{k} $$

components, and the total size of the array (the number of components per grid point) is the product of these counts over all variables in the name. In the spin-polarized case, $\rho$, $\nabla^2\rho$ and $\tau$ have $s=2$ while $\sigma$ has $s=3$; in the spin-unpolarized case every count is $1$, so all arrays have a single component per point. The full array length is this count multiplied by np.

For example, v4rho2sigma2 differentiates twice with respect to $\rho$ ($s=2$, $k=2 \Rightarrow 3$) and twice with respect to $\sigma$ ($s=3$, $k=2 \Rightarrow 6$), giving $3\times 6 = 18$ components in the polarized case.

The per-point component counts for the LDA and GGA arrays are:

The meta-GGA arrays involving $\nabla^2\rho$ and $\tau$ follow the same rule; for instance v4sigma2lapltau has $6\times 3\times 2 = 36$ components in the polarized case.

Any derivative output for which you do not want the result may be passed as a NULL pointer, and the corresponding quantity will not be computed. Before requesting a given derivative order, you can check the relevant XC_FLAGS_HAVE_EXC, XC_FLAGS_HAVE_VXC, XC_FLAGS_HAVE_FXC, XC_FLAGS_HAVE_KXC, or XC_FLAGS_HAVE_LXC flag in func.info->flags to verify that the functional supports it.

In addition to the individual routines above and the all-in-one xc_lda, xc_gga, and xc_mgga evaluators, Libxc also provides combined routines that return a contiguous range of derivative orders, e.g.

void xc_lda_exc_vxc_fxc(const xc_func_type *p, size_t np, const double *rho,
     double *exc, double *vrho, double *v2rho2);
void xc_gga_exc_vxc_fxc(const xc_func_type *p, size_t np, const double *rho,
     const double *sigma, double *exc, double *vrho, double *vsigma,
     double *v2rho2, double *v2rhosigma, double *v2sigma2);

with the analogous xc_*_vxc_fxc and xc_*_vxc_fxc_kxc variants for each family.

Accessor functions

For hybrid and non-local-correlation functionals, Libxc evaluates only the semilocal part of the functional. The exact-exchange (Hartree–Fock) admixture of a hybrid and the non-local van der Waals correlation (VV10) of a non-local functional are not computed by Libxc and must be added by the calling program. The flags XC_FLAGS_HYB_CAM, XC_FLAGS_HYB_CAMY, XC_FLAGS_HYB_LC, XC_FLAGS_HYB_LCY, and XC_FLAGS_VV10 indicate which contributions a given functional expects, and the accessor functions below return the parameters needed to construct them.

Some functionals provide accessor functions that give parameters needed for calculating quantities outside of libxc:

The mixing parameter for global hybrids:

double xc_hyb_exx_coef(const xc_func_type *p);

Parameters for CAM-type hybrids, where omega is the range-separation parameter, alpha is the fraction of Hartree-Fock exchange, and beta is the fraction of short-range exchange:

void xc_hyb_cam_coef(const xc_func_type *p, double *omega, double *alpha, double *beta);

Non-local correlation parameters for functionals such as WB97X_V, VV10, LC_VV10, and WB97M_V:

void  xc_nlc_coef(const xc_func_type *p, double *nlc_b, double *nlc_C);

Asymptotic parameter for AK13:

double xc_gga_ak13_get_asymptotic(double homo);

Destruction

If you no longer need your xc functional, you should destroy it with the routine

void xc_func_end (xc_func_type *p);

input:
  p: structure holding the functional to destroy

Other useful functions

xc_family_from_id

If you have the identifier of a functional (that was read, for example from an input file), and want to find out to which family it belongs to, you can use the routine

int xc_family_from_id(int functional);

input:
  functional: identifier of the functional
returns:
  XC_FAMILY_UNKNOWN: could not find the family
  XC_FAMILY_LDA, XC_FAMILY_GGA, etc.

Analogously for Fortran:

  integer function xc_family_from_id (functional)

xc_functional_get_name

If you know the identifier of a functional and want to know its name without having to call xc_func_init, you can use the function

char *xc_functional_get_name(int number);

input
  number: identifier of the functional
returns:
  The name of the functional

xc_functional_get_number

If you know the name of a functional and want to know its identifier, you can use the function

int xc_functional_get_number(const char *name);

input
  name: the name of the functional in the form "xc_gga_x_pbe" or "gga_x_pbe". Note that the name is case insensitive.
returns:
  The functional identifier

xc_number_of_functionals

This function returns the total number of functionals available in the library. It is particularly useful if one wants to construct a list of all the available functional names and/or identification numbers (see xc_available_functional_names and xc_available_functional_names).

int xc_number_of_functionals();

returns:
  Number of functionals

xc_maximum_name_length

This function returns the maximum string length of all the functional names available in the library. It is particularly useful if one wants to construct a list of all the available functional names (see xc_available_functional_names).

int xc_maximum_name_length();

returns:
  Maximum name length.

xc_available_functional_numbers

Returns a list of all known functional identifiers.

void xc_available_functional_numbers(int *list);

output:
  list: list of all the known functional identifiers.

xc_available_functional_names

Returns a list of all known functional names.

void xc_available_functional_names(char **list);

output:
  list: list of all the known functional names.

Density threshold

When the density is smaller than a certain threshold, the library will skip the evaluation of the functional and instead set all the output quantities to zero. This is to avoid underflows and/or overflows. Each functional is defined with a default threshold that should be adquate for most cases, but the range of densities in which a functional can be safely evaluated might depend on the computer architecture or the compiler and compiler flags used. Therefore, the library provides a function to change the threshold

void xc_func_set_dens_threshold(xc_func_type *p, double dens_threshold);

input
  p: structure holding the functional
  dens_threshold: the new density threshold

Analogous setters are available for the other reduced variables, which can be used to safeguard the evaluation of GGA and meta-GGA functionals:

void xc_func_set_zeta_threshold (xc_func_type *p, double t_zeta);
void xc_func_set_sigma_threshold(xc_func_type *p, double t_sigma);
void xc_func_set_tau_threshold  (xc_func_type *p, double t_tau);

In order to help the user choose an appropriate value for the density threshold, the xc-threshold utility is provided to evaluate a functional for several different trial densities.

Version

The actual version of libxc being used can be obtained by calling

void xc_version (int *major, int *minor, int *micro);

output:
  major: the major version number
  minor: the minor version number
  micro: the micro version number

This information is also available in the xc_version.h header file, where XC_VERSION, XC_MAJOR_VERSION, XC_MINOR_VERSION, and XC_MICRO_VERSION symbols are defined.

Linking via GNU autotools

An m4 macro for use in a configure script for a project linking to the Fortran 90 Libxc interface is available in the distribution, taken from Octopus.