Libxc:manual

From OctopusWiki

Jump to: navigation, search

Contents

[edit]

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
[edit]

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
[edit]

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)
[edit]

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);
[edit]

Evaluation

The routines that should be used depend on the functional family.

[edit]

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


[edit]

GGA

This is done by the routines 


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

[edit]

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);
[edit]

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


[edit]

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
[edit]

Available functionals

This is a list of the functionals the library currently knows about. Exc, Vxc, Fxc and Kxc implementation audit status.

[edit]

LDA Functionals

[edit]
LDA Exchange
[edit]
LDA Correlation
[edit]
LDA Exchange-Correlation
[edit]

GGA Functionals

[edit]
GGA Exchange
[edit]
GGA Correlation
[edit]
GGA Exchange-Correlation
[edit]

Hybrid Functionals

TODO

  • XC_HYB_GGA_B97_1_REV
[edit]

MetaGGA Functionals

[edit]
MetaGGA Exchange
[edit]
MetaGGA Correlation
[edit]

LCA Functionals


Back to libxc

Retrieved from "http://www.tddft.org/programs/octopus/wiki/index.php/Libxc:manual"
Personal tools