XClib - libxc initialization

PW/src/init_run.f90

@@ -32,6 +32,8 @@ SUBROUTINE init_run()
   USE esm,                ONLY : do_comp_esm, esm_init
   USE tsvdw_module,       ONLY : tsvdw_initialize
   USE Coul_cut_2D,        ONLY : do_cutoff_2D, cutoff_fact 
+  USE lsda_mod,           ONLY : nspin
+  USE xc_lib,             ONLY : xclib_dft_is_libxc, xclib_init_libxc
   !
   IMPLICIT NONE
   !
@@ -109,6 +111,8 @@ SUBROUTINE init_run()
   CALL allocate_wfc_k()
   CALL openfil()
   !
+  IF (xclib_dft_is_libxc('ANY')) CALL xclib_init_libxc( nspin )
+  !
   CALL hinit0()
   !
   CALL potinit()

PW/src/make.depend

@@ -1012,6 +1012,7 @@ init_run.o : ../../Modules/recvec_subs.o
 init_run.o : ../../Modules/tsvdw.o
 init_run.o : ../../Modules/wannier_new.o
 init_run.o : ../../upflib/uspp.o
+init_run.o : ../../XClib/xc_lib.o
 init_run.o : Coul_cut_2D.o
 init_run.o : bp_mod.o
 init_run.o : dynamics_module.o
@@ -1791,6 +1792,7 @@ read_file_new.o : ../../Modules/recvec.o
 read_file_new.o : ../../Modules/recvec_subs.o
 read_file_new.o : ../../Modules/wavefunctions.o
 read_file_new.o : ../../upflib/uspp.o
+read_file_new.o : ../../XClib/xc_lib.o
 read_file_new.o : Coul_cut_2D.o
 read_file_new.o : buffers.o
 read_file_new.o : esm.o

PW/src/read_file_new.f90

@@ -71,9 +71,10 @@ SUBROUTINE read_file_new ( needwf )
   USE gvect,          ONLY : ngm, g
   USE gvecw,          ONLY : gcutw
   USE klist,          ONLY : nkstot, nks, xk, wk
-  USE lsda_mod,       ONLY : isk
+  USE lsda_mod,       ONLY : isk, nspin
   USE wvfct,          ONLY : nbnd, et, wg
   USE pw_restart_new, ONLY : read_xml_file
+  USE xc_lib,         ONLY : xclib_dft_is_libxc, xclib_init_libxc
   !
   IMPLICIT NONE
   !
@@ -88,6 +89,10 @@ SUBROUTINE read_file_new ( needwf )
   !
   CALL read_xml_file ( wfc_is_collected )
   !
+  ! ... initialize Libxc if needed
+  !
+  IF (xclib_dft_is_libxc('ANY')) CALL xclib_init_libxc( nspin )
+  !
   ! ... more initializations: pseudopotentials / G-vectors / FFT arrays /
   ! ... charge density / potential / ... , but not KS orbitals
   !

XClib/Makefile

@@ -25,8 +25,7 @@ xc_wrapper_d_lda_lsda.o \
 xc_wrapper_gga.o \
 xc_wrapper_lda_lsda.o \
 xc_wrapper_mgga.o \
-xc_lib.o \
+xc_lib.o
-xclib_test.o
 
 BEEF = \
 pbecor.o \
@@ -40,7 +39,7 @@ xc_lib.a: $(XCL) $(BEEF)
 
 xclib_test.x : xclib_test.o xc_lib.a 
 	$(LD) $(LD_LIBS) $(LDFLAGS) $(DFLAGS) -o $@ \
-		xclib_test.o xc_lib.a $(BLAS_LIBS)
+		xclib_test.o xc_lib.a $(BLAS_LIBS) $(LD_LIBS) $(LIBXC_LIBS)
 	- ( cd ../bin ; ln -fs ../XClib/$@ . )
 
 clean :

XClib/dft_mod.f90

@@ -26,9 +26,10 @@ MODULE dft_mod
             xclib_get_dft_short, xclib_get_dft_long,  &
             xclib_get_exx_fraction, xclib_get_finite_size_cell_volume, &
             get_screening_parameter, get_gau_parameter
-  PUBLIC :: xclib_dft_is, xclib_dft_is_libxc,  &
+  PUBLIC :: xclib_dft_is, xclib_dft_is_libxc, xclib_init_libxc,  &
             start_exx, stop_exx, dft_has_finite_size_correction, &
-            exx_is_active, igcc_is_lyp, xclib_reset_dft, dft_force_hybrid
+            exx_is_active, igcc_is_lyp, xclib_reset_dft, dft_force_hybrid, &
+            xclib_finalize_libxc
 #if defined(__LIBXC)
   PUBLIC :: get_libxc_flags_exc
 #endif
@@ -1013,17 +1014,20 @@ CONTAINS
      LOGICAL :: xclib_dft_is_libxc
      CHARACTER(len=*), INTENT(IN) :: family
      !! LDA, GGA or MGGA
-     CHARACTER(len=*), INTENT(IN) :: kindf
+     CHARACTER(len=*), INTENT(IN), OPTIONAL :: kindf
      !! EXCH or CORR
      !
-     CHARACTER(len=4) :: cfamily, ckindf
+     CHARACTER(len=4) :: cfamily='', ckindf
      INTEGER :: i, ln
      !
+     xclib_dft_is_libxc = .FALSE.
+     !
      ln = LEN_TRIM(family)
      !
      DO i = 1, ln
        cfamily(i:i) = capital(family(i:i))
      ENDDO
+     IF ( PRESENT(kindf) ) THEN
        DO i = 1, 4
          ckindf(i:i) = capital(kindf(i:i))
        ENDDO
@@ -1041,6 +1045,9 @@ CONTAINS
        CASE DEFAULT
          CALL xclib_error( 'xclib_dft_is_libxc', 'input not recognized', 1 )
        END SELECT
+     ELSE
+       IF (TRIM(cfamily)=='ANY'.AND.ANY(is_libxc(:))) xclib_dft_is_libxc=.TRUE.
+     ENDIF
      !
      RETURN
      !
@@ -1213,6 +1220,71 @@ CONTAINS
      IF (is_present) volume = finite_size_cell_volume
   END SUBROUTINE xclib_get_finite_size_cell_volume
   !
+  !--------------------------------------------------------------------------
+  SUBROUTINE xclib_init_libxc( xclib_nspin )
+    !------------------------------------------------------------------------
+    !! Initialize Libxc functionals, if present.
+    USE dft_par_mod,  ONLY: iexch, icorr, igcx, igcc, imeta, imetac, &
+                            is_libxc, libxc_initialized
+#if defined(__LIBXC)
+    USE dft_par_mod,  ONLY: n_ext_params, xc_func, xc_info, par_list
+#endif
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: xclib_nspin
+    !! 1: unpolarized case; 2: polarized
+    INTEGER :: iid, ip
+    INTEGER :: id_vec(6)
+    !
+#if defined(__LIBXC)
+    id_vec(1)=iexch ; id_vec(2)=icorr
+    id_vec(3)=igcx  ; id_vec(4)=igcc
+    id_vec(5)=imeta ; id_vec(6)=imetac
+    !
+    DO iid = 1, 6
+      IF (libxc_initialized(iid)) THEN
+        CALL xc_f03_func_end( xc_func(iid) )
+        libxc_initialized(iid) = .FALSE.
+      ENDIF
+      IF (is_libxc(iid)) THEN
+        CALL xc_f03_func_init( xc_func(iid), id_vec(iid), xclib_nspin )
+        xc_info(iid) = xc_f03_func_get_info( xc_func(iid) )
+        n_ext_params(iid) = xc_f03_func_info_get_n_ext_params( xc_info(iid) )
+        DO ip = 1, n_ext_params(iid)
+          par_list(iid,ip) = xc_f03_func_info_get_ext_params_default_value( &
+                                                           xc_info(iid), ip )
+        ENDDO
+        libxc_initialized(iid) = .TRUE.
+      ENDIF
+    ENDDO
+#endif
+    RETURN
+  END SUBROUTINE xclib_init_libxc
+  !
+  !--------------------------------------------------------------------------
+  SUBROUTINE xclib_finalize_libxc()
+    !------------------------------------------------------------------------
+    !! Finalize Libxc functionals, if present.
+    USE dft_par_mod,  ONLY: iexch, icorr, igcx, igcc, imeta, imetac, &
+                            is_libxc
+#if defined(__LIBXC)
+    USE dft_par_mod,  ONLY: xc_func
+#endif
+    IMPLICIT NONE
+    INTEGER :: iid
+    INTEGER :: id_vec(6)
+    !
+#if defined(__LIBXC)
+    id_vec(1)=iexch ; id_vec(2)=icorr
+    id_vec(3)=igcx  ; id_vec(4)=igcc
+    id_vec(5)=imeta ; id_vec(6)=imetac
+    !
+    DO iid = 1, 6
+      IF (is_libxc(iid)) CALL xc_f03_func_end( xc_func(iid) )
+    ENDDO
+#endif
+    RETURN
+  END SUBROUTINE xclib_finalize_libxc
+  !  
 #if defined(__LIBXC)
   !------------------------------------------------------------------------
   SUBROUTINE get_libxc_flags_exc( xc_info, eflag )

XClib/xc_input_params_mod.f90

@@ -11,6 +11,9 @@ MODULE dft_par_mod
     !! Parameters that define the XC functionals.
     !
     USE kind_l, ONLY: DP
+#if defined(__LIBXC)
+    USE xc_f03_lib_m
+#endif
     !
     IMPLICIT NONE
     !
@@ -22,6 +25,20 @@ MODULE dft_par_mod
     LOGICAL :: is_libxc(6) = .FALSE.
     !! \(\text{is_libxc(i)}=TRUE\) if the i-th term of the input 
     !! functional is from Libxc
+    !
+    LOGICAL :: libxc_initialized(6) = .FALSE.
+    !! TRUE if libxc functionals have been initialized
+    !
+#if defined(__LIBXC)
+    TYPE(xc_f03_func_t) :: xc_func(6)
+    !! pointers to libxc functional structs
+    TYPE(xc_f03_func_info_t) :: xc_info(6)
+    !! pointers to libxc info structs
+    INTEGER  :: n_ext_params(6) = 0._DP
+    !! number of external parameters for each functional
+    REAL(DP) :: par_list(6,10)
+    !! list of external parameters
+#endif
     !
     LOGICAL  :: exx_started = .FALSE.
     !! TRUE if Exact Exchange is active
@@ -69,7 +86,7 @@ MODULE dft_par_mod
     LOGICAL :: has_finite_size_correction = .FALSE.
     !! TRUE if finite size correction is present
     LOGICAL :: finite_size_cell_volume_set = .FALSE.
-    !! TRUE if the cell volume hase been set for finite size correction.
+    !! TRUE if the cell volume has been set for finite size correction.
     LOGICAL :: ismeta      = .FALSE.
     !! TRUE if the functional is MGGA
     LOGICAL :: ishybrid    = .FALSE.

XClib/xc_lib.f90

@@ -45,6 +45,8 @@ MODULE xc_lib
   !
   PUBLIC :: xclib_dft_is,                      &
             xclib_dft_is_libxc,                &
+            xclib_init_libxc,                  &
+            xclib_finalize_libxc,              &
             start_exx, stop_exx,               &
             dft_has_finite_size_correction,    &
             exx_is_active,                     &

XClib/xc_wrapper_d_gga.f90

@@ -23,6 +23,7 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
 #if defined(__LIBXC)
 #include "xc_version.h"
   USE xc_f03_lib_m
+  USE dft_par_mod,      ONLY: xc_func, xc_info
 #endif
   !
   IMPLICIT NONE
@@ -44,8 +45,6 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
   REAL(DP), ALLOCATABLE :: vrrc(:,:), vsrc(:,:), vssc(:), vrzc(:,:) 
   !
 #if defined(__LIBXC)
-  TYPE(xc_f03_func_t) :: xc_func
-  TYPE(xc_f03_func_info_t) :: xc_info2
   INTEGER :: fkind
   REAL(DP), ALLOCATABLE :: rho_lbxc(:)
   REAL(DP), ALLOCATABLE :: v2rho2_x(:), v2rhosigma_x(:), v2sigma2_x(:) 
@@ -119,10 +118,8 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
     ! ... DERIVATIVE FOR EXCHANGE
     v2rho2_x = 0._DP ;  v2rhosigma_x = 0._DP ;  v2sigma2_x = 0._DP
     IF (igcx /= 0) THEN 
-      CALL xc_f03_func_init( xc_func, igcx, sp )
+      CALL xc_f03_gga_fxc( xc_func(3), lengthxc, rho_lbxc(1), sigma(1), v2rho2_x(1), &
-       CALL xc_f03_gga_fxc( xc_func, lengthxc, rho_lbxc(1), sigma(1), v2rho2_x(1), &
                            v2rhosigma_x(1), v2sigma2_x(1) )
-      CALL xc_f03_func_end( xc_func )
     ENDIF
   ENDIF
   !
@@ -132,12 +129,9 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
     ! ... DERIVATIVE FOR CORRELATION
     v2rho2_c = 0._DP ;  v2rhosigma_c = 0._DP ;  v2sigma2_c = 0._DP
     IF (igcc /= 0) THEN 
-      CALL xc_f03_func_init( xc_func, igcc, sp )
+      fkind = xc_f03_func_info_get_kind( xc_info(4) )
-       xc_info2 = xc_f03_func_get_info( xc_func )
+      CALL xc_f03_gga_fxc( xc_func(4), lengthxc, rho_lbxc(1), sigma(1), v2rho2_c(1), &
-       fkind = xc_f03_func_info_get_kind( xc_info2 )
-       CALL xc_f03_gga_fxc( xc_func, lengthxc, rho_lbxc(1), sigma(1), v2rho2_c(1), &
                            v2rhosigma_c(1), v2sigma2_c(1) )
-      CALL xc_f03_func_end( xc_func )
     ENDIF
   ENDIF
   !

XClib/xc_wrapper_d_lda_lsda.f90

@@ -18,6 +18,7 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
 #if defined(__LIBXC)
 #include "xc_version.h"
   USE xc_f03_lib_m
+  USE dft_par_mod,     ONLY: xc_func, xc_info
 #endif
   !
   IMPLICIT NONE
@@ -34,8 +35,6 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
   ! ... local variables
   !
 #if defined(__LIBXC)
-  TYPE(xc_f03_func_t) :: xc_func
-  TYPE(xc_f03_func_info_t) :: xc_info2
   INTEGER :: pol_unpol, fkind_x
   REAL(DP), ALLOCATABLE :: rho_lxc(:)
   REAL(DP), ALLOCATABLE :: dmex_lxc(:), dmcr_lxc(:)
@@ -96,9 +95,7 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
     ! ... DERIVATIVE FOR EXCHANGE
     dmex_lxc(:) = 0.0_DP
     IF (iexch /= 0) THEN    
-       CALL xc_f03_func_init( xc_func, iexch, pol_unpol )
+       CALL xc_f03_lda_fxc( xc_func(1), lengthxc, rho_lxc(1), dmex_lxc(1) )
-        CALL xc_f03_lda_fxc( xc_func, lengthxc, rho_lxc(1), dmex_lxc(1) )
-       CALL xc_f03_func_end( xc_func )
     ENDIF  
   ENDIF
   !
@@ -107,11 +104,8 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
     ! ... DERIVATIVE FOR CORRELATION
     dmcr_lxc(:) = 0.0_DP
     IF (icorr /= 0) THEN
-       CALL xc_f03_func_init( xc_func, icorr, pol_unpol )
+       fkind_x  = xc_f03_func_info_get_kind( xc_info(2) )
-        xc_info2 = xc_f03_func_get_info( xc_func )
+       CALL xc_f03_lda_fxc( xc_func(2), lengthxc, rho_lxc(1), dmcr_lxc(1) )
-        fkind_x  = xc_f03_func_info_get_kind( xc_info2 )
-        CALL xc_f03_lda_fxc( xc_func, lengthxc, rho_lxc(1), dmcr_lxc(1) )
-       CALL xc_f03_func_end( xc_func )
     ENDIF
   ENDIF
   !

XClib/xc_wrapper_gga.f90

@@ -17,6 +17,7 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
 #if defined(__LIBXC)
 #include "xc_version.h"
   USE xc_f03_lib_m
+  USE dft_par_mod,   ONLY: xc_func, xc_info 
 #endif
   !
   USE kind_l,        ONLY: DP
@@ -52,8 +53,6 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
   ! ... local variables
   !
 #if defined(__LIBXC)
-  TYPE(xc_f03_func_t) :: xc_func
-  TYPE(xc_f03_func_info_t) :: xc_info1, xc_info2
   REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:)
   REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
   REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:)
@@ -151,12 +150,9 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
   !
   IF ( is_libxc(4) ) THEN  !lda part of LYP not present in libxc (still so? - check)
     !
-    CALL xc_f03_func_init( xc_func, igcc, pol_unpol )
+    CALL xc_f03_func_set_dens_threshold( xc_func(4), rho_threshold_gga )
-     xc_info2 = xc_f03_func_get_info( xc_func )
+    fkind_x  = xc_f03_func_info_get_kind( xc_info(4) )
-     CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_gga )
+    CALL xc_f03_gga_exc_vxc( xc_func(4), lengthxc, rho_lxc(1), sigma(1), ec_lxc(1), vc_rho(1), vc_sigma(1) )
-     fkind_x  = xc_f03_func_info_get_kind( xc_info2 )
-     CALL xc_f03_gga_exc_vxc( xc_func, lengthxc, rho_lxc(1), sigma(1), ec_lxc(1), vc_rho(1), vc_sigma(1) )
-    CALL xc_f03_func_end( xc_func )
     !
     IF (.NOT. POLARIZED) THEN
       DO k = 1, length
@@ -238,11 +234,8 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
   !
   IF ( is_libxc(3) ) THEN
     !
-    CALL xc_f03_func_init( xc_func, igcx, pol_unpol )
+    CALL xc_f03_func_set_dens_threshold( xc_func(3), rho_threshold_gga )
-     xc_info1 = xc_f03_func_get_info( xc_func )
+    CALL xc_f03_gga_exc_vxc( xc_func(3), lengthxc, rho_lxc(1), sigma(1), ex_lxc(1), vx_rho(1), vx_sigma(1) )
-     CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_gga )
-     CALL xc_f03_gga_exc_vxc( xc_func, lengthxc, rho_lxc(1), sigma(1), ex_lxc(1), vx_rho(1), vx_sigma(1) )
-    CALL xc_f03_func_end( xc_func )
     !
     IF (.NOT. POLARIZED) THEN
       DO k = 1, length

XClib/xc_wrapper_lda_lsda.f90

@@ -14,6 +14,7 @@ SUBROUTINE xc( length, sr_d, sv_d, rho_in, ex_out, ec_out, vx_out, vc_out )
 #if defined(__LIBXC)
 #include "xc_version.h"
   USE xc_f03_lib_m
+  USE dft_par_mod,       ONLY: xc_func, xc_info
 #endif
   !
   USE kind_l,            ONLY: DP
@@ -43,8 +44,6 @@ SUBROUTINE xc( length, sr_d, sv_d, rho_in, ex_out, ec_out, vx_out, vc_out )
   ! ... local variables
   !
 #if defined(__LIBXC)
-  TYPE(xc_f03_func_t) :: xc_func
-  TYPE(xc_f03_func_info_t) :: xc_info1, xc_info2
   INTEGER :: fkind_x
   REAL(DP) :: amag
   REAL(DP), ALLOCATABLE :: rho_lxc(:)
@@ -105,21 +104,15 @@ SUBROUTINE xc( length, sr_d, sv_d, rho_in, ex_out, ec_out, vx_out, vc_out )
   !
   ! ... EXCHANGE
   IF ( is_libxc(1) ) THEN
-     CALL xc_f03_func_init( xc_func, iexch, sv_d )
+     CALL xc_f03_func_set_dens_threshold( xc_func(1), rho_threshold_lda )
-       xc_info1 = xc_f03_func_get_info( xc_func )
+     CALL xc_f03_lda_exc_vxc( xc_func(1), lengthxc, rho_lxc(1), ex_out(1), vx_lxc(1) )
-       CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_lda )
-       CALL xc_f03_lda_exc_vxc( xc_func, lengthxc, rho_lxc(1), ex_out(1), vx_lxc(1) )
-     CALL xc_f03_func_end( xc_func )
   ENDIF
   !
   ! ... CORRELATION
   IF ( is_libxc(2) ) THEN
-     CALL xc_f03_func_init( xc_func, icorr, sv_d )
+     CALL xc_f03_func_set_dens_threshold( xc_func(2), rho_threshold_lda )
-      xc_info2 = xc_f03_func_get_info( xc_func )
+     fkind_x = xc_f03_func_info_get_kind( xc_info(2) )
-      CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_lda )
+     CALL xc_f03_lda_exc_vxc( xc_func(2), lengthxc, rho_lxc(1), ec_out(1), vc_lxc(1) )
-      fkind_x  = xc_f03_func_info_get_kind( xc_info2 )
-      CALL xc_f03_lda_exc_vxc( xc_func, lengthxc, rho_lxc(1), ec_out(1), vc_lxc(1) )
-     CALL xc_f03_func_end( xc_func )
   ENDIF
   !
   IF ( ((.NOT.is_libxc(1)) .OR. (.NOT.is_libxc(2))) &

XClib/xc_wrapper_mgga.f90

@@ -15,6 +15,7 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
 #include "xc_version.h"
   USE xc_f03_lib_m
   USE dft_mod,       ONLY: get_libxc_flags_exc
+  USE dft_par_mod,   ONLY: xc_func, xc_info
 #endif 
   !
   USE kind_l,        ONLY: DP
@@ -60,9 +61,6 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
   REAL(DP), ALLOCATABLE :: grho2(:,:)
   !
 #if defined(__LIBXC)
-  TYPE(xc_f03_func_t) :: xc_func
-  TYPE(xc_f03_func_info_t) :: xc_info1, xc_info2
-  !
   REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:), tau_lxc(:)
   REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
   REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:), vx_tau(:)
@@ -149,18 +147,15 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
   ! META EXCHANGE
   !
   IF ( is_libxc(5) ) THEN
-     CALL xc_f03_func_init( xc_func, imeta, pol_unpol )
+    CALL xc_f03_func_set_dens_threshold( xc_func(5), rho_threshold_mgga )
-     xc_info1 = xc_f03_func_get_info( xc_func )
+    CALL get_libxc_flags_exc( xc_info(5), eflag )
-     CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_mgga )
-     CALL get_libxc_flags_exc( xc_info1, eflag )
     IF (eflag==1) THEN
-       CALL xc_f03_mgga_exc_vxc( xc_func, lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
+      CALL xc_f03_mgga_exc_vxc( xc_func(5), lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
                                 ex_lxc(1), vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
     ELSE
-       CALL xc_f03_mgga_vxc( xc_func, lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
+      CALL xc_f03_mgga_vxc( xc_func(5), lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
                             vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
     ENDIF
-    CALL xc_f03_func_end( xc_func )
     !
     IF (.NOT. POLARIZED) THEN
       DO k = 1, length
@@ -197,12 +192,9 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
   !
   IF ( is_libxc(6) ) THEN
     !
-    CALL xc_f03_func_init( xc_func, imetac, pol_unpol )
+    CALL xc_f03_func_set_dens_threshold( xc_func(6), rho_threshold_mgga )
-    xc_info1 = xc_f03_func_get_info( xc_func )
+    CALL xc_f03_mgga_exc_vxc( xc_func(6), lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
-    CALL xc_f03_func_set_dens_threshold( xc_func, rho_threshold_mgga )
-    CALL xc_f03_mgga_exc_vxc( xc_func, lengthxc, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
                               ec_lxc(1), vc_rho(1), vc_sigma(1), vlapl_rho(1), vc_tau(1) )
-    CALL xc_f03_func_end( xc_func )
     !
     IF (.NOT. POLARIZED) THEN
        DO k = 1, length

XClib/xclib_test.f90

@@ -29,10 +29,12 @@ PROGRAM xclib_test
   USE constants_l, ONLY: pi
   USE xc_lib,      ONLY: xclib_set_dft_from_name, xclib_set_exx_fraction, &
                          xclib_get_ID, xclib_reset_dft, xc_gcx,           &
-                         xclib_dft_is_libxc
+                         xclib_dft_is_libxc, xclib_init_libxc,            &
+                         xclib_finalize_libxc
   USE xclib_parallel_include
 #if defined(__LIBXC)
-  USE xc_f90_lib_m
+  USE xc_f03_lib_m
+  USE dft_par_mod, ONLY: xc_func, xc_info
 #endif
   !
   IMPLICIT NONE
@@ -44,8 +46,6 @@ PROGRAM xclib_test
 #endif
   !
 #if defined(__LIBXC)
-  TYPE(xc_f90_func_t) :: xc_func
-  TYPE(xc_f90_func_info_t) :: xc_info
   CHARACTER(LEN=120) :: lxc_kind, lxc_family
   INTEGER :: n_ext, id(6)
 #endif
@@ -314,6 +314,9 @@ PROGRAM xclib_test
     !
     !
 #if defined(__LIBXC)
+    !
+    IF (xclib_dft_is_libxc('ANY')) CALL xclib_init_libxc( 1 )
+    !
     WRITE(stdout,*) CHAR(10)//"LIBXC functional infos:"
     !
     id(1) = iexch1 ; id(2) = icorr1
@@ -324,11 +327,7 @@ PROGRAM xclib_test
       IF (is_libxc(i)) THEN
         WRITE(stdout,*) CHAR(10)//"Functional with ID:", id(i)
         !
-        CALL xc_f90_func_init( xc_func, id(i), 1 )
+        SELECT CASE( xc_f03_func_info_get_kind(xc_info(i)) )
-        !
-        xc_info = xc_f90_func_get_info(xc_func)
-        !
-        SELECT CASE( xc_f90_func_info_get_kind(xc_info) )
         CASE( XC_EXCHANGE )
           WRITE(lxc_kind, '(a)') 'Exchange functional'
         CASE( XC_CORRELATION )
@@ -342,7 +341,7 @@ PROGRAM xclib_test
           WRITE(lxc_kind, '(a)') 'Unknown kind'
         END SELECT
         !
-        SELECT CASE( xc_f90_func_info_get_family(xc_info) )
+        SELECT CASE( xc_f03_func_info_get_family(xc_info(i)) )
         CASE( XC_FAMILY_LDA )
           WRITE(lxc_family,'(a)') "LDA"
         CASE( XC_FAMILY_GGA )
@@ -359,31 +358,32 @@ PROGRAM xclib_test
         !
         WRITE(*,'("The functional ''", a, "'' is an ", a, ", it belongs to &
                &the ''", a, "'' family and is defined in the reference(s): &
-               &")') TRIM(xc_f90_func_info_get_name(xc_info)), TRIM(lxc_kind)&
+               &")') TRIM(xc_f03_func_info_get_name(xc_info(i))), TRIM(lxc_kind)&
                ,TRIM(lxc_family)
         ii = 0
         DO WHILE( ii >= 0 )
-         WRITE(*,'(a,i1,2a)') '[',ii+1,'] ',TRIM(xc_f90_func_reference_get_ref( &
+         WRITE(*,'(a,i1,2a)') '[',ii+1,'] ',TRIM(xc_f03_func_reference_get_ref( &
-                                   xc_f90_func_info_get_references(xc_info, ii)))
+                                   xc_f03_func_info_get_references(xc_info(i), ii)))
         ENDDO
         !
         WRITE(stdout,*)
-        n_ext = xc_f90_func_info_get_n_ext_params( xc_info )
+        n_ext = xc_f03_func_info_get_n_ext_params( xc_info(i) )
         WRITE(stdout,*) 'Number of external parameters: ', n_ext
         !
         IF ( n_ext/=0 ) THEN
           DO ii = 0, n_ext-1
             WRITE(stdout,*) &
-              TRIM(xc_f90_func_info_get_ext_params_description(xc_info, ii))
+              TRIM(xc_f03_func_info_get_ext_params_description(xc_info(i), ii))
             WRITE(stdout,*) 'Default value: ', &
-                   xc_f90_func_info_get_ext_params_default_value(xc_info, ii)
+                   xc_f03_func_info_get_ext_params_default_value(xc_info(i), ii)
           ENDDO
         ENDIF
         !
-        CALL xc_f90_func_end( xc_func )
-        !
       ENDIF
     ENDDO
+    !
+    IF (xclib_dft_is_libxc('ANY')) CALL xclib_finalize_libxc()
+    !
 #endif    
     !
     121 FORMAT('Exch: ',I3,' is libxc: ',L1,';  Corr: ',I3,' is libxc: ',L1 )
@@ -475,6 +475,8 @@ PROGRAM xclib_test
   np = 1
   IF (ns==2) np = 3
   !
+  IF (xclib_dft_is_libxc('ANY')) CALL xclib_init_libxc( ns )
+  !
   !==========================================================================
   ! ALLOCATIONS OF XC I/O ARRAYS
   !==========================================================================
@@ -767,6 +769,7 @@ PROGRAM xclib_test
     IF (test == 'dft-comparison') THEN
       CALL xclib_reset_dft()
       CALL xclib_set_dft_from_name( dft2 )
+      IF (xclib_dft_is_libxc('ANY')) CALL xclib_init_libxc( ns )
     ENDIF
     IF (test(1:4)=='gen-') dft2 = dft1
     !
@@ -1323,6 +1326,8 @@ PROGRAM xclib_test
   ! FINALIZE
   !==========================================================================
   !
+  IF (xclib_dft_is_libxc('ANY')) CALL xclib_finalize_libxc()
+  !
   DEALLOCATE( rho, rho_tz )
   !
   IF ( GGA .OR. MGGA ) DEALLOCATE( grho )