Libxc:manual

From OctopusWiki

Jump to: navigation, search

Contents

Compatibility notice

  • lda_x_set_params: libxc v1.3 uses now three parameters instead of two, the last one specifies the omega used for range separation.

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

#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 zk[5];
  int i, func_id = 1;

  if(xc_func_init(&func, func_id, XC_UNPOLARIZED) != 0){
    fprintf(stderr, "Functional '%d' not found\n", func_id);
    return 1;
  }

  switch(func.info->family)
    {
    case XC_FAMILY_LDA:
      xc_lda_exc(&func, 5, rho, zk);
      break;
    case XC_FAMILY_GGA:
    case XC_FAMILY_HYB_GGA:
      xc_gga_exc(&func, 5, rho, sigma, zk);
      break;
    }

  for(i=0; i<5; i+=1){
    printf("%lf %lf\n", rho[i], zk[i]);
  }

  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 unit volume (zk). Finally the functions xc_func_end cleans up.


Fortran 90 bindings are also included in libxc. These can be found in the file libxc_master.F90 (show). In general, calling libxc from Fortran is as simple as from C. Here is the previous example in Fortran:

program lxctest
   use xc_f90_types_m
   use xc_f90_lib_m

   implicit none

   TYPE(xc_f90_pointer_t) :: xc_func
   TYPE(xc_f90_pointer_t) :: xc_info
   real(8) :: rho(5) = (/0.1, 0.2, 0.3, 0.4, 0.5/)
   real(8) :: sigma(5) = (/0.2, 0.3, 0.4, 0.5, 0.6/)
   real(8) :: zk(5)
   integer :: i, func_id = 1

   call xc_f90_func_init(xc_func, xc_info, func_id, XC_UNPOLARIZED)

   select case (xc_f90_info_family(xc_info))
   case(XC_FAMILY_LDA)
     call xc_f90_lda_exc(xc_func, 5, rho(1), zk(1))
   case(XC_FAMILY_GGA, XC_FAMILY_HYB_GGA)
     call xc_f90_gga_exc(xc_func, 5, rho(1), sigma(1), zk(1))
   end select

   do i = 1, 5
     write(*,"(F8.6,1X,F9.6)") rho(i), zk(i)
   end do

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

/* flags that can be used in info.flags */
#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_STABLE           (1 <<  9) /*  512 */
#define XC_FLAGS_DEVELOPMENT      (1 << 10) /* 1024 */

typedef struct{
  int   number;   /* indentifier number */
  int   kind;     /* XC_EXCHANGE or XC_CORRELATION */

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

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

  ...
} 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,
  "PAM Dirac, Proceedings of the Cambridge Philosophical Society 26, 376 (1930)"
  "F Bloch, Zeitschrift fuer Physik 57, 545 (1929)",
  XC_FLAGS_3D | XC_FLAGS_HAVE_EXC | XC_FLAGS_HAVE_VXC | XC_FLAGS_HAVE_FXC | XC_FLAGS_HAVE_KXC, 

  ...
};

Note that the references are separated by a newline. The user of the library can access the information in the following way:

#include <stdio.h>

#include <xc.h>

int main()
{
  xc_func_type func;

  xc_func_init(&func, XC_GGA_X_B88, XC_UNPOLARIZED);

  printf("The functional '%s' is defined in the reference(s):\n%s\n", func.info->name, func.info->refs);

  xc_func_end(&func);
}


Several routines are available to access the information about the functional from Fortran. You can take a look at the following example for help

program lxctest
   use xc_f90_types_m
   use xc_f90_lib_m

   implicit none

   TYPE(xc_f90_pointer_t) :: xc_func
   TYPE(xc_f90_pointer_t) :: xc_info
   integer :: i
   character(len=120) :: s1, s2
   type(xc_f90_pointer_t) :: str

   call xc_f90_func_init(xc_func, xc_info, XC_GGA_C_PBE, XC_UNPOLARIZED)

   select case(xc_f90_info_kind(xc_info))
   case(XC_EXCHANGE)
     write(*, '(a)') 'Exchange'
   case(XC_CORRELATION)
     write(*, '(a)') 'Correlation'
   case(XC_EXCHANGE_CORRELATION)
     write(*, '(a)') 'Exchange-correlation'
   end select

   call xc_f90_info_name(xc_info, s1)
   select case(xc_f90_info_family(xc_info))
   case (XC_FAMILY_LDA);      write(s2,'(a)') "LDA"
   case (XC_FAMILY_GGA);      write(s2,'(a)') "GGA"
   case (XC_FAMILY_HYB_GGA);  write(s2,'(a)') "Hybrid GGA"
   case (XC_FAMILY_MGGA);     write(s2,'(a)') "MGGA"
   case (XC_FAMILY_LCA);      write(s2,'(a)') "LCA"
   end select
   write(*, '(4a)') trim(s1), ' (', trim(s2), ')'
      
   i = 0
   call xc_f90_info_refs(xc_info, i, str, s1)
   do while(i >= 0)
     write(*, '(a,i1,2a)') '[', i, '] ', trim(s1)
     call xc_f90_info_refs(xc_info, i, str, s1)
   end do

   call xc_f90_func_end(xc_func)

end program lxctest

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)

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)

There are some functionals that require extra information. These parameters are set by default to some value and changing them is done using the xc_XXX_set_params routines. For example, by default the non-relativistic LDA exchange is used. To use the relativistic version, one should do:

 xc_lda_x_set_params(&func,  XC_RELATIVISTIC, 0.0);

[The last number is omega (ω) which is used for range-separated LDAs, e.g. in the the HSE functionals; for standard LDA, omega is 0.]

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(xc_lda_type *p, int np, double *rho, double *exc);
void xc_lda_vxc(xc_lda_type *p, int np, double *rho, double *vxc);
void xc_lda_vxc_exc(xc_lda_type *p, inp np, double *rho, double *exc, double *vxc);
void xc_lda_fxc(xc_lda_type *p, int np, double *rho, double *fxc);
void xc_lda_kxc(xc_lda_type *p, int np, double *rho, double *kxc);
void xc_lda(xc_lda_type *p, int np, double *rho, double *exc, double *vxc, double *fxc, double *kxc);
input:
  p: structure obtained from calling xc_func_init
  np: number of points
  rho[]: the density
output:
  exc: the energy per unit particle
  vxc[]: first derivative of the energy per unit volume
  fxc[]: second derivative of the energy per unit volume
  kxc[]: third derivative of the energy per unit volume

The derivatives are defined as


  {\rm vxc}_\alpha = \frac{d \epsilon}{d n_\alpha}


  {\rm fxc}_{\alpha\beta} = \frac{d^2 \epsilon}{d n_\alpha d n_\beta}


  {\rm kxc}_{\alpha\beta} = \frac{d^3 \epsilon}{d n_\alpha d n_\beta d n_\gamma}

where ε is the energy per unit volume.

If the functional was initialized with nspin=XC_UNPOLARIZED, the spin indices should be dropped from the previous expressions. Otherwise, the parameters have dimensions rho[2], vxc[2], fxc[3] and kxc[4]. The components of fxc are


  f_{\rm xc}[1] = f_{{\rm xc} \uparrow\uparrow}
  \qquad
  f_{\rm xc}[2] = f_{{\rm xc} \uparrow\downarrow}
  \qquad
  f_{\rm xc}[3] = f_{{\rm xc} \downarrow\downarrow}
  \qquad

and the ones of kxc are


  k_{\rm xc}[1] = k_{{\rm xc} \uparrow\uparrow\uparrow}
  \qquad
  k_{\rm xc}[2] = k_{{\rm xc} \uparrow\uparrow\downarrow}
  \qquad
  k_{\rm xc}[3] = k_{{\rm xc} \uparrow\downarrow\downarrow}
  \qquad
  k_{\rm xc}[4] = k_{{\rm xc} \downarrow\downarrow\downarrow}
  \qquad


GGA

This is done by the routines


void xc_gga_exc(xc_func_type *p, int np, double *rho, double *sigma, double *exc)
void xc_gga_vxc(xc_func_type *p, int np, double *rho, double *sigma, double double *vrho, double *vsigma)
void xc_gga_exc_vxc(xc_func_type *p, int np, double *rho, double *sigma, double *exc, double *vrho, double *vsigma)
void xc_gga_fxc(xc_func_type *p, int np, double *rho, double *sigma, double *v2rho2, double *v2rhosigma, double *v2sigma2)
void xc_gga(xc_func_type *p, int np, double *rho, double *sigma, double *exc, double *vrho, double *vsigma, double *v2rho2, double *v2rhosigma, double *v2sigma2)
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

The derivatives are defined as


  {\rm vrho}_{\alpha} = \frac{\partial \epsilon}{\partial n_\alpha}
  \qquad
  {\rm vsigma}_{\alpha} = \frac{\partial \epsilon}{\partial \sigma_\alpha}


  {\rm v2rho2}_{\alpha\beta} = \frac{d^2 \epsilon}{\partial n_\alpha \partial n_\beta}
  \qquad
  {\rm v2rhosigma}_{\alpha\beta} = \frac{\partial \epsilon}{\partial n_\alpha \partial \sigma_\beta}
  \qquad
  {\rm v2sigma2}_{\alpha\beta} = \frac{d^2 \epsilon}{\partial \sigma_\alpha \partial \sigma_\beta}


where ε is the energy per unit volume.

If the functional was initialized with nspin=XC_UNPOLARIZED, the spin indices should be dropped from the previous expressions. Otherwise, the parameters have dimensions rho[2], vxc[2], sigma[3], vsigma[3], v2rho2[3], v2rhosigma[6], v2sigma2[6]. The components of sigma are


  \sigma[1] = \nabla n_\uparrow \cdot \nabla n_\uparrow
  \qquad
  \sigma[2] = \nabla n_\uparrow \cdot \nabla n_\downarrow
  \qquad
  \sigma[3] = \nabla n_\downarrow \cdot \nabla n_\downarrow
  \qquad

Hybrid GGA

In this case the evaluation of the functional is done in the same way as the GGAs. Furthermore, there is a routine that returns the mixing parameter of the hybrid

double xc_hyb_gga_exx_coef(xc_func_type *p);

MetaGGA

This is done by the routines

void xc_mgga_exc(xc_func_type *p, int np, double *rho, double *sigma, double *lapl_rho, double *tau, double *zk)
void xc_mgga_vxc(xc_func_type *p, int np, double *rho, double *sigma, double *lapl_rho, double *tau, double *vrho, double *vsigma, double *vlapl_rho, double *vtau)
void xc_mgga_exc_vxc(xc_func_type *p, int np, double *rho, double *sigma, double *lapl_rho, double *tau, 
                     double * zk, double *vrho, double *vsigma, double *vlapl_rho, double *vtau)
void xc_mgga_fxc(xc_func_type *p, int np, double *rho, double *sigma, double *lapl_rho, double *tau,
                 double *v2rho2, double *v2rhosigma, double *v2sigma2, double *v2rhotau, double *v2tausigma, double *v2tau2)
void xc_mgga(xc_func_type *p, int np, double *rho, double *sigma, double *lapl_rho, double *tau,
             double *zk, double *vrho, double *vsigma, double *vlapl_rho, double *vtau,
             double *v2rho2, double *v2rhosigma, double *v2sigma2, double *v2rhotau, double *v2tausigma, double *v2tau2)
input:
  p: structure obtained from calling xc_func_init
  np: number of points
  rho[]: the density
  sigma[]: contracted gradients of the density
  lapl_rho[]: the laplacian of the density
  tau[]: the kinetic energy density
output:
  zk: 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_rho[]: 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

The derivatives are defined as


  {\rm vrho}_{\alpha} = \frac{\partial \epsilon}{\partial n_\alpha}
  \qquad
  {\rm vsigma}_{\alpha} = \frac{\partial \epsilon}{\partial \sigma_\alpha}
  \qquad
  {\rm vlapl\_rho}_{\alpha} = \frac{\partial \epsilon}{\partial \nabla^2 n_\alpha}
  \qquad
  {\rm vtau}_{\alpha} = \frac{\partial \epsilon}{\partial \tau_\alpha}

Note that the kinetic energy density is defined without the 1/2 factor


  {\tau}_{\alpha} = \sum_i | \nabla \psi_i |^2

Destruction

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

void xc_func_end (xc_func_type *p);
input:
  p: structure holding the functional to destroy

Available functionals

This is a list of the functionals the library currently knows about.

LDA Functionals

LDA Exchange
LDA Correlation
LDA Exchange-Correlation
LDA Kinetic Energy

GGA Functionals

GGA Exchange
GGA Correlation
GGA Exchange-Correlation
GGA Kinetic Energy

(For a review of lambda and gamma, see: E. V. Ludeña and V. V. Karasiev, Kinetic energy functionals: history, challenges and prospects. In: Reviews of Modern Quantum Chemistry: A Celebration of the Contributions of Robert Parr, Vol. 1, edited by K. D. Sen (World Scientific, Singapore, 2002), pp. 612–665.)

Hybrid Functionals

Hybrid Exchange

Hybrid Exchange-Correlation

TODO

  • XC_HYB_GGA_B97_1_REV

MetaGGA Functionals

MetaGGA Exchange
MetaGGA Correlation

Back to libxc

Personal tools