Merge lp:~alesster/siesta/trunk_new_xc into lp:siesta

[Aleksandr V. Terentjev]

There were realized 4 new functionals:
SG4, GGA functional, keyword for XC.authors in fdf: SG4, sg4
Ref. L. A. Constantin et al, Phys. Rev. B 93, 045126 (2016)

SG4+VV10m, VDW functional, keyword for XC.authors in fdf: SG4VV, sg4vv
Ref. A. V. Terentjev et al, Computation 6, 7 (2018)

PBE+VV10, VDW functional, keyword for XC.authors in fdf: PBEVV, pbevv, PBEvv, pbeVV
Ref. H. Peng and J. P. Perdew, Phys. Rev. B 95, 0811105 (2017)

PBEsol+VV10s, VDW functional, keyword for XC.authors in fdf: PVESVVS, pbesvvs, PBESVVs, PBESvvs
Ref. A. V. Terentjev et al, Phys. Rev. B 98, 214108 (2018)

IMPORTANT NOTES:
During the testing two test examples were crashed with version 781:
h2o_filteret_basis and si001-diags

With new version 782 all tests were finished successfully exepting h2o_filteret_basis and si001-diags.

Check has shown the differences bwtween tested and refrenced output files. But it also gives the differences for the previous version 781.

[Alberto Garcia (albertog)]

Sorry about the delay in handling your merge request for new vdW options in SiestaXC.
The delay was motivated by two issues:

1. The Gitlab transition.
2. The development of a stand-alone library, ligGridXC, based on the code in SiestaXC.

LibGridXC is being developed at

    https://gitlab.com/siesta-project/libraries/libgridxc

and is "the future", in the sense that it will be required as a external library for Siesta very soon, once we merge in the next few weeks the PSML branch.

LibGridXC has the same functionality of SiestaXC, with some streamlining, and, in addition, offers an interface to LibXC. This interface is optional now, but very likely it will be made compulsory in the future.

As I understand it, your merge request:

* Adds a few more vdW flavors.
* Implements a new semilocal SG4 functional that is needed by some of the new vdW flavors.

I have no problem with the content of your merge request. I think it should definitely go into libGridXC, but I am not sure whether to integrate it also into the current code base (based on SiestaXC) which is "going away" soon.

In any case, I would request that you make new merge requests within Gitlab:

1. One for libgridXC.
2. (If you think it is worth it) Another one for the "master" branch of Siesta:

             https://gitlab.com/siesta-project/siesta

Note that, as part of the Gitlab transition, the Launchpad branch that you proposed for merging is now available in git form in:

             https://gitlab.com/siesta-project/siesta-legacy/alesster_trunk_new_xc

As a further comment, note that the SG4 functional is already provided by libXC. However, the current way in which the semilocal pieces of the vdW functionals is handled in ligGridXC needs to be checked for compatibility with libXC. It is not as easy as passing the appropriate LIBXC-XXXX string to ggaxc, since the functional object needs to be initialized in advance. So I would keep the explicit implementation of SG4 for now.

Unmerged revisions

782. By Aleksandr V. Terentjev on 2019-09-06

-------------- This line and the following will be ignored --------------

modified:
  Docs/siesta.tex
  Src/SiestaXC/ggaxc.f
  Src/SiestaXC/vdwxc.F90
  Src/SiestaXC/vv_vdwxc.F90
  Src/atom.F
  Tests/Makefile
  version.info
unknown:
  Docs/branch-changes/NOTES.trunk_new_xc
  Tests/ag_sg4/
  Tests/graphite_c6_sg4/
  Tests/graphite_vdw_pbesvvs/
  Tests/graphite_vdw_pbevv/
  Tests/graphite_vdw_sg4vv/
  Tests/Reference/ag_sg4.out
  Tests/Reference/graphite_c6_sg4.out
  Tests/Reference/graphite_vdw_pbesvvs.out
  Tests/Reference/graphite_vdw_pbevv.out
  Tests/Reference/graphite_vdw_sg4vv.out

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2019-08-22 09:05:46 +0000
+++ Docs/siesta.tex	2019-09-06 14:47:07 +0000
@@ -3723,6 +3723,10 @@
     implemented in \siesta)
 
     \item%
+    \fdf*{SG4}: GGA functional based on semiclassical atom theory of Lucian A. Constantin, Aleksandrs Terentjevs, Fabio Della Sala, Pietro Cortona and Eduardo Fabiano, 
+    Phys. Rev. B \textbf{93}, 045126 (2016).
+
+    \item%
     \fdf*{DRSLL} (equivalent to \fdf*{DF1}): \xcidx{vdW-DF1}
     \xcidx{DRSLL}%
     van der Waals \xcidx{vdW} density functional (vdW-DF) \xcidx{vdW-DF}
@@ -3757,6 +3761,19 @@
     \fdf*{VV}: \xcidx{VV}%
     vdW-DF functional of O. A. Vydrov and T. Van Voorhis, 
     J. Chem. Phys. \textbf{133}, 244103 (2010)
+
+    \item%
+    \fdf*{PBEVV}: vdW-DF functional (PBE+VV10L) with PBE exchange and correlation functionals and with non-local correlation term of O. A. Vydrov and T. Van Voorhis (VV).    PBE+VV10L is the modification of PBE+rVV10L  for \siesta\  of H. Peng and J. P. Perdew, Phys. Rev. B  \textbf{95}, 0811105 (2017).
+    
+    \item%
+    \fdf*{SG4VV}: vdW-DF functional (SG4+VV10m) with SG4 exchange and correlation functionals and with non-local correlation term of O. A. Vydrov and T. Van Voorhis (VV).
+    SG4+VV10m is the modification of SG4+rVV10m functional for \siesta\ of Aleksandr V. Terentjev, Pietro Cortona, Lucian A. Constantin, José M. Pitarke 
+    Fabio Della Sala and Eduardo Fabiano, Computation \textbf{6}, 7 (2018).
+    
+    \item%
+    \fdf*{PBESVVS}: vdW-DF functional (PBEsol+VV10s) with PBEsol exchange and correlation functionals and with non-local correlation term of O. A. Vydrov and 
+    T. Van Voorhis (VV). PBEsol+VV10s is the modification of PBEsol+rVV10s functional for \siesta\ of Aleksandr V. Terentjev, Lucian A. Constantin and J. M. Pitarke, 
+    Phys. Rev. B \textbf{98}, 214108 (2018).
     
   \end{itemize}
 
 
=== modified file 'Src/SiestaXC/ggaxc.f'
--- Src/SiestaXC/ggaxc.f	2018-04-19 10:08:47 +0000
+++ Src/SiestaXC/ggaxc.f	2019-09-06 14:47:07 +0000
@@ -24,8 +24,9 @@
 ! b88x,        ! Becke, PRA 38, 3098 (1988) (exchange only)
 ! b88kbmx,     ! Klimes et al, JPCM 22, 022201 (2009) (exchange only)
 ! c09x,        ! Cooper, PRB 81, 161104 (2010) (exchange only)
-! bhx          ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exchange only)
-!
+! bhx,         ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exchange only)     
+! sg4xc        ! GGA Constantin et al, PRB, 93, 045126 (2016)
+!     
 !   PUBLIC parameters, types, and variables available from this module:
 ! none
 !
@@ -76,8 +77,9 @@
      .  b88x,         ! Becke, PRA 38, 3098 (1988) (exchange only)
      .  b88kbmx,      ! Klimes et al, JPCM 22, 022201 (2009) (exchange only)
      .  c09x,         ! Cooper, PRB 81, 161104 (2010) (exchange only)
-     .  bhx           ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exch. only)
-
+     .  bhx,          ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exch. only)
+     .  sg4xc         ! GGA Constantin et al, PRB, 93, 045126 (2016)
+      
       PRIVATE  ! Nothing is declared public beyond this point
 
       CONTAINS
@@ -255,6 +257,10 @@
         CALL PBESOLXC( IREL, NS, DD, GDD,
      .                 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
 
+      ELSEIF (AUTHOR.EQ.'SG4'.OR.AUTHOR.EQ.'sg4') THEN
+        CALL SG4XC( IREL, NS, DD, GDD,
+     .        EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )      
+        
       ELSEIF (AUTHOR.EQ.'AM05' .OR. AUTHOR.EQ.'am05') THEN
         CALL AM05XC( IREL, NS, DD, GDD,
      .               EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD )
@@ -2756,4 +2762,263 @@
 
       END SUBROUTINE C09X
 
-      END MODULE m_ggaxc
+!      END MODULE m_ggaxc
+
+c     
+c the SG4 subroutine
+c 
+
+      SUBROUTINE SG4XC( IREL, nspin, Dens, GDens,
+     .                  EX, EC, DEXDD, DECDD, DEXDGD, DECDGD )
+
+C *********************************************************************
+C Implements GGA based on Semiclassical atom theory
+C Ref: L.A.Constantin et al, PRB, 93, 045126 (2016)  
+C ******** INPUT ******************************************************
+C INTEGER IREL           : Relativistic-exchange switch (0=No, 1=Yes)
+C INTEGER nspin          : Number of spin polarizations (1 or 2)
+C REAL*8  Dens(nspin)    : Total electron density (if nspin=1) or
+C                           spin electron density (if nspin=2)
+C REAL*8  GDens(3,nspin) : Total or spin density gradient
+C ******** OUTPUT *****************************************************
+C REAL*8  EX             : Exchange energy density
+C REAL*8  EC             : Correlation energy density
+C REAL*8  DEXDD(nspin)   : Partial derivative
+C                           d(DensTot*Ex)/dDens(ispin),
+C                           where DensTot = Sum_ispin( Dens(ispin) )
+C                          For a constant density, this is the
+C                          exchange potential
+C REAL*8  DECDD(nspin)   : Partial derivative
+C                           d(DensTot*Ec)/dDens(ispin),
+C                           where DensTot = Sum_ispin( Dens(ispin) )
+C                          For a constant density, this is the
+C                          correlation potential
+C REAL*8  DEXDGD(3,nspin): Partial derivative
+C                           d(DensTot*Ex)/d(GradDens(i,ispin))
+C REAL*8  DECDGD(3,nspin): Partial derivative
+C                           d(DensTot*Ec)/d(GradDens(i,ispin))
+C ********* UNITS ****************************************************
+C Lengths in Bohr
+C Densities in electrons per Bohr**3
+C Energies in Hartrees
+C Gradient vectors in cartesian coordinates
+C ********* ROUTINES CALLED ******************************************
+C EXCHNG, PW92C
+C ********************************************************************
+
+      use precision, only : dp
+
+      implicit          none
+      INTEGER           IREL, nspin
+      real(dp)          Dens(nspin), DECDD(nspin), DECDGD(3,nspin),
+     .                  DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin) 
+      
+C Internal variables
+      INTEGER
+     .  IS, IX
+
+      real(dp)
+     .  A, BETA, D(2), DADD, DECUDD, DENMIN,
+     .  DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD,
+     .  DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2),
+     .  DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD,
+     .  DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2),
+     .  EC, ECUNIF, EX, EXUNIF,
+     .  F, F1, F2, F3, F4, FC, FX, FOUTHD,
+     .  GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3),
+     .  H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S,
+     .  T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA, 
+c Lucian: new variables
+     .  beta0, sigma, alpha, amu1, amu2, ak1, ak2, term1, dbetadd,
+     .  dterm1dd, dbetadgd, dadgd, dterm1dgd, xx1, xx2, dxx1ds, dxx2ds,
+     .  dfds
+
+
+C Lower bounds of density and its gradient to avoid divisions by zero
+      PARAMETER ( DENMIN = 1.D-12 )
+      PARAMETER ( GDMIN  = 1.D-12 )
+
+C Fix some numerical parameters
+      PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0,
+     .            THD=1.D0/3.D0, THRHLF=1.5D0,
+     .            TWO=2.D0, TWOTHD=2.D0/3.D0 )
+
+C Fix some more numerical constants
+c Lucian: constants for SG4 correlation and exchange
+      PI = 4 * ATAN(1.D0)
+      GAMMA = (1 - LOG(TWO)) / PI**2
+      BETA0=0.07903051884257d0
+      sigma=0.07d0
+      alpha=0.8d0
+c
+      amu1=0.042d0
+      amu2=0.218d0     
+      ak1=0.56028717948717954d0     
+      ak2=0.24371282051282048d0
+C Translate density and its gradient to new variables
+      IF (nspin .EQ. 1) THEN
+         if(Dens(1).le.1.0d-12) Dens(1)=1.0d-12
+        D(1) = HALF * Dens(1)
+        D(2) = D(1)
+        DT = MAX( DENMIN, Dens(1) )
+        DO 10 IX = 1,3
+          GD(IX,1) = HALF * GDens(IX,1)
+          GD(IX,2) = GD(IX,1)
+          GDT(IX) = GDens(IX,1)
+   10   CONTINUE
+      ELSE
+         if(Dens(1).le.1.0d-12) Dens(1)=1.0d-12
+         if(Dens(2).le.1.0d-12) Dens(2)=1.0d-12
+        D(1) = Dens(1)
+        D(2) = Dens(2)
+        DT = MAX( DENMIN, Dens(1)+Dens(2) )
+        DO 20 IX = 1,3
+          GD(IX,1) = GDens(IX,1)
+          GD(IX,2) = GDens(IX,2)
+          GDT(IX) = GDens(IX,1) + GDens(IX,2)
+   20   CONTINUE
+      ENDIF
+      GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 )
+      GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 )
+      GDMT   = SQRT( GDT(1)**2  + GDT(2)**2  + GDT(3)**2  )
+      GDMT = MAX( GDMIN, GDMT )
+
+C Find local correlation energy and potential
+      CALL PW92C( 2, D, ECUNIF, VCUNIF )
+
+C Find total correlation energy
+      RS = ( 3 / (4*PI*DT) )**THD
+      KF = (3 * PI**2 * DT)**THD
+      KS = SQRT( 4 * KF / PI )
+      ZETA = ( D(1) - D(2) ) / DT
+      ZETA = MAX( -1.D0+DENMIN, ZETA )
+      ZETA = MIN(  1.D0-DENMIN, ZETA )
+      PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD )
+      T = GDMT / (2 * PHI * KS * DT)
+c Define beta
+      if(t.ge.100.d0)t=100.d0      
+      beta=beta0+sigma*t*(1.d0-dexp(-rs**2))
+c
+      F1 = ECUNIF / GAMMA / PHI**3
+      F2 = EXP(-F1)
+      A = BETA / GAMMA / (F2-1)
+      F3 = T**2 + A * T**4
+      F4 = BETA/GAMMA * F3 / (1 + A*F3)
+c Define phi**(alpha*t**3) term
+      term1 = phi ** ( alpha * t **3 )
+      H = term1 * GAMMA * PHI**3 * LOG( 1 + F4 )
+c
+      FC = ECUNIF + H
+
+C Find correlation energy derivatives
+      DRSDD = - (THD * RS / DT)
+      DKFDD =   THD * KF / DT
+      DKSDD = HALF * KS * DKFDD / KF
+      DZDD(1) =   1 / DT - ZETA / DT
+      DZDD(2) = - (1 / DT) - ZETA / DT
+      DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD )
+      DO 40 IS = 1,2
+        DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT
+        DPDD = DPDZ * DZDD(IS)
+        DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT )
+c Derivative of beta
+        dbetadd = sigma * DTDD * ( 1.d0 - dexp ( - rs ** 2 ) ) +
+     .            sigma * t * 2.d0 * rs * dexp ( - rs ** 2 ) * DRSDD
+c
+        DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI )
+        DF2DD = (- F2) * DF1DD
+c Derivative of A
+        DADD = (- A) * DF2DD / (F2-1) + dbetadd / GAMMA / (F2-1)
+c
+        DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4
+        DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) )+
+     .          dbetadd/GAMMA * F3 / (1 + A*F3)
+
+c Derivative of term1
+        dterm1dd = term1 * ( alpha * 3.d0 * t * t * dtdd * dlog(phi) +
+     .             alpha * t * t * t / phi * dpdd )
+c
+c Derivative of H 
+        DHDD= 3 * phi**2 * term1 * GAMMA * LOG( 1 + F4 ) * DPDD +
+     .  dterm1dd * GAMMA * PHI**3 * LOG( 1 + F4 ) +
+     .  term1 * GAMMA * PHI**3 * DF4DD / (1+F4)
+c        
+        DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD
+        DO 30 IX = 1,3
+          DTDGD = (T / GDMT) * GDT(IX) / GDMT
+c Derivative of beta with respect to gradient
+          dbetadgd = sigma * DTDGD * ( 1.d0 - dexp ( - rs**2 ) )
+c Derivative of A with respect to gradient
+          dadgd = dbetadgd / GAMMA / (F2-1)
+c
+          DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 ) +
+     .             dadgd * T**4
+c   
+          DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) )+
+     .             dbetadgd /GAMMA * F3 / (1 + A*F3) -
+     .             BETA /GAMMA * F3 / (1 + A*F3)**2 * dadgd *f3
+c   
+          dterm1dgd = term1 * alpha * 3.d0 * t * t * dtdgd * dlog(phi)          
+c 
+          DHDGD = term1 * GAMMA * PHI**3 * DF4DGD / (1+F4) +
+     .            dterm1dgd * GAMMA * PHI**3 * LOG( 1 + F4 )
+c       
+          DFCDGD(IX,IS) = DT * DHDGD
+   30   CONTINUE
+   40 CONTINUE
+
+C Find exchange energy and potential
+      FX = 0
+      DO 60 IS = 1,2
+        DS(IS)   = MAX( DENMIN, 2 * D(IS) )
+        GDMS = MAX( GDMIN, 2 * GDM(IS) )
+        KFS = (3 * PI**2 * DS(IS))**THD
+        S = GDMS / (2 * KFS * DS(IS))
+c Exchange enhancement factor 
+        xx1=amu1*s*s/ak1
+        xx2=amu2*s*s/ak2
+        F=1+ak1+ak2-ak1/(1.d0+xx1+xx1**2+xx1**3+xx1**4)-
+     .    ak2/(1.d0+xx2)
+c
+c       Note nspin=1 in call to exchng...
+c
+        CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) )
+        FX = FX + DS(IS) * EXUNIF * F
+
+        DKFDD = THD * KFS / DS(IS)
+        DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) )
+c Derivative of exchange enhancement factor
+        dxx1ds=amu1*2.d0*s/ak1
+        dxx2ds=amu2*2.d0*s/ak2
+        dfds=ak1/(1.d0+xx1+xx1**2+xx1**3+xx1**4)**2 * 
+     .       dxx1ds*(1.d0+2.d0*xx1+3.d0*xx1**2+4.d0*xx1**3) +
+     .       ak2/(1.d0+xx2)**2 * dxx2ds
+        dfdd=dfds*dsdd
+        DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD
+
+        DO 50 IX = 1,3
+          GDS = 2 * GD(IX,IS)
+          DSDGD = (S / GDMS) * GDS / GDMS
+c  
+          DFDGD = dfds * DSDGD
+c
+          DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD
+   50   CONTINUE
+   60 CONTINUE
+      FX = HALF * FX / DT
+
+C Set output arguments
+      EX = FX
+      EC = FC
+      DO 90 IS = 1,nspin
+        DEXDD(IS) = DFXDD(IS)
+        DECDD(IS) = DFCDD(IS)
+        DO 80 IX = 1,3
+          DEXDGD(IX,IS) = DFXDGD(IX,IS)
+          DECDGD(IX,IS) = DFCDGD(IX,IS)
+   80   CONTINUE
+   90 CONTINUE
+
+      END SUBROUTINE SG4XC
+      
+      END MODULE m_ggaxc      
 
=== modified file 'Src/SiestaXC/vdwxc.F90'
--- Src/SiestaXC/vdwxc.F90	2015-05-08 09:34:16 +0000
+++ Src/SiestaXC/vdwxc.F90	2019-09-06 14:47:07 +0000
@@ -257,12 +257,33 @@
   use m_vv_vdwxc,  only: vv_vdw_get_kmesh  ! Size and values of (kf,kg) mesh
   use m_vv_vdwxc,  only: vv_vdw_phi        ! Interpolates rho1*rho2*phi(k1,k2,k)
   use m_vv_vdwxc,  only: vv_vdw_set_kcut   ! Sets the planewave cutoff kc
-
+!! New vdw 
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_beta       ! Parameter beta of VV2010 functional
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_theta      ! Func. theta of VV2010 functional
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_get_kmesh  ! Size and values of (kf,kg) mesh
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_phi        ! Interpolates rho1*rho2*phi(k1,k2,k)
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_set_kcut   ! Sets the planewave cutoff kc
+
+  use m_pbevv_vdwxc,  only: pbevv_vdw_beta       ! Parameter beta of VV2010 functional
+  use m_pbevv_vdwxc,  only: pbevv_vdw_theta      ! Func. theta of VV2010 functional
+  use m_pbevv_vdwxc,  only: pbevv_vdw_get_kmesh  ! Size and values of (kf,kg) mesh
+  use m_pbevv_vdwxc,  only: pbevv_vdw_phi        ! Interpolates rho1*rho2*phi(k1,k2,k)
+  use m_pbevv_vdwxc,  only: pbevv_vdw_set_kcut   ! Sets the planewave cutoff kc
+
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_beta       ! Parameter beta of VV2010 functional
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_theta      ! Func. theta of VV2010 functional
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_get_kmesh  ! Size and values of (kf,kg) mesh
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_phi        ! Interpolates rho1*rho2*phi(k1,k2,k)
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_set_kcut   ! Sets the planewave cutoff kc
 ! Used module parameters
   use precision,   only: dp                ! Real double precision type
 
 #ifdef DEBUG_XC
   use m_vv_vdwxc,  only: vv_vdw_phiofr     ! Interpolates rho1*rho2*phi(k1,k2,r)
+!! New use of vdw_phiofr
+  use m_sg4vv_vdwxc,  only: sg4vv_vdw_phiofr
+  use m_pbevv_vdwxc,  only: pbevv_vdw_phiofr
+  use m_pbesvvs_vdwxc,  only: pbesvvs_vdw_phiofr
   use debugXC,     only: udebug     ! File unit for debug output
 !  use plot_module, only: plot
 #endif /* DEBUG_XC */
@@ -348,7 +369,7 @@
   real(dp),parameter:: ytol = 1.e-15_dp     ! Tol. for saturated q
 
   ! Private module variables and arrays
-  character(len=5),save:: vdw_author='DRSLL' ! Functional 'flavour' name
+  character(len=7),save:: vdw_author='PBESVVS' ! Functional 'flavour' name
   real(dp),save:: dmesh(nd)                ! Mesh points for phi(d1,d2) table
   real(dp),save:: qmesh(mq)                ! Mesh points for phi(q1,q2,r)
   real(dp),save:: phi_table(0:3,0:3,nd,nd) ! Coefs. for bicubic interpolation
@@ -748,8 +769,18 @@
   ! Trap VV-version exception
 #ifdef DEBUG_XC
   if (vdw_author=='VV') then
-    call vv_vdw_phiofr( r, phi )
-    return
+     call vv_vdw_phiofr( r, phi )
+    return
+!! New vdw_phiofr
+  elseif (vdw_author=='SG4VV') then
+     call sg4vv_vdw_phiofr( r, phi )
+    return
+  elseif (vdw_author=='PBEVV') then
+     call pbevv_vdw_phiofr( r, phi )
+    return
+  elseif (vdw_author=='PBESVVS') then
+     call pbesvvs_vdw_phiofr( r, phi )
+    return 
   end if
 #endif /* DEBUG_XC */
 
@@ -1362,6 +1393,22 @@
     dedrho = eps
     dedgrho = 0
     return
+!! New vdw_beta
+  elseif (vdw_author=='SG4VV') then
+    eps = sg4vv_vdw_beta()
+    dedrho = eps
+    dedgrho = 0
+    return
+  elseif (vdw_author=='PBEVV') then
+    eps = pbevv_vdw_beta()
+    dedrho = eps
+    dedgrho = 0
+    return
+  elseif (vdw_author=='PBESVVS') then
+    eps = pbesvvs_vdw_beta()
+    dedrho = eps
+    dedgrho = 0
+    return 
   end if
 
   if (.not.initialized) then
@@ -1424,11 +1471,30 @@
 
   ! Trap VV-version exception
   if (vdw_author=='VV') then
-    call vv_vdw_get_kmesh( nkf, nkg )
+     call vv_vdw_get_kmesh( nkf, nkg )
     n = nkf*nkg
     if (present(q)) &
       call die('vdw_get_qmesh: ERROR q-mesh not available for author=VV')
     return
+!! New vdw_get_kmesh
+  elseif (vdw_author=='SG4VV') then
+     call sg4vv_vdw_get_kmesh( nkf, nkg )
+    n = nkf*nkg
+    if (present(q)) &
+      call die('vdw_get_qmesh: ERROR q-mesh not available for author=SG4VV')
+    return
+  elseif (vdw_author=='PBEVV') then
+     call pbevv_vdw_get_kmesh( nkf, nkg )
+    n = nkf*nkg
+    if (present(q)) &
+      call die('vdw_get_qmesh: ERROR q-mesh not available for author=PBEVV')
+    return
+  elseif (vdw_author=='PBESVVS') then
+     call pbesvvs_vdw_get_kmesh( nkf, nkg )
+    n = nkf*nkg
+    if (present(q)) &
+      call die('vdw_get_qmesh: ERROR q-mesh not available for author=PBESVVS')
+    return
   end if
 
   if (.not.qmesh_set) call set_qmesh()
@@ -1514,6 +1580,33 @@
                 dEXdD, dEdDaux, dEXdGD, dEdGDaux )
     call GGAxc( 'PBE', iRel, nSpin, D, GD, epsAux, epsC, &
                 dEdDaux, dECdD, dEdGDaux, dECdGD )
+!! New vdw functionals
+  else if (vdw_author=='SG4VV' .or. vdw_author=='sg4vv') then
+    ! A. V. Terentjev et al, Computation 6, 7 (2018)
+    ! SG4 GGA for both exchange and correlation
+    call GGAxc( 'SG4', iRel, nSpin, D, GD, epsX, epsAux, &
+                dEXdD, dEdDaux, dEXdGD, dEdGDaux )
+    call GGAxc( 'SG4', iRel, nSpin, D, GD, epsAux, epsC, &
+                dEdDaux, dECdD, dEdGDaux, dECdGD )
+
+ else if (vdw_author=='PBEVV' .or. vdw_author=='pbevv' .or. &
+          vdw_author=='PBEvv' .or. vdw_author=='pbeVV' ) then
+    ! H. Peng and J. P. Perdew, Phys. Rev. B  95, 0811105 (2017)
+    ! PBE GGA for both exchange and correlation
+    call GGAxc( 'PBE', iRel, nSpin, D, GD, epsX, epsAux, &
+                dEXdD, dEdDaux, dEXdGD, dEdGDaux )
+    call GGAxc( 'PBE', iRel, nSpin, D, GD, epsAux, epsC, &
+                dEdDaux, dECdD, dEdGDaux, dECdGD )
+
+ else if (vdw_author=='PBESVVS' .or. vdw_author=='pbesvvs' .or. &
+          vdw_author=='PBESVVs' .or. vdw_author=='PBESvvs' ) then
+    ! A. V. Terentjev et al, Phys. Rev. B 98, 214108 (2018)
+    ! PBESOL GGA for both exchange and correlation
+    call GGAxc( 'PBESOL', iRel, nSpin, D, GD, epsX, epsAux, &
+                dEXdD, dEdDaux, dEXdGD, dEdGDaux )
+    call GGAxc( 'PBESOL', iRel, nSpin, D, GD, epsAux, epsC, &
+                dEdDaux, dECdD, dEdGDaux, dECdGD )
+
   else
     stop 'vdw_exchng ERROR: unknown author'
   end if
@@ -1542,6 +1635,16 @@
   if (vdw_author=='VV') then
     call vv_vdw_phi( k, phi, dphidk )
     return
+!! New vdw_phi
+  elseif (vdw_author=='SG4VV') then
+    call sg4vv_vdw_phi( k, phi, dphidk )
+    return
+  elseif (vdw_author=='PBEVV') then
+    call pbevv_vdw_phi( k, phi, dphidk )
+    return
+  elseif (vdw_author=='PBESVVS') then
+    call pbesvvs_vdw_phi( k, phi, dphidk )
+    return  
   end if
 
   if (.not.kcut_set) stop 'vdw_phi: ERROR: kcut must be previously set'
@@ -1605,6 +1708,13 @@
     zab = -0.8491_dp
   else if (author=='VV') then
     ! Vydrov and Van Voorhis, JCP 133, 244103 (2010)
+!! New vdw functionals
+  else if (author=='SG4VV') then
+     ! A. V. Terentjev et al, Computation 6, 7 (2018)
+  else if (author=='PBEVV') then
+     ! H. Peng and J. P. Perdew, Phys. Rev. B  95, 0811105 (2017)
+  else if (author=='PBESVVS') then
+     ! A. V. Terentjev et al, Phys. Rev. B  98, 214108 (2018)
   else
     stop 'vdw_set_author: ERROR: author not known'
   end if
@@ -1630,8 +1740,18 @@
 
   ! Trap VV-version exception
   if (vdw_author=='VV') then
-    call vv_vdw_set_kcut( kc )
+     call vv_vdw_set_kcut( kc )
     return
+!! New vdw_set_kcut
+  elseif (vdw_author=='SG4VV') then
+     call sg4vv_vdw_set_kcut( kc )
+     return
+  elseif (vdw_author=='PBEVV') then
+     call pbevv_vdw_set_kcut( kc )
+     return
+  elseif (vdw_author=='PBESVVS') then
+     call pbesvvs_vdw_set_kcut( kc )
+     return  
   end if
 
   if (kc == kcut) return   ! Work alredy done
@@ -1780,6 +1900,16 @@
   if (vdw_author=='VV') then
     call vv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
     return
+!! New vdw_theta
+  elseif (vdw_author=='SG4VV') then
+    call sg4vv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+    return
+  elseif (vdw_author=='PBEVV') then
+    call pbevv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+    return
+  elseif (vdw_author=='PBESVVS') then
+    call pbesvvs_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+    return 
   end if
 
   ns = min(nspin,2)
 
=== modified file 'Src/SiestaXC/vv_vdwxc.F90'
--- Src/SiestaXC/vv_vdwxc.F90	2016-05-10 10:45:36 +0000
+++ Src/SiestaXC/vv_vdwxc.F90	2019-09-06 14:47:07 +0000
@@ -993,3 +993,2520 @@
 end subroutine vv_vdw_theta
 
 END MODULE m_vv_vdwxc
+
+!!!==================SG4VV====================
+
+MODULE m_sg4vv_vdwxc
+
+! Used module procedures
+  use sys,          only: die              ! Termination routine
+  use mesh1D,       only: get_mesh         ! Returns the mesh points
+  use m_radfft,     only: radfft           ! Radial fast Fourier transform
+  use alloc,        only: re_alloc         ! Re-allocation routine
+  use mesh1D,       only: set_mesh         ! Sets a 1D mesh
+  use interpolation,only: spline           ! Sets spline interpolation
+  use interpolation,only: splint           ! Calculates spline interpolation
+
+! Used module parameters
+  use precision,    only: dp               ! Real double precision type
+
+#ifdef DEBUG_XC
+  use debugXC,      only: udebug           ! File unit for debug output
+!  use plot_module, only: plot
+#endif /* DEBUG_XC */
+
+  implicit none
+
+! Called by m_vdwxc
+PUBLIC ::            &
+  sg4vv_vdw_beta,       &! Returns parameter beta of the VV2010 functional
+  sg4vv_vdw_theta,      &! Finds function theta_ik(rho,grad_rho)
+  sg4vv_vdw_get_kmesh,  &! Returns size and values of (kf,kg) mesh
+  sg4vv_vdw_phi,        &! Finds and interpolates rho1*rho2*phi(k1,k2,k)
+  sg4vv_vdw_set_kcut     ! Sets the planewave cutoff kc of the integration grid
+
+#ifdef DEBUG_XC
+! Called by debugging test programs
+PUBLIC :: &
+  sg4vv_vdw_phi_val,    &! Kernel phi(kf1,kf2,kg1,kg2,r12) with no interpolation
+  sg4vv_vdw_phiofr       ! Kernel phi(k1,k2,r) at tabulated k1,k2-mesh values
+#endif /* DEBUG_XC */
+
+PRIVATE  ! Nothing is declared public beyond this point
+
+!! Set parameters of the SG4+rVV10m functional
+!! Ref. A. V. Terentjev et al, Computation 6, 7 (2018)
+  real(dp),parameter:: vv_C = 0.0001
+  real(dp),parameter:: vv_b = 12.0
+!!
+  real(dp),parameter:: vv_beta = 0.00497
+
+  ! Mesh parameters for table of phi(k1,k2,r) and its Fourier transform
+  character(len=*),parameter:: kernelPrefactor='rho' ! Prefactor of the
+                                             ! nonlocal kernel for interpolation
+                                             ! ('kf'|'kf2'|'sqr_rho'|'rho')
+  integer, parameter:: nr = 2048             ! Radial mesh points (power of 2)
+  real(dp),parameter:: rcut = 100._dp        ! Radial cutoff: r>rcut => phi=0
+  real(dp),parameter:: rmin = 1.e-6_dp       ! Min. radius as denominator
+  integer, parameter:: nkf =  7              ! Number of Fermi wavevectors
+  integer, parameter:: nkg =  5              ! Num. of grad(n)/n wavevectors
+  integer, parameter:: nkfg = nkf*nkg        ! Num. of (kf,kg) mesh points
+  real(dp),parameter:: kfcut = 5.0_dp        ! Max. Fermi wavevec.
+  real(dp),parameter:: kgcut = 5.0_dp        ! Max. grad(n)/n
+  real(dp),parameter:: dkfmaxdkfmin = 20.0_dp ! Last/first kf mesh interval
+  real(dp),parameter:: dkgmaxdkgmin = 2.0_dp  ! Last/first kg mesh interval
+  real(dp),parameter:: ktol = 1.e-12_dp      ! 'out of range' tolerance
+  real(dp),parameter:: amin = 1.e-12_dp      ! tiny denominator to avoid /0
+  real(dp),parameter:: rhoMin = 1.e-10_dp    ! Min. density for nonzero Vxc
+
+  ! Parameters for cutoff function, used in radial Fourier transforms of phi
+  integer, parameter:: ncut1 =  8      ! cutoff(x)=(1-x**ncut1)**ncut2
+  integer, parameter:: ncut2 =  4
+
+  ! Parameters for saturate function, used to enforce that q<qcut
+  integer, parameter:: nsat = 15 ! h(x,xc)=1/(1/x**nsat+1/xc**nsat)**(1/nsat)
+                                 ! nsat=15 approximately corresponds to 
+                                 ! nsat=12 for eq.(5) of Roman-Soler PRL
+
+  ! Private module variables and arrays
+  logical, save:: kcut_set=.false.         ! Has kcut been set?
+  logical, save:: kmesh_set=.false.        ! Has (kf,kg) mesh been set?
+  logical, save:: phi_table_set=.false.    ! Has phi_table been set?
+  integer, save:: nk                       ! Num. of radial mesh k-points
+  real(dp),save:: dr                       ! r-mesh interval
+  real(dp),save:: dk                       ! k-mesh interval
+  real(dp),save:: kcut                     ! Planewave cutoff: k>kcut => phi=0
+  real(dp),save:: kmax                     ! Max. k vector in Fourier transforms
+  real(dp),save:: kfmesh(nkf)              ! Mesh points for Fermi wavevector
+  real(dp),save:: kgmesh(nkg)              ! Mesh points for grad(n)/n
+  real(dp),pointer,save:: &
+                  phir(:,:,:)=>null(),    &! Table of phi(r)
+                  phik(:,:,:)=>null(),    &! Table of phi(k)
+                  d2phidr2(:,:,:)=>null(),&! Table of d^2(phi)/dr^2
+                  d2phidk2(:,:,:)=>null()  ! Table of d^2(phi)/dk^2
+
+CONTAINS
+
+! -----------------------------------------------------------------------------
+
+real(dp) function cutoff( x )
+
+! Defines a smooth cutoff function that falls from y(0)=1 to y(1)=0
+
+  implicit none
+  real(dp),intent(in):: x
+
+  if (x<=0._dp) then
+    cutoff = 1
+  else if (x>=1._dp) then
+    cutoff = 0
+  else
+    cutoff = (1-x**ncut1)**ncut2
+  end if
+
+end function cutoff
+
+!-----------------------------------------------------------------------------
+
+subroutine iofk( kf, kg, ikf, ikg )
+
+! Finds indexes ikf and ikg such that kfmesh(ikf) <= kf < kfmesh(ikf+1)
+! and kgmesh(ikg) <= kg < kgmesh(ikg+1) in logarithmic meshes of the form
+! k(ik) = b*(exp((ik-1)*a)-1)
+
+  implicit none
+  real(dp), intent(in) :: kf, kg
+  integer,  intent(out):: ikf, ikg
+
+  real(dp),save:: akf, akg, bkf, bkg
+  logical, save:: first_call = .true.
+  integer:: ik
+
+  ! Mesh data initializations
+  if (first_call) then
+    ! Find logarithmic-mesh a & b parameters: x(j)=x(1)+b*(exp(a*(j-1))-1)
+    akf = log( (kfmesh(nkf)-kfmesh(nkf-1)) / (kfmesh(2)-kfmesh(1)) ) / (nkf-2)
+    akg = log( (kgmesh(nkg)-kgmesh(nkg-1)) / (kgmesh(2)-kgmesh(1)) ) / (nkg-2)
+    akf = max( akf, amin )
+    akg = max( akg, amin )
+    bkf = (kfmesh(nkf) - kfmesh(1)) / (exp(akf*(nkf-1)) - 1)
+    bkg = (kgmesh(nkg) - kgmesh(1)) / (exp(akg*(nkg-1)) - 1)
+    ! Check that meshes are indeed logarithmic
+    do ik = 1,nkf
+      if (abs(kfmesh(ik)-kfmesh(1)-bkf*(exp(akf*(ik-1))-1))>ktol) &
+         call die('sg4vv_vdw_iofk ERROR: kfmesh not logarithmic')
+    end do
+    do ik = 1,nkg
+      if (abs(kgmesh(ik)-kgmesh(1)-bkg*(exp(akg*(ik-1))-1))>ktol) &
+         call die('sg4vv_vdw_iofk ERROR: kgmesh not logarithmic')
+    end do
+    first_call = .false.
+  end if
+
+  ! Check that kf and kg are within interpolation range
+  if (kf<kfmesh(1)-ktol .or. kf>kfmesh(nkf)+ktol .or. &
+      kg<kgmesh(1)-ktol .or. kg>kgmesh(nkg)+ktol) then
+     call die('sg4vv_vdw_iofk ERROR: (kf,kg) out of range')
+  endif
+
+  ! Find interpolation mesh intervals
+  ikf = 1 + int( log( 1 + (kf-kfmesh(1))/bkf ) / akf )
+  ikg = 1 + int( log( 1 + (kg-kgmesh(1))/bkg ) / akg )
+  ikf = max( 1, ikf )
+  ikg = max( 1, ikg )
+  ikf = min( nkf-1, ikf )
+  ikg = min( nkg-1, ikg )
+
+end subroutine iofk
+
+! -----------------------------------------------------------------------------
+
+subroutine kofn( n, gn, kf, kg, dkfdn, dkgdn, dkfdgn, dkgdgn )
+
+! Finds Fermi and gradient wavevectors from density and gradient
+
+  implicit none
+  real(dp), intent(in) :: n          ! Electron density
+  real(dp), intent(in) :: gn(3)      ! Density gradient
+  real(dp), intent(out):: kf         ! Local Fermi wavevector
+  real(dp), intent(out):: kg         ! |grad(n)|/n
+  real(dp), intent(out):: dkfdn      ! dkf/dn
+  real(dp), intent(out):: dkgdn      ! dkg/dn
+  real(dp), intent(out):: dkfdgn(3)  ! dkf/dgn
+  real(dp), intent(out):: dkgdgn(3)  ! dkg/dgn
+
+  real(dp):: gn2, pi
+
+! Trap exception for zero density
+  if (n <= 1.e-30_dp) then
+    kf = 0
+    kg = 0
+    dkfdn = 0
+    dkfdgn = 0
+    dkgdn = 0
+    dkgdgn = 0
+    return
+  end if
+
+! Find kf and kg
+  pi = acos(-1._dp)
+  kf = (3*pi**2*n)**(1._dp/3)
+  gn2 = sum(gn**2)
+  kg = sqrt(gn2)/n
+
+! Find derivatives
+  dkfdn = kf/n/3
+  dkfdgn = 0
+  dkgdn = -kg/n
+  dkgdgn = kg*gn/gn2
+
+end subroutine kofn
+
+!-----------------------------------------------------------------------------
+
+subroutine pofk( kf, kg, p, dpdkf, dpdkg )
+
+! Finds the values and derivatives, at (kf,kg), of the bicubic polynomials 
+! p_i(kf,kg) such that
+!    y(kf,kg) = Sum_i p_i(kf,kg) * y_i
+! is the bicubic spline interpolation at (kf,kg) of (any) function y with
+! values y_i at mesh points (kfmesh,kgmesh)_i
+
+  implicit none
+  real(dp),intent(in) :: kf, kg ! point at which the polynomials are required
+  real(dp),intent(out):: p(nkfg)      ! polynomial values at (kf,kg)
+  real(dp),intent(out):: dpdkf(nkfg)  ! dp/dkf at (kf,kg)
+  real(dp),intent(out):: dpdkg(nkfg)  ! dp/dkg at (kf,kg)
+
+  integer :: ikf, ikg, ikfg, ikf0, ikg0
+  real(dp):: a, b, dk, dkf0dkf, dkg0dkg, kf0, kg0
+  real(dp):: pkf0(nkf), dpkfdkf0(nkf), pkg0(nkg), dpkgdkg0(nkg)
+  logical, save :: first_call=.true.
+  real(dp),save :: pkf(nkf,nkf), d2pkfdkf2(nkf,nkf)
+  real(dp),save :: pkg(nkg,nkg), d2pkgdkg2(nkg,nkg)
+
+! Set up spline polynomial basis
+  if (first_call) then
+    do ikf = 1,nkf
+      pkf(:,ikf) = 0    ! ikf'th polynomial pkf(:,ikf) is one at kfmesh(ikf) 
+      pkf(ikf,ikf) = 1  ! and zero at all other points
+      call spline( kfmesh, pkf(:,ikf), nkf, 1.e30_dp, 1.e30_dp, &
+                   d2pkfdkf2(:,ikf) )
+!      call spline( kfmesh, pkf(:,ikf), nkf, 0._dp, 0._dp, d2pkfdkf2(:,ikf) )
+    end do
+    do ikg = 1,nkg
+      pkg(:,ikg) = 0
+      pkg(ikg,ikg) = 1
+      call spline( kgmesh, pkg(:,ikg), nkg, 1.e30_dp, 1.e30_dp, &
+                   d2pkgdkg2(:,ikg) )
+!      call spline( kgmesh, pkg(:,ikg), nkg, 0._dp, 0._dp, d2pkgdkg2(:,ikg) )
+    end do
+
+!   DEBUG
+    open(22,file='pkf.dat')
+    do ikf = 1,nkf
+      write(22,'(20e15.6)') kfmesh(ikf), pkf(:,ikf), d2pkfdkf2(:,ikf)
+    end do
+    close(22)
+    open(22,file='pkg.dat')
+    do ikg = 1,nkg
+      write(22,'(20e15.6)') kgmesh(ikg), pkg(:,ikg), d2pkgdkg2(:,ikg)
+    end do
+    close(22)
+!   END DEBUG
+
+    first_call = .false.
+  end if
+
+! 'Saturate' (kf,kg) values to bring them to interpolation range
+  call saturate( kf, kfcut, kf0, dkf0dkf )
+  call saturate( kg, kgcut, kg0, dkg0dkg )
+
+! Find interval of k mesh in which (kf0,kg0) point is included
+  call iofk( kf0, kg0, ikf0, ikg0 )
+
+! Evaluate pkf polynomials of spline basis at kf0
+! The splint code is in-lined here because it is a hot point
+  dk = kfmesh(ikf0+1) - kfmesh(ikf0)
+  a = (kfmesh(ikf0+1) - kf0) / dk   ! dadkf0 = -1/dk
+  b = (kf0 - kfmesh(ikf0)) / dk     ! dbdkf0 = +1/dk
+  do ikf = 1,nkf
+    pkf0(ikf) = a*pkf(ikf0,ikf) + b*pkf(ikf0+1,ikf) &
+              + ((a**3-a)*d2pkfdkf2(ikf0,ikf) +     &
+                 (b**3-b)*d2pkfdkf2(ikf0+1,ikf)) * dk**2/6
+    dpkfdkf0(ikf) = - (pkf(ikf0,ikf) - pkf(ikf0+1,ikf)) / dk &
+                    - ((3*a**2-1)*d2pkfdkf2(ikf0,ikf) -     &
+                       (3*b**2-1)*d2pkfdkf2(ikf0+1,ikf)) * dk/6
+  end do
+
+! Evaluate pkg polynomials of spline basis at kg0
+  dk = kgmesh(ikg0+1) - kgmesh(ikg0)
+  a = (kgmesh(ikg0+1) - kg0) / dk   ! dadkg0 = -1/dk
+  b = (kg0 - kgmesh(ikg0)) / dk     ! dbdkg0 = +1/dk
+  do ikg = 1,nkg
+    pkg0(ikg) = a*pkg(ikg0,ikg) + b*pkg(ikg0+1,ikg) &
+              + ((a**3-a)*d2pkgdkg2(ikg0,ikg) +     &
+                 (b**3-b)*d2pkgdkg2(ikg0+1,ikg)) * dk**2/6
+    dpkgdkg0(ikg) = - (pkg(ikg0,ikg) - pkg(ikg0+1,ikg)) / dk &
+                    - ((3*a**2-1)*d2pkgdkg2(ikg0,ikg) -     &
+                       (3*b**2-1)*d2pkgdkg2(ikg0+1,ikg)) * dk/6
+  end do
+
+! Evaluate pkf*pkg polynomials at (kf0,kg0)
+  ikfg = 0
+  do ikg = 1,nkg
+    do ikf = 1,nkf
+      ikfg = ikfg+1
+      p(ikfg) = pkf0(ikf) * pkg0(ikg)
+      dpdkf(ikfg) = dpkfdkf0(ikf)*dkf0dkf * pkg0(ikg)
+      dpdkg(ikfg) = pkf0(ikf) * dpkgdkg0(ikg)*dkg0dkg
+    end do
+  end do
+
+end subroutine pofk
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate( x, xc, y, dydx )
+
+  ! Defines a function y(x,xc) = 1/(1/x^nsat+1/xc^nsat)^(1/nsat), where nsat
+  ! is an integer set in the module header. This function is approximately
+  ! equal to x for x<xc and it saturates to xc when x->infinity
+
+  implicit none
+  real(dp),intent(in) :: x     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: y     ! Function value
+  real(dp),intent(out):: dydx  ! Derivative dy/dx
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('sg4vv_vdwxc_saturate ERROR: xc<=0')
+  else if (x<0._dp) then
+    call die('sg4vv_vdwxc_saturate ERROR: x<0')
+  else if (x==0._dp) then
+    y = 0
+    dydx = 1
+  else
+    z = 1/x**nsat + 1/xc**nsat
+    y = 1/z**(1._dp/nsat)
+    dydx = y/z/x**(nsat+1)
+  end if
+
+end subroutine saturate
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate_inverse( y, xc, x, dxdy )
+
+! Finds the inverse of the function defined in saturate subroutine:
+!   y=1/(1/x^n+1/xc^n)^(1/n)  =>  x=1/(1/y^n-1/xc^n)^(1/n)
+
+  implicit none
+  real(dp),intent(in) :: y     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: x     ! Inverse function value
+  real(dp),intent(out):: dxdy  ! Derivative dx/dy
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('sg4vv_vdwxc_saturate_inverse ERROR: xc<=0')
+  else if (y<0._dp .or. y>=xc) then
+    call die('sg4vv_vdwxc_saturate_inverse ERROR: y out of range')
+  else if (y==0._dp) then
+    x = 0
+    dxdy = 1
+  else
+    z = 1/y**nsat - 1/xc**nsat
+    x = 1/z**(1._dp/nsat)
+    dxdy = x/z/y**(nsat+1)
+  end if
+
+end subroutine saturate_inverse
+
+!-----------------------------------------------------------------------------
+
+subroutine set_kmesh()
+
+! Sets mesh of q values
+
+  implicit none
+  integer :: mkf, mkg
+
+  if (.not.kmesh_set) then
+    call set_mesh( nkf, xmax=kfcut, dxndx1=dkfmaxdkfmin )
+    call get_mesh( nkf, mkf, kfmesh )
+    call set_mesh( nkg, xmax=kgcut, dxndx1=dkgmaxdkgmin )
+    call get_mesh( nkg, mkg, kgmesh )
+    kmesh_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(/,a,/,(10f8.4))') 'sg4vv_vdw_set_kmesh: kfmesh =', kfmesh
+    write(udebug,'(/,a,/,(10f8.4))') 'sg4vv_vdw_set_kmesh: kgmesh =', kgmesh
+#endif /* DEBUG_XC */
+  end if
+
+end subroutine set_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine set_phi_table()
+
+! Finds and stores in memory the interpolation table (mesh points and 
+! function values) for the kernel phi(k1,k2,k).
+
+  implicit none
+  character(len=*),parameter:: myName = 'sg4vv_vdwxc/set_phi_table '
+  integer :: ik, ikf1, ikf2, ikg1, ikg2, ik1, ik2, ir
+  real(dp):: dkdk0, dphidk0, dphidkmax, dphidr0, dphidrmax, &
+             k, kf1, kf2, kg1, kg2, pi, r(0:nr)
+
+! Check that table was not set yet
+  if (phi_table_set) return
+
+! Check that kf, kg, r, and k meshes have been set
+  if (.not.kmesh_set) call set_kmesh()
+  if (.not.kcut_set) call die('sg4vv_vdw_set_phi_table ERROR: kcut not set')
+  forall(ir=0:nr) r(ir) = ir*dr
+  pi = acos(-1.0_dp)
+
+! Allocate arrays
+  call re_alloc( phir,     0,nr, 1,nkfg, 1,nkfg, myName//'phir' )
+  call re_alloc( phik,     0,nr, 1,nkfg, 1,nkfg, myName//'phik' )
+  call re_alloc( d2phidr2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidr2' )
+  call re_alloc( d2phidk2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidk2' )
+
+! Loop on (k1,k2) mesh points
+  do ikg2 = 1,nkg                    ! loop on kg2
+    do ikf2 = 1,nkf                  ! loop on kf2
+      ik2 = ikf2 + nkf*(ikg2-1)      ! combined (ikf2,ikg2) index
+      do ikg1 = 1,nkg                ! loop on kg1
+        do ikf1 = 1,nkf              ! loop on kf1
+          ik1 = ikf1 + nkf*(ikg1-1)  ! combined (ikf1,ikg1) index
+          if (ik1>ik2) cycle         ! since we will symmetrize at the end
+
+          ! Saturated (kf,kg) values
+          kf1 = kfmesh(ikf1)
+          kf2 = kfmesh(ikf2)
+          kg1 = kgmesh(ikg1)
+          kg2 = kgmesh(ikg2)
+
+          ! Find original (unsaturated) kf anf kg values
+          ! call saturate_inverse( kfmesh(ikf1), kfcut, kf1, dkdk0 )
+          ! call saturate_inverse( kfmesh(ikf2), kfcut, kf2, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg1), kgcut, kg1, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg2), kgcut, kg2, dkdk0 )
+
+          ! Find kernel as a function of r12
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = sg4vv_vdw_phi_val( kf1, kf2, kg1, kg2, r(ir) )
+          end do
+
+          ! Kill kernel smoothly at r=rcut
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = phir(ir,ik1,ik2) * cutoff( r(ir)/rcut )
+          end do
+
+          ! Find kernel in reciprocal space
+          call radfft( 0, nr, rcut, phir(:,ik1,ik2), phik(:,ik1,ik2) )
+          phik(:,ik1,ik2) = phik(:,ik1,ik2) * (2*pi)**1.5_dp
+
+          ! Filter out above kcut
+          phik(nk:nr,ik1,ik2) = 0
+
+          ! Soft filter below kcut
+          do ik = 1,nk
+            k = ik * dk
+            phik(ik,ik1,ik2) = phik(ik,ik1,ik2) * cutoff(k/kcut)
+          end do
+
+          ! Find filtered kernel in real space
+          call radfft( 0, nr, kmax, phik(:,ik1,ik2), phir(:,ik1,ik2) )
+          phir(:,ik1,ik2) = phir(:,ik1,ik2) / (2*pi)**1.5_dp
+
+          ! Set up spline interpolation tables
+          dphidr0 = 0      ! since phi(k1,k2,r) is even in r
+          dphidk0 = 0      ! and therefore phi(k1,k2,k) is also even in k
+          dphidrmax = 0    ! since phi->0 for r->infty
+          dphidkmax = 0    ! and also when k->infty
+          call spline( dr, phir(:,ik1,ik2), nr+1, dphidr0, dphidrmax, &
+                       d2phidr2(:,ik1,ik2) )
+          call spline( dk, phik(:,ik1,ik2), nk+1, dphidk0, dphidkmax, &
+                       d2phidk2(:,ik1,ik2) )
+
+          ! Fill symmetric elements
+          phir(:,ik2,ik1) = phir(:,ik1,ik2)
+          phik(:,ik2,ik1) = phik(:,ik1,ik2)
+          d2phidr2(:,ik2,ik1) = d2phidr2(:,ik1,ik2)
+          d2phidk2(:,ik2,ik1) = d2phidk2(:,ik1,ik2)
+
+!#ifdef DEBUG_XC
+!          if (.false. .and. ik1==ik2) then
+!            print*, 'sg4vv_vdw_set_kcut: ik1,ik2=', ik1, ik2
+!            call window( 0._dp, 5._dp, -1._dp, 4._dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 1._dp )
+!            call plot( nr+1, r, phi, phir(:,ik1,ik2) )
+!            call window( 0._dp, 10._dp, -0.05_dp, 0.15_dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 0.05_dp )
+!            call plot( nr+1, r, r**2*phi, r**2*phir(:,ik1,ik2))
+!            call show()
+!          end if
+!#endif /* DEBUG_XC */
+
+        end do ! ikf1
+      end do ! ikg1
+    end do ! ikf2
+  end do ! ikg2
+
+!#ifdef DEBUG_XC
+!  open(17,file='vv_phi.table')
+!  write(17,'(3a6,2a12,/,(3i6,2f15.9))') &
+!    'ik1', 'ik2', 'ir', 'phi', 'd2phi/dk2', &
+!!    (((ik1,ik2,ir,phir(ir,ik1,ik2),d2phidr2(ir,ik1,ik2), &
+!    (((ik1,ik2,ir,phik(ir,ik1,ik2),d2phidk2(ir,ik1,ik2), &
+!       ir=0,100),ik2=1,nkfg),ik1=1,nkfg)
+!  close(17)
+!#endif /* DEBUG_XC */
+
+! Mark table as set
+  phi_table_set = .true.
+
+end subroutine set_phi_table
+
+! -----------------------------------------------------------------------------
+
+real(dp) function sg4vv_vdw_phi_val( kf1, kf2, kg1, kg2, r12 )
+
+! vdW energy kernel of Vydrov-vanVoorhis, eq.(2) of JCP 133, 244103 (2010)
+! This subroutine calculates the 'raw' kernel directly, without interpolarions
+! In practice, it returns n1*n2*phi(k1,k2,r12), which is smooth for kf->0
+! Input:
+!   real(dp):: kf1, kf2 ! Fermi wavevectors at points 1 and 2
+!   real(dp):: kg1, kg2 ! |grad(n)|/n at points 1 and 2
+!   real(dp):: r12      ! distance between points 1 and 2
+
+! Arguments
+  implicit none
+  real(dp),intent(in) :: kf1, kf2, kg1, kg2, r12
+
+! Internal variables
+  real(dp):: g1, g2, kappa1, kappa2, kf1m, kf2m, n1, n2, &
+             phi, pi, w01, w02, wg1, wg2, wp1, wp2
+
+! Avoid dividing by zero when kf=0 
+  kf1m = max(kf1,ktol)
+  kf2m = max(kf2,ktol)
+
+! Find kernel
+  pi = acos(-1.0_dp)
+  n1 = kf1m**3/(3*pi**2)          ! electron density at point 1
+  n2 = kf2m**3/(3*pi**2)          ! electron density at point 2
+  wp1 = sqrt(4*pi*n1)             ! local plasma frequency at point 1
+  wp2 = sqrt(4*pi*n2)             ! local plasma frequency at point 2
+  wg1 = sqrt(vv_C*kg1**4)         ! local band gap at point 1
+  wg2 = sqrt(vv_C*kg2**4)         ! local band gap at point 2
+  kappa1 = vv_b*kf1m**2/wp1       ! local VV kappa variable (eq.(9)) at point 1
+  kappa2 = vv_b*kf2m**2/wp2       ! kappa variable at point 2
+  w01 = sqrt(wg1**2+wp1**2/3)     ! local w0 frequency (eq.(5)) at point 1
+  w02 = sqrt(wg2**2+wp2**2/3)     ! local w0 frequency at point 2
+  g1 = w01*r12**2 + kappa1        ! local g variable (eq.(3)) at point 1
+  g2 = w02*r12**2 + kappa2        ! local g variable at point 2
+  phi = -1.5_dp/g1/g2/(g1+g2)     ! VV kernel phi (eq.(2))
+
+  if (kernelPrefactor=='rho') then
+    ! Return whole integrand of eq.(1)
+    sg4vv_vdw_phi_val = n1*n2*phi
+  else if (kernelPrefactor=='kf') then
+    ! Return kf1*kf2*phi
+    sg4vv_vdw_phi_val = kf1*kf2*phi
+  else if (kernelPrefactor=='kf2') then
+    ! Return (kf1*kf2)**2*phi
+    sg4vv_vdw_phi_val = (kf1*kf2)**2*phi
+  else if (kernelPrefactor=='sqr_rho') then
+    ! Find and return sqrt(n1*n2)*phi
+    sg4vv_vdw_phi_val = sqrt(n1*n2)*phi
+  else
+    call die('sg4vv_vdw_phi_val ERROR: unknown kernelPrefactor')
+  end if
+
+end function sg4vv_vdw_phi_val
+
+! -----------------------------------------------------------------------------
+
+real(dp) function sg4vv_vdw_beta()
+
+! Returns parameter beta of VV functional
+
+  implicit none
+  sg4vv_vdw_beta = vv_beta
+
+end function sg4vv_vdw_beta
+
+!-----------------------------------------------------------------------------
+
+subroutine sg4vv_vdw_get_kmesh( mkf, mkg, kf, kg )
+
+! Returns size and values of (kf,kg) mesh
+
+  implicit none
+  integer,          intent(out) :: mkf    ! Number of kf mesh points
+  integer,          intent(out) :: mkg    ! Number of kg mesh points
+  real(dp),optional,intent(out) :: kf(:)  ! Values of kf mesh points
+  real(dp),optional,intent(out) :: kg(:)  ! Values of kg mesh points
+  integer:: nmax
+  if (.not.kmesh_set) call set_kmesh()
+  mkf = nkf
+  mkg = nkg
+  if (present(kf)) then
+    nmax = max( nkf, size(kf) )
+    kf(1:nmax) = kfmesh(1:nmax)
+  end if
+  if (present(kg)) then
+    nmax = max( nkg, size(kg) )
+    kg(1:nmax) = kgmesh(1:nmax)
+  end if
+end subroutine sg4vv_vdw_get_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine sg4vv_vdw_phi( k, phi, dphidk )
+
+! Finds by interpolation phi(k1,k2,k) (Fourier transform of phi(k1,k2,r)),
+! with k1=(kf1,kg1), k2=(kf2,kg2) for all mesh values of kf (Fermi 
+! wavevector) and kg (grad(n)/n). If the interpolation table does not exist,
+! it is calculated in the first call to sg4vv_vdw_phi. It requires a previous 
+! call to sg4vv_vdw_set_kc to set k mesh.
+
+  implicit none
+  real(dp),intent(in) :: k            ! Modulus of actual k vector
+  real(dp),intent(out):: phi(:,:)     ! phi(k1,k2,k) at given k
+                                      ! for all k1,k2 in (kf,kg) mesh
+  real(dp),intent(out):: dphidk(:,:)  ! dphi(k1,k2,k)/dk at given k
+
+  integer :: ik, ik1, ik2
+  real(dp):: a, a2, a3, b, b2, b3
+
+! Set interpolation table
+  if (.not.phi_table_set) call set_phi_table()
+
+! Check argument sizes
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    call die('sg4vv_vdw_phi: ERROR: size(phi) too small')
+
+! Find phi values at point k
+  if (k >= kcut) then
+    phi(:,:) = 0
+    dphidk(:,:) = 0
+  else
+    ! Expand interpolation inline since this is the hottest point in VdW
+    ik = int(k/dk)
+    a = ((ik+1)*dk-k)/dk
+    b = 1 - a
+    a2 = (3*a**2 -1) * dk / 6
+    b2 = (3*b**2 -1) * dk / 6
+    a3 = (a**3 - a) * dk**2 / 6
+    b3 = (b**3 - b) * dk**2 / 6
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+!        call splint( dk, phik(:,ik1,ik2), d2phidk2(:,ik1,ik2), &
+!                     nk+1, k, phi(ik1,ik2), dphidk(ik1,ik2), pr )
+        phi(ik1,ik2) = a*phik(ik,ik1,ik2) + b*phik(ik+1,ik1,ik2) &
+                + a3*d2phidk2(ik,ik1,ik2) + b3*d2phidk2(ik+1,ik1,ik2)
+        dphidk(ik1,ik2) = (-phik(ik,ik1,ik2) &
+                               +phik(ik+1,ik1,ik2) )/dk &
+                - a2*d2phidk2(ik,ik1,ik2) + b2*d2phidk2(ik+1,ik1,ik2)
+        phi(ik2,ik1) = phi(ik1,ik2)
+        dphidk(ik2,ik1) = dphidk(ik1,ik2)
+      end do
+    end do
+  end if
+
+end subroutine sg4vv_vdw_phi
+
+!-----------------------------------------------------------------------------
+
+subroutine sg4vv_vdw_phiofr( r, phi )
+
+! Finds phi(k1,k2,r) with k1=(kf1,kg1), k2=(kf2,kg2) for mesh values of
+! kf's (Fermi wavevectors) and kg's (grad(n)/n)
+
+  implicit none
+  real(dp),intent(in) :: r
+  real(dp),intent(out):: phi(:,:)
+
+  integer :: ikf1, ikf2, ikg1, ikg2, ik1, ik2
+  real(dp):: dphidr
+
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    stop 'vv_phiofr: ERROR: size(phi) too small'
+  if (.not.phi_table_set) call set_phi_table()
+
+  if (r >= rcut) then
+    phi(:,:) = 0
+  else
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+        call splint( dr, phir(:,ik1,ik2), d2phidr2(:,ik1,ik2), &
+                     nr+1, r, phi(ik1,ik2), dphidr )
+        phi(ik2,ik1) = phi(ik1,ik2)
+      end do ! ik1
+    end do ! ik2
+  end if ! (r>=rcut)
+
+end subroutine sg4vv_vdw_phiofr
+
+!-----------------------------------------------------------------------------
+
+subroutine sg4vv_vdw_set_kcut( kc )
+
+! Sets the reciprocal planewave cutoff kc of the integration grid, and finds 
+! the interpolation table to be used by vdw_phi to obtain the vdW kernel phi
+! at the reciprocal mesh wavevectors.
+
+  implicit none
+  real(dp),intent(in):: kc  ! Planewave cutoff: k>kcut => phi=0
+
+  real(dp):: pi
+
+  ! Set kcut and radial mesh parameters, if not already set
+  if (.not. kc==kcut) then
+    kcut = kc
+    pi = acos(-1._dp)
+    dr = rcut / nr
+    dk = pi / rcut
+    kmax = pi / dr
+    nk = int(kcut/dk) + 1
+    if (nk>nr) stop 'sg4vv_vdw_set_kcut: ERROR: nk>nr'
+    kcut_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(a,5f8.3)') 'sg4vv_vdw_set_kcut: kfcut,kgcut,rcut,kcut,kmax=', &
+      kfcut, kgcut, rcut, kc, kmax
+#endif /* DEBUG_XC */
+  end if
+
+  ! Set (kf,kg) mesh and phi table
+  call set_kmesh()
+  call set_phi_table()
+
+end subroutine sg4vv_vdw_set_kcut
+
+!-----------------------------------------------------------------------------
+
+subroutine sg4vv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+
+! Finds the value and derivatives of theta_i(rho,grad_rho) of eq.(8)
+! of Roman-Soler PRL 2009. In practice, theta_i=p_i rather than rho*p_i,
+! beacuse p_i is used here to expand rho1*rho2*phi rather than only phi.
+
+  implicit none
+  integer, intent(in) :: nspin                 ! Number of spin components
+  real(dp),intent(in) :: rhos(nspin)           ! Electron spin density
+  real(dp),intent(in) :: grhos(3,nspin)        ! Spin density gradient
+  real(dp),intent(out):: theta(nkfg)           ! Exp. polynomials at (kf,kg)
+  real(dp),intent(out):: dtdrho(nkfg,nspin)    ! dtheta(ik)/drhos
+  real(dp),intent(out):: dtdgrho(3,nkfg,nspin) ! dtheta(ik)/dgrhos(ix)
+
+  integer :: is, ix, ns
+  real(dp):: rho, grho(3), dpdkf(nkfg), dpdkg(nkfg), dkfdrho, dkfdgrho(3), &
+             dkgdrho, dkgdgrho(3), kf, kg, p(nkfg)
+
+  ! Sum spin components of electron density
+  ns = min(nspin,2)     ! num. of diagonal spin components (if noncollinear)
+  rho = sum(rhos(1:ns))            ! local electron density
+  grho = sum(grhos(:,1:ns),dim=2)  ! local density gradient
+
+  ! If density is between threshold, return zero
+  if (rho<rhoMin) then
+    theta = 0
+    dtdrho = 0
+    dtdgrho = 0
+    return
+  end if
+
+  ! Find local Fermi and gradient wavevectors, and their derivatives
+  call kofn( rho, grho, kf, kg, dkfdrho, dkgdrho, dkfdgrho, dkgdgrho )
+
+  ! Find expansion polynomials of integrand kernel of nonlocal vdW energy
+  call pofk( kf, kg, p, dpdkf, dpdkg )
+
+  ! Find theta functions and their derivatives with respect to rho and grho
+  if (kernelPrefactor=='rho') then
+    ! This is the right code if sg4vv_vdw_phi_val returns rho1*rho2*phi
+    theta(1:nkfg) = p(1:nkfg)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = dpdkf(1:nkfg)*dkfdrho + dpdkg(1:nkfg)*dkgdrho
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix)
+      end do
+    end do
+  else if (kernelPrefactor=='kf') then
+    ! This is the right code if sg4vv_vdw_phi_val returns kf1*kf2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf &
+                        + p(1:nkfg)*(1/kf-rho/kf**2*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf
+      end do
+    end do
+  else if (kernelPrefactor=='kf2') then
+    ! This is the right code if sg4vv_vdw_phi_val returns (kf1*kf2)**2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf**2
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf**2 &
+                        + p(1:nkfg)*(1/kf**2-2*rho/kf**3*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf**2
+      end do
+    end do
+  else if (kernelPrefactor=='sqr_rho') then
+    ! This is the right code if sg4vv_vdw_phi_val returns sqrt(rho1*rho2)*phi
+    rho = max(rho,1.e-12_dp)  ! to avoid division by zero
+    theta(1:nkfg) = p(1:nkfg) * sqrt(rho)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = ( dpdkf(1:nkfg)*dkfdrho &
+                          + dpdkg(1:nkfg)*dkgdrho ) * sqrt(rho) &
+                        + p(1:nkfg) * 0.5/sqrt(rho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) * sqrt(rho) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix) * sqrt(rho)
+      end do
+    end do
+  else
+    call die('sg4vv_vdw_theta ERROR: unknown kernelPrefactor')
+  end if ! (kernelPrefactor=='rho')
+
+end subroutine sg4vv_vdw_theta
+
+END MODULE m_sg4vv_vdwxc
+
+!!!=============== PBE+VV10 ========================
+
+MODULE m_pbevv_vdwxc
+
+! Used module procedures
+  use sys,          only: die              ! Termination routine
+  use mesh1D,       only: get_mesh         ! Returns the mesh points
+  use m_radfft,     only: radfft           ! Radial fast Fourier transform
+  use alloc,        only: re_alloc         ! Re-allocation routine
+  use mesh1D,       only: set_mesh         ! Sets a 1D mesh
+  use interpolation,only: spline           ! Sets spline interpolation
+  use interpolation,only: splint           ! Calculates spline interpolation
+
+! Used module parameters
+  use precision,    only: dp               ! Real double precision type
+
+#ifdef DEBUG_XC
+  use debugXC,      only: udebug           ! File unit for debug output
+!  use plot_module, only: plot
+#endif /* DEBUG_XC */
+
+  implicit none
+
+! Called by m_vdwxc
+PUBLIC ::            &
+  pbevv_vdw_beta,       &! Returns parameter beta of the VV2010 functional
+  pbevv_vdw_theta,      &! Finds function theta_ik(rho,grad_rho)
+  pbevv_vdw_get_kmesh,  &! Returns size and values of (kf,kg) mesh
+  pbevv_vdw_phi,        &! Finds and interpolates rho1*rho2*phi(k1,k2,k)
+  pbevv_vdw_set_kcut     ! Sets the planewave cutoff kc of the integration grid
+
+#ifdef DEBUG_XC
+! Called by debugging test programs
+PUBLIC :: &
+  pbevv_vdw_phi_val,    &! Kernel phi(kf1,kf2,kg1,kg2,r12) with no interpolation
+  pbevv_vdw_phiofr       ! Kernel phi(k1,k2,r) at tabulated k1,k2-mesh values
+#endif /* DEBUG_XC */
+
+PRIVATE  ! Nothing is declared public beyond this point
+
+  ! Set parameters of the PBE+rVV10,
+  ! Ref. H. Peng and J. P. Perdew, Phys. Rev. B  95, 0811105 (2017)
+  real(dp),parameter:: vv_C = 0.0093
+  real(dp),parameter:: vv_b = 10.0
+  real(dp),parameter:: vv_beta = 0.00497
+
+  ! Mesh parameters for table of phi(k1,k2,r) and its Fourier transform
+  character(len=*),parameter:: kernelPrefactor='rho' ! Prefactor of the
+                                             ! nonlocal kernel for interpolation
+                                             ! ('kf'|'kf2'|'sqr_rho'|'rho')
+  integer, parameter:: nr = 2048             ! Radial mesh points (power of 2)
+  real(dp),parameter:: rcut = 100._dp        ! Radial cutoff: r>rcut => phi=0
+  real(dp),parameter:: rmin = 1.e-6_dp       ! Min. radius as denominator
+  integer, parameter:: nkf =  7              ! Number of Fermi wavevectors
+  integer, parameter:: nkg =  5              ! Num. of grad(n)/n wavevectors
+  integer, parameter:: nkfg = nkf*nkg        ! Num. of (kf,kg) mesh points
+  real(dp),parameter:: kfcut = 5.0_dp        ! Max. Fermi wavevec.
+  real(dp),parameter:: kgcut = 5.0_dp        ! Max. grad(n)/n
+  real(dp),parameter:: dkfmaxdkfmin = 20.0_dp ! Last/first kf mesh interval
+  real(dp),parameter:: dkgmaxdkgmin = 2.0_dp  ! Last/first kg mesh interval
+  real(dp),parameter:: ktol = 1.e-12_dp      ! 'out of range' tolerance
+  real(dp),parameter:: amin = 1.e-12_dp      ! tiny denominator to avoid /0
+  real(dp),parameter:: rhoMin = 1.e-10_dp    ! Min. density for nonzero Vxc
+
+  ! Parameters for cutoff function, used in radial Fourier transforms of phi
+  integer, parameter:: ncut1 =  8      ! cutoff(x)=(1-x**ncut1)**ncut2
+  integer, parameter:: ncut2 =  4
+
+  ! Parameters for saturate function, used to enforce that q<qcut
+  integer, parameter:: nsat = 15 ! h(x,xc)=1/(1/x**nsat+1/xc**nsat)**(1/nsat)
+                                 ! nsat=15 approximately corresponds to 
+                                 ! nsat=12 for eq.(5) of Roman-Soler PRL
+
+  ! Private module variables and arrays
+  logical, save:: kcut_set=.false.         ! Has kcut been set?
+  logical, save:: kmesh_set=.false.        ! Has (kf,kg) mesh been set?
+  logical, save:: phi_table_set=.false.    ! Has phi_table been set?
+  integer, save:: nk                       ! Num. of radial mesh k-points
+  real(dp),save:: dr                       ! r-mesh interval
+  real(dp),save:: dk                       ! k-mesh interval
+  real(dp),save:: kcut                     ! Planewave cutoff: k>kcut => phi=0
+  real(dp),save:: kmax                     ! Max. k vector in Fourier transforms
+  real(dp),save:: kfmesh(nkf)              ! Mesh points for Fermi wavevector
+  real(dp),save:: kgmesh(nkg)              ! Mesh points for grad(n)/n
+  real(dp),pointer,save:: &
+                  phir(:,:,:)=>null(),    &! Table of phi(r)
+                  phik(:,:,:)=>null(),    &! Table of phi(k)
+                  d2phidr2(:,:,:)=>null(),&! Table of d^2(phi)/dr^2
+                  d2phidk2(:,:,:)=>null()  ! Table of d^2(phi)/dk^2
+
+CONTAINS
+
+! -----------------------------------------------------------------------------
+
+real(dp) function cutoff( x )
+
+! Defines a smooth cutoff function that falls from y(0)=1 to y(1)=0
+
+  implicit none
+  real(dp),intent(in):: x
+
+  if (x<=0._dp) then
+    cutoff = 1
+  else if (x>=1._dp) then
+    cutoff = 0
+  else
+    cutoff = (1-x**ncut1)**ncut2
+  end if
+
+end function cutoff
+
+!-----------------------------------------------------------------------------
+
+subroutine iofk( kf, kg, ikf, ikg )
+
+! Finds indexes ikf and ikg such that kfmesh(ikf) <= kf < kfmesh(ikf+1)
+! and kgmesh(ikg) <= kg < kgmesh(ikg+1) in logarithmic meshes of the form
+! k(ik) = b*(exp((ik-1)*a)-1)
+
+  implicit none
+  real(dp), intent(in) :: kf, kg
+  integer,  intent(out):: ikf, ikg
+
+  real(dp),save:: akf, akg, bkf, bkg
+  logical, save:: first_call = .true.
+  integer:: ik
+
+  ! Mesh data initializations
+  if (first_call) then
+    ! Find logarithmic-mesh a & b parameters: x(j)=x(1)+b*(exp(a*(j-1))-1)
+    akf = log( (kfmesh(nkf)-kfmesh(nkf-1)) / (kfmesh(2)-kfmesh(1)) ) / (nkf-2)
+    akg = log( (kgmesh(nkg)-kgmesh(nkg-1)) / (kgmesh(2)-kgmesh(1)) ) / (nkg-2)
+    akf = max( akf, amin )
+    akg = max( akg, amin )
+    bkf = (kfmesh(nkf) - kfmesh(1)) / (exp(akf*(nkf-1)) - 1)
+    bkg = (kgmesh(nkg) - kgmesh(1)) / (exp(akg*(nkg-1)) - 1)
+    ! Check that meshes are indeed logarithmic
+    do ik = 1,nkf
+      if (abs(kfmesh(ik)-kfmesh(1)-bkf*(exp(akf*(ik-1))-1))>ktol) &
+         call die('pbevv_vdw_iofk ERROR: kfmesh not logarithmic')
+    end do
+    do ik = 1,nkg
+      if (abs(kgmesh(ik)-kgmesh(1)-bkg*(exp(akg*(ik-1))-1))>ktol) &
+         call die('pbevv_vdw_iofk ERROR: kgmesh not logarithmic')
+    end do
+    first_call = .false.
+  end if
+
+  ! Check that kf and kg are within interpolation range
+  if (kf<kfmesh(1)-ktol .or. kf>kfmesh(nkf)+ktol .or. &
+      kg<kgmesh(1)-ktol .or. kg>kgmesh(nkg)+ktol) then
+     call die('pbevv_vdw_iofk ERROR: (kf,kg) out of range')
+  endif
+
+  ! Find interpolation mesh intervals
+  ikf = 1 + int( log( 1 + (kf-kfmesh(1))/bkf ) / akf )
+  ikg = 1 + int( log( 1 + (kg-kgmesh(1))/bkg ) / akg )
+  ikf = max( 1, ikf )
+  ikg = max( 1, ikg )
+  ikf = min( nkf-1, ikf )
+  ikg = min( nkg-1, ikg )
+
+end subroutine iofk
+
+! -----------------------------------------------------------------------------
+
+subroutine kofn( n, gn, kf, kg, dkfdn, dkgdn, dkfdgn, dkgdgn )
+
+! Finds Fermi and gradient wavevectors from density and gradient
+
+  implicit none
+  real(dp), intent(in) :: n          ! Electron density
+  real(dp), intent(in) :: gn(3)      ! Density gradient
+  real(dp), intent(out):: kf         ! Local Fermi wavevector
+  real(dp), intent(out):: kg         ! |grad(n)|/n
+  real(dp), intent(out):: dkfdn      ! dkf/dn
+  real(dp), intent(out):: dkgdn      ! dkg/dn
+  real(dp), intent(out):: dkfdgn(3)  ! dkf/dgn
+  real(dp), intent(out):: dkgdgn(3)  ! dkg/dgn
+
+  real(dp):: gn2, pi
+
+! Trap exception for zero density
+  if (n <= 1.e-30_dp) then
+    kf = 0
+    kg = 0
+    dkfdn = 0
+    dkfdgn = 0
+    dkgdn = 0
+    dkgdgn = 0
+    return
+  end if
+
+! Find kf and kg
+  pi = acos(-1._dp)
+  kf = (3*pi**2*n)**(1._dp/3)
+  gn2 = sum(gn**2)
+  kg = sqrt(gn2)/n
+
+! Find derivatives
+  dkfdn = kf/n/3
+  dkfdgn = 0
+  dkgdn = -kg/n
+  dkgdgn = kg*gn/gn2
+
+end subroutine kofn
+
+!-----------------------------------------------------------------------------
+
+subroutine pofk( kf, kg, p, dpdkf, dpdkg )
+
+! Finds the values and derivatives, at (kf,kg), of the bicubic polynomials 
+! p_i(kf,kg) such that
+!    y(kf,kg) = Sum_i p_i(kf,kg) * y_i
+! is the bicubic spline interpolation at (kf,kg) of (any) function y with
+! values y_i at mesh points (kfmesh,kgmesh)_i
+
+  implicit none
+  real(dp),intent(in) :: kf, kg ! point at which the polynomials are required
+  real(dp),intent(out):: p(nkfg)      ! polynomial values at (kf,kg)
+  real(dp),intent(out):: dpdkf(nkfg)  ! dp/dkf at (kf,kg)
+  real(dp),intent(out):: dpdkg(nkfg)  ! dp/dkg at (kf,kg)
+
+  integer :: ikf, ikg, ikfg, ikf0, ikg0
+  real(dp):: a, b, dk, dkf0dkf, dkg0dkg, kf0, kg0
+  real(dp):: pkf0(nkf), dpkfdkf0(nkf), pkg0(nkg), dpkgdkg0(nkg)
+  logical, save :: first_call=.true.
+  real(dp),save :: pkf(nkf,nkf), d2pkfdkf2(nkf,nkf)
+  real(dp),save :: pkg(nkg,nkg), d2pkgdkg2(nkg,nkg)
+
+! Set up spline polynomial basis
+  if (first_call) then
+    do ikf = 1,nkf
+      pkf(:,ikf) = 0    ! ikf'th polynomial pkf(:,ikf) is one at kfmesh(ikf) 
+      pkf(ikf,ikf) = 1  ! and zero at all other points
+      call spline( kfmesh, pkf(:,ikf), nkf, 1.e30_dp, 1.e30_dp, &
+                   d2pkfdkf2(:,ikf) )
+!      call spline( kfmesh, pkf(:,ikf), nkf, 0._dp, 0._dp, d2pkfdkf2(:,ikf) )
+    end do
+    do ikg = 1,nkg
+      pkg(:,ikg) = 0
+      pkg(ikg,ikg) = 1
+      call spline( kgmesh, pkg(:,ikg), nkg, 1.e30_dp, 1.e30_dp, &
+                   d2pkgdkg2(:,ikg) )
+!      call spline( kgmesh, pkg(:,ikg), nkg, 0._dp, 0._dp, d2pkgdkg2(:,ikg) )
+    end do
+
+!   DEBUG
+    open(22,file='pkf.dat')
+    do ikf = 1,nkf
+      write(22,'(20e15.6)') kfmesh(ikf), pkf(:,ikf), d2pkfdkf2(:,ikf)
+    end do
+    close(22)
+    open(22,file='pkg.dat')
+    do ikg = 1,nkg
+      write(22,'(20e15.6)') kgmesh(ikg), pkg(:,ikg), d2pkgdkg2(:,ikg)
+    end do
+    close(22)
+!   END DEBUG
+
+    first_call = .false.
+  end if
+
+! 'Saturate' (kf,kg) values to bring them to interpolation range
+  call saturate( kf, kfcut, kf0, dkf0dkf )
+  call saturate( kg, kgcut, kg0, dkg0dkg )
+
+! Find interval of k mesh in which (kf0,kg0) point is included
+  call iofk( kf0, kg0, ikf0, ikg0 )
+
+! Evaluate pkf polynomials of spline basis at kf0
+! The splint code is in-lined here because it is a hot point
+  dk = kfmesh(ikf0+1) - kfmesh(ikf0)
+  a = (kfmesh(ikf0+1) - kf0) / dk   ! dadkf0 = -1/dk
+  b = (kf0 - kfmesh(ikf0)) / dk     ! dbdkf0 = +1/dk
+  do ikf = 1,nkf
+    pkf0(ikf) = a*pkf(ikf0,ikf) + b*pkf(ikf0+1,ikf) &
+              + ((a**3-a)*d2pkfdkf2(ikf0,ikf) +     &
+                 (b**3-b)*d2pkfdkf2(ikf0+1,ikf)) * dk**2/6
+    dpkfdkf0(ikf) = - (pkf(ikf0,ikf) - pkf(ikf0+1,ikf)) / dk &
+                    - ((3*a**2-1)*d2pkfdkf2(ikf0,ikf) -     &
+                       (3*b**2-1)*d2pkfdkf2(ikf0+1,ikf)) * dk/6
+  end do
+
+! Evaluate pkg polynomials of spline basis at kg0
+  dk = kgmesh(ikg0+1) - kgmesh(ikg0)
+  a = (kgmesh(ikg0+1) - kg0) / dk   ! dadkg0 = -1/dk
+  b = (kg0 - kgmesh(ikg0)) / dk     ! dbdkg0 = +1/dk
+  do ikg = 1,nkg
+    pkg0(ikg) = a*pkg(ikg0,ikg) + b*pkg(ikg0+1,ikg) &
+              + ((a**3-a)*d2pkgdkg2(ikg0,ikg) +     &
+                 (b**3-b)*d2pkgdkg2(ikg0+1,ikg)) * dk**2/6
+    dpkgdkg0(ikg) = - (pkg(ikg0,ikg) - pkg(ikg0+1,ikg)) / dk &
+                    - ((3*a**2-1)*d2pkgdkg2(ikg0,ikg) -     &
+                       (3*b**2-1)*d2pkgdkg2(ikg0+1,ikg)) * dk/6
+  end do
+
+! Evaluate pkf*pkg polynomials at (kf0,kg0)
+  ikfg = 0
+  do ikg = 1,nkg
+    do ikf = 1,nkf
+      ikfg = ikfg+1
+      p(ikfg) = pkf0(ikf) * pkg0(ikg)
+      dpdkf(ikfg) = dpkfdkf0(ikf)*dkf0dkf * pkg0(ikg)
+      dpdkg(ikfg) = pkf0(ikf) * dpkgdkg0(ikg)*dkg0dkg
+    end do
+  end do
+
+end subroutine pofk
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate( x, xc, y, dydx )
+
+  ! Defines a function y(x,xc) = 1/(1/x^nsat+1/xc^nsat)^(1/nsat), where nsat
+  ! is an integer set in the module header. This function is approximately
+  ! equal to x for x<xc and it saturates to xc when x->infinity
+
+  implicit none
+  real(dp),intent(in) :: x     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: y     ! Function value
+  real(dp),intent(out):: dydx  ! Derivative dy/dx
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('pbevv_vdwxc_saturate ERROR: xc<=0')
+  else if (x<0._dp) then
+    call die('pbevv_vdwxc_saturate ERROR: x<0')
+  else if (x==0._dp) then
+    y = 0
+    dydx = 1
+  else
+    z = 1/x**nsat + 1/xc**nsat
+    y = 1/z**(1._dp/nsat)
+    dydx = y/z/x**(nsat+1)
+  end if
+
+end subroutine saturate
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate_inverse( y, xc, x, dxdy )
+
+! Finds the inverse of the function defined in saturate subroutine:
+!   y=1/(1/x^n+1/xc^n)^(1/n)  =>  x=1/(1/y^n-1/xc^n)^(1/n)
+
+  implicit none
+  real(dp),intent(in) :: y     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: x     ! Inverse function value
+  real(dp),intent(out):: dxdy  ! Derivative dx/dy
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('pbevv_vdwxc_saturate_inverse ERROR: xc<=0')
+  else if (y<0._dp .or. y>=xc) then
+    call die('pbevv_vdwxc_saturate_inverse ERROR: y out of range')
+  else if (y==0._dp) then
+    x = 0
+    dxdy = 1
+  else
+    z = 1/y**nsat - 1/xc**nsat
+    x = 1/z**(1._dp/nsat)
+    dxdy = x/z/y**(nsat+1)
+  end if
+
+end subroutine saturate_inverse
+
+!-----------------------------------------------------------------------------
+
+subroutine set_kmesh()
+
+! Sets mesh of q values
+
+  implicit none
+  integer :: mkf, mkg
+
+  if (.not.kmesh_set) then
+    call set_mesh( nkf, xmax=kfcut, dxndx1=dkfmaxdkfmin )
+    call get_mesh( nkf, mkf, kfmesh )
+    call set_mesh( nkg, xmax=kgcut, dxndx1=dkgmaxdkgmin )
+    call get_mesh( nkg, mkg, kgmesh )
+    kmesh_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(/,a,/,(10f8.4))') 'pbevv_vdw_set_kmesh: kfmesh =', kfmesh
+    write(udebug,'(/,a,/,(10f8.4))') 'pbevv_vdw_set_kmesh: kgmesh =', kgmesh
+#endif /* DEBUG_XC */
+  end if
+
+end subroutine set_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine set_phi_table()
+
+! Finds and stores in memory the interpolation table (mesh points and 
+! function values) for the kernel phi(k1,k2,k).
+
+  implicit none
+  character(len=*),parameter:: myName = 'pbevv_vdwxc/set_phi_table '
+  integer :: ik, ikf1, ikf2, ikg1, ikg2, ik1, ik2, ir
+  real(dp):: dkdk0, dphidk0, dphidkmax, dphidr0, dphidrmax, &
+             k, kf1, kf2, kg1, kg2, pi, r(0:nr)
+
+! Check that table was not set yet
+  if (phi_table_set) return
+
+! Check that kf, kg, r, and k meshes have been set
+  if (.not.kmesh_set) call set_kmesh()
+  if (.not.kcut_set) call die('pbevv_vdw_set_phi_table ERROR: kcut not set')
+  forall(ir=0:nr) r(ir) = ir*dr
+  pi = acos(-1.0_dp)
+
+! Allocate arrays
+  call re_alloc( phir,     0,nr, 1,nkfg, 1,nkfg, myName//'phir' )
+  call re_alloc( phik,     0,nr, 1,nkfg, 1,nkfg, myName//'phik' )
+  call re_alloc( d2phidr2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidr2' )
+  call re_alloc( d2phidk2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidk2' )
+
+! Loop on (k1,k2) mesh points
+  do ikg2 = 1,nkg                    ! loop on kg2
+    do ikf2 = 1,nkf                  ! loop on kf2
+      ik2 = ikf2 + nkf*(ikg2-1)      ! combined (ikf2,ikg2) index
+      do ikg1 = 1,nkg                ! loop on kg1
+        do ikf1 = 1,nkf              ! loop on kf1
+          ik1 = ikf1 + nkf*(ikg1-1)  ! combined (ikf1,ikg1) index
+          if (ik1>ik2) cycle         ! since we will symmetrize at the end
+
+          ! Saturated (kf,kg) values
+          kf1 = kfmesh(ikf1)
+          kf2 = kfmesh(ikf2)
+          kg1 = kgmesh(ikg1)
+          kg2 = kgmesh(ikg2)
+
+          ! Find original (unsaturated) kf anf kg values
+          ! call saturate_inverse( kfmesh(ikf1), kfcut, kf1, dkdk0 )
+          ! call saturate_inverse( kfmesh(ikf2), kfcut, kf2, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg1), kgcut, kg1, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg2), kgcut, kg2, dkdk0 )
+
+          ! Find kernel as a function of r12
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = pbevv_vdw_phi_val( kf1, kf2, kg1, kg2, r(ir) )
+          end do
+
+          ! Kill kernel smoothly at r=rcut
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = phir(ir,ik1,ik2) * cutoff( r(ir)/rcut )
+          end do
+
+          ! Find kernel in reciprocal space
+          call radfft( 0, nr, rcut, phir(:,ik1,ik2), phik(:,ik1,ik2) )
+          phik(:,ik1,ik2) = phik(:,ik1,ik2) * (2*pi)**1.5_dp
+
+          ! Filter out above kcut
+          phik(nk:nr,ik1,ik2) = 0
+
+          ! Soft filter below kcut
+          do ik = 1,nk
+            k = ik * dk
+            phik(ik,ik1,ik2) = phik(ik,ik1,ik2) * cutoff(k/kcut)
+          end do
+
+          ! Find filtered kernel in real space
+          call radfft( 0, nr, kmax, phik(:,ik1,ik2), phir(:,ik1,ik2) )
+          phir(:,ik1,ik2) = phir(:,ik1,ik2) / (2*pi)**1.5_dp
+
+          ! Set up spline interpolation tables
+          dphidr0 = 0      ! since phi(k1,k2,r) is even in r
+          dphidk0 = 0      ! and therefore phi(k1,k2,k) is also even in k
+          dphidrmax = 0    ! since phi->0 for r->infty
+          dphidkmax = 0    ! and also when k->infty
+          call spline( dr, phir(:,ik1,ik2), nr+1, dphidr0, dphidrmax, &
+                       d2phidr2(:,ik1,ik2) )
+          call spline( dk, phik(:,ik1,ik2), nk+1, dphidk0, dphidkmax, &
+                       d2phidk2(:,ik1,ik2) )
+
+          ! Fill symmetric elements
+          phir(:,ik2,ik1) = phir(:,ik1,ik2)
+          phik(:,ik2,ik1) = phik(:,ik1,ik2)
+          d2phidr2(:,ik2,ik1) = d2phidr2(:,ik1,ik2)
+          d2phidk2(:,ik2,ik1) = d2phidk2(:,ik1,ik2)
+
+!#ifdef DEBUG_XC
+!          if (.false. .and. ik1==ik2) then
+!            print*, 'pbevv_vdw_set_kcut: ik1,ik2=', ik1, ik2
+!            call window( 0._dp, 5._dp, -1._dp, 4._dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 1._dp )
+!            call plot( nr+1, r, phi, phir(:,ik1,ik2) )
+!            call window( 0._dp, 10._dp, -0.05_dp, 0.15_dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 0.05_dp )
+!            call plot( nr+1, r, r**2*phi, r**2*phir(:,ik1,ik2))
+!            call show()
+!          end if
+!#endif /* DEBUG_XC */
+
+        end do ! ikf1
+      end do ! ikg1
+    end do ! ikf2
+  end do ! ikg2
+
+!#ifdef DEBUG_XC
+!  open(17,file='vv_phi.table')
+!  write(17,'(3a6,2a12,/,(3i6,2f15.9))') &
+!    'ik1', 'ik2', 'ir', 'phi', 'd2phi/dk2', &
+!!    (((ik1,ik2,ir,phir(ir,ik1,ik2),d2phidr2(ir,ik1,ik2), &
+!    (((ik1,ik2,ir,phik(ir,ik1,ik2),d2phidk2(ir,ik1,ik2), &
+!       ir=0,100),ik2=1,nkfg),ik1=1,nkfg)
+!  close(17)
+!#endif /* DEBUG_XC */
+
+! Mark table as set
+  phi_table_set = .true.
+
+end subroutine set_phi_table
+
+! -----------------------------------------------------------------------------
+
+real(dp) function pbevv_vdw_phi_val( kf1, kf2, kg1, kg2, r12 )
+
+! vdW energy kernel of Vydrov-vanVoorhis, eq.(2) of JCP 133, 244103 (2010)
+! This subroutine calculates the 'raw' kernel directly, without interpolarions
+! In practice, it returns n1*n2*phi(k1,k2,r12), which is smooth for kf->0
+! Input:
+!   real(dp):: kf1, kf2 ! Fermi wavevectors at points 1 and 2
+!   real(dp):: kg1, kg2 ! |grad(n)|/n at points 1 and 2
+!   real(dp):: r12      ! distance between points 1 and 2
+
+! Arguments
+  implicit none
+  real(dp),intent(in) :: kf1, kf2, kg1, kg2, r12
+
+! Internal variables
+  real(dp):: g1, g2, kappa1, kappa2, kf1m, kf2m, n1, n2, &
+             phi, pi, w01, w02, wg1, wg2, wp1, wp2
+
+! Avoid dividing by zero when kf=0 
+  kf1m = max(kf1,ktol)
+  kf2m = max(kf2,ktol)
+
+! Find kernel
+  pi = acos(-1.0_dp)
+  n1 = kf1m**3/(3*pi**2)          ! electron density at point 1
+  n2 = kf2m**3/(3*pi**2)          ! electron density at point 2
+  wp1 = sqrt(4*pi*n1)             ! local plasma frequency at point 1
+  wp2 = sqrt(4*pi*n2)             ! local plasma frequency at point 2
+  wg1 = sqrt(vv_C*kg1**4)         ! local band gap at point 1
+  wg2 = sqrt(vv_C*kg2**4)         ! local band gap at point 2
+  kappa1 = vv_b*kf1m**2/wp1       ! local VV kappa variable (eq.(9)) at point 1
+  kappa2 = vv_b*kf2m**2/wp2       ! kappa variable at point 2
+  w01 = sqrt(wg1**2+wp1**2/3)     ! local w0 frequency (eq.(5)) at point 1
+  w02 = sqrt(wg2**2+wp2**2/3)     ! local w0 frequency at point 2
+  g1 = w01*r12**2 + kappa1        ! local g variable (eq.(3)) at point 1
+  g2 = w02*r12**2 + kappa2        ! local g variable at point 2
+  phi = -1.5_dp/g1/g2/(g1+g2)     ! VV kernel phi (eq.(2))
+
+  if (kernelPrefactor=='rho') then
+    ! Return whole integrand of eq.(1)
+    pbevv_vdw_phi_val = n1*n2*phi
+  else if (kernelPrefactor=='kf') then
+    ! Return kf1*kf2*phi
+    pbevv_vdw_phi_val = kf1*kf2*phi
+  else if (kernelPrefactor=='kf2') then
+    ! Return (kf1*kf2)**2*phi
+    pbevv_vdw_phi_val = (kf1*kf2)**2*phi
+  else if (kernelPrefactor=='sqr_rho') then
+    ! Find and return sqrt(n1*n2)*phi
+    pbevv_vdw_phi_val = sqrt(n1*n2)*phi
+  else
+    call die('pbevv_vdw_phi_val ERROR: unknown kernelPrefactor')
+  end if
+
+end function pbevv_vdw_phi_val
+
+! -----------------------------------------------------------------------------
+
+real(dp) function pbevv_vdw_beta()
+
+! Returns parameter beta of VV functional
+
+  implicit none
+  pbevv_vdw_beta = vv_beta
+
+end function pbevv_vdw_beta
+
+!-----------------------------------------------------------------------------
+
+subroutine pbevv_vdw_get_kmesh( mkf, mkg, kf, kg )
+
+! Returns size and values of (kf,kg) mesh
+
+  implicit none
+  integer,          intent(out) :: mkf    ! Number of kf mesh points
+  integer,          intent(out) :: mkg    ! Number of kg mesh points
+  real(dp),optional,intent(out) :: kf(:)  ! Values of kf mesh points
+  real(dp),optional,intent(out) :: kg(:)  ! Values of kg mesh points
+  integer:: nmax
+  if (.not.kmesh_set) call set_kmesh()
+  mkf = nkf
+  mkg = nkg
+  if (present(kf)) then
+    nmax = max( nkf, size(kf) )
+    kf(1:nmax) = kfmesh(1:nmax)
+  end if
+  if (present(kg)) then
+    nmax = max( nkg, size(kg) )
+    kg(1:nmax) = kgmesh(1:nmax)
+  end if
+end subroutine pbevv_vdw_get_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine pbevv_vdw_phi( k, phi, dphidk )
+
+! Finds by interpolation phi(k1,k2,k) (Fourier transform of phi(k1,k2,r)),
+! with k1=(kf1,kg1), k2=(kf2,kg2) for all mesh values of kf (Fermi 
+! wavevector) and kg (grad(n)/n). If the interpolation table does not exist,
+! it is calculated in the first call to pbevv_vdw_phi. It requires a previous 
+! call to pbevv_vdw_set_kc to set k mesh.
+
+  implicit none
+  real(dp),intent(in) :: k            ! Modulus of actual k vector
+  real(dp),intent(out):: phi(:,:)     ! phi(k1,k2,k) at given k
+                                      ! for all k1,k2 in (kf,kg) mesh
+  real(dp),intent(out):: dphidk(:,:)  ! dphi(k1,k2,k)/dk at given k
+
+  integer :: ik, ik1, ik2
+  real(dp):: a, a2, a3, b, b2, b3
+
+! Set interpolation table
+  if (.not.phi_table_set) call set_phi_table()
+
+! Check argument sizes
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    call die('pbevv_vdw_phi: ERROR: size(phi) too small')
+
+! Find phi values at point k
+  if (k >= kcut) then
+    phi(:,:) = 0
+    dphidk(:,:) = 0
+  else
+    ! Expand interpolation inline since this is the hottest point in VdW
+    ik = int(k/dk)
+    a = ((ik+1)*dk-k)/dk
+    b = 1 - a
+    a2 = (3*a**2 -1) * dk / 6
+    b2 = (3*b**2 -1) * dk / 6
+    a3 = (a**3 - a) * dk**2 / 6
+    b3 = (b**3 - b) * dk**2 / 6
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+!        call splint( dk, phik(:,ik1,ik2), d2phidk2(:,ik1,ik2), &
+!                     nk+1, k, phi(ik1,ik2), dphidk(ik1,ik2), pr )
+        phi(ik1,ik2) = a*phik(ik,ik1,ik2) + b*phik(ik+1,ik1,ik2) &
+                + a3*d2phidk2(ik,ik1,ik2) + b3*d2phidk2(ik+1,ik1,ik2)
+        dphidk(ik1,ik2) = (-phik(ik,ik1,ik2) &
+                               +phik(ik+1,ik1,ik2) )/dk &
+                - a2*d2phidk2(ik,ik1,ik2) + b2*d2phidk2(ik+1,ik1,ik2)
+        phi(ik2,ik1) = phi(ik1,ik2)
+        dphidk(ik2,ik1) = dphidk(ik1,ik2)
+      end do
+    end do
+  end if
+
+end subroutine pbevv_vdw_phi
+
+!-----------------------------------------------------------------------------
+
+subroutine pbevv_vdw_phiofr( r, phi )
+
+! Finds phi(k1,k2,r) with k1=(kf1,kg1), k2=(kf2,kg2) for mesh values of
+! kf's (Fermi wavevectors) and kg's (grad(n)/n)
+
+  implicit none
+  real(dp),intent(in) :: r
+  real(dp),intent(out):: phi(:,:)
+
+  integer :: ikf1, ikf2, ikg1, ikg2, ik1, ik2
+  real(dp):: dphidr
+
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    stop 'vv_phiofr: ERROR: size(phi) too small'
+  if (.not.phi_table_set) call set_phi_table()
+
+  if (r >= rcut) then
+    phi(:,:) = 0
+  else
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+        call splint( dr, phir(:,ik1,ik2), d2phidr2(:,ik1,ik2), &
+                     nr+1, r, phi(ik1,ik2), dphidr )
+        phi(ik2,ik1) = phi(ik1,ik2)
+      end do ! ik1
+    end do ! ik2
+  end if ! (r>=rcut)
+
+end subroutine pbevv_vdw_phiofr
+
+!-----------------------------------------------------------------------------
+
+subroutine pbevv_vdw_set_kcut( kc )
+
+! Sets the reciprocal planewave cutoff kc of the integration grid, and finds 
+! the interpolation table to be used by vdw_phi to obtain the vdW kernel phi
+! at the reciprocal mesh wavevectors.
+
+  implicit none
+  real(dp),intent(in):: kc  ! Planewave cutoff: k>kcut => phi=0
+
+  real(dp):: pi
+
+  ! Set kcut and radial mesh parameters, if not already set
+  if (.not. kc==kcut) then
+    kcut = kc
+    pi = acos(-1._dp)
+    dr = rcut / nr
+    dk = pi / rcut
+    kmax = pi / dr
+    nk = int(kcut/dk) + 1
+    if (nk>nr) stop 'pbevv_vdw_set_kcut: ERROR: nk>nr'
+    kcut_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(a,5f8.3)') 'pbevv_vdw_set_kcut: kfcut,kgcut,rcut,kcut,kmax=', &
+      kfcut, kgcut, rcut, kc, kmax
+#endif /* DEBUG_XC */
+  end if
+
+  ! Set (kf,kg) mesh and phi table
+  call set_kmesh()
+  call set_phi_table()
+
+end subroutine pbevv_vdw_set_kcut
+
+!-----------------------------------------------------------------------------
+
+subroutine pbevv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+
+! Finds the value and derivatives of theta_i(rho,grad_rho) of eq.(8)
+! of Roman-Soler PRL 2009. In practice, theta_i=p_i rather than rho*p_i,
+! beacuse p_i is used here to expand rho1*rho2*phi rather than only phi.
+
+  implicit none
+  integer, intent(in) :: nspin                 ! Number of spin components
+  real(dp),intent(in) :: rhos(nspin)           ! Electron spin density
+  real(dp),intent(in) :: grhos(3,nspin)        ! Spin density gradient
+  real(dp),intent(out):: theta(nkfg)           ! Exp. polynomials at (kf,kg)
+  real(dp),intent(out):: dtdrho(nkfg,nspin)    ! dtheta(ik)/drhos
+  real(dp),intent(out):: dtdgrho(3,nkfg,nspin) ! dtheta(ik)/dgrhos(ix)
+
+  integer :: is, ix, ns
+  real(dp):: rho, grho(3), dpdkf(nkfg), dpdkg(nkfg), dkfdrho, dkfdgrho(3), &
+             dkgdrho, dkgdgrho(3), kf, kg, p(nkfg)
+
+  ! Sum spin components of electron density
+  ns = min(nspin,2)     ! num. of diagonal spin components (if noncollinear)
+  rho = sum(rhos(1:ns))            ! local electron density
+  grho = sum(grhos(:,1:ns),dim=2)  ! local density gradient
+
+  ! If density is between threshold, return zero
+  if (rho<rhoMin) then
+    theta = 0
+    dtdrho = 0
+    dtdgrho = 0
+    return
+  end if
+
+  ! Find local Fermi and gradient wavevectors, and their derivatives
+  call kofn( rho, grho, kf, kg, dkfdrho, dkgdrho, dkfdgrho, dkgdgrho )
+
+  ! Find expansion polynomials of integrand kernel of nonlocal vdW energy
+  call pofk( kf, kg, p, dpdkf, dpdkg )
+
+  ! Find theta functions and their derivatives with respect to rho and grho
+  if (kernelPrefactor=='rho') then
+    ! This is the right code if pbevv_vdw_phi_val returns rho1*rho2*phi
+    theta(1:nkfg) = p(1:nkfg)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = dpdkf(1:nkfg)*dkfdrho + dpdkg(1:nkfg)*dkgdrho
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix)
+      end do
+    end do
+  else if (kernelPrefactor=='kf') then
+    ! This is the right code if pbevv_vdw_phi_val returns kf1*kf2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf &
+                        + p(1:nkfg)*(1/kf-rho/kf**2*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf
+      end do
+    end do
+  else if (kernelPrefactor=='kf2') then
+    ! This is the right code if pbevv_vdw_phi_val returns (kf1*kf2)**2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf**2
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf**2 &
+                        + p(1:nkfg)*(1/kf**2-2*rho/kf**3*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf**2
+      end do
+    end do
+  else if (kernelPrefactor=='sqr_rho') then
+    ! This is the right code if pbevv_vdw_phi_val returns sqrt(rho1*rho2)*phi
+    rho = max(rho,1.e-12_dp)  ! to avoid division by zero
+    theta(1:nkfg) = p(1:nkfg) * sqrt(rho)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = ( dpdkf(1:nkfg)*dkfdrho &
+                          + dpdkg(1:nkfg)*dkgdrho ) * sqrt(rho) &
+                        + p(1:nkfg) * 0.5/sqrt(rho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) * sqrt(rho) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix) * sqrt(rho)
+      end do
+    end do
+  else
+    call die('pbevv_vdw_theta ERROR: unknown kernelPrefactor')
+  end if ! (kernelPrefactor=='rho')
+
+end subroutine pbevv_vdw_theta
+
+END MODULE m_pbevv_vdwxc
+
+!!!=============== PBEsol+VV10s =========================
+
+MODULE m_pbesvvs_vdwxc
+
+! Used module procedures
+  use sys,          only: die              ! Termination routine
+  use mesh1D,       only: get_mesh         ! Returns the mesh points
+  use m_radfft,     only: radfft           ! Radial fast Fourier transform
+  use alloc,        only: re_alloc         ! Re-allocation routine
+  use mesh1D,       only: set_mesh         ! Sets a 1D mesh
+  use interpolation,only: spline           ! Sets spline interpolation
+  use interpolation,only: splint           ! Calculates spline interpolation
+
+! Used module parameters
+  use precision,    only: dp               ! Real double precision type
+
+#ifdef DEBUG_XC
+  use debugXC,      only: udebug           ! File unit for debug output
+!  use plot_module, only: plot
+#endif /* DEBUG_XC */
+
+  implicit none
+
+! Called by m_vdwxc
+PUBLIC ::            &
+  pbesvvs_vdw_beta,       &! Returns parameter beta of the VV2010 functional
+  pbesvvs_vdw_theta,      &! Finds function theta_ik(rho,grad_rho)
+  pbesvvs_vdw_get_kmesh,  &! Returns size and values of (kf,kg) mesh
+  pbesvvs_vdw_phi,        &! Finds and interpolates rho1*rho2*phi(k1,k2,k)
+  pbesvvs_vdw_set_kcut     ! Sets the planewave cutoff kc of the integration grid
+
+#ifdef DEBUG_XC
+! Called by debugging test programs
+PUBLIC :: &
+  pbesvvs_vdw_phi_val,    &! Kernel phi(kf1,kf2,kg1,kg2,r12) with no interpolation
+  pbesvvs_vdw_phiofr       ! Kernel phi(k1,k2,r) at tabulated k1,k2-mesh values
+#endif /* DEBUG_XC */
+
+PRIVATE  ! Nothing is declared public beyond this point
+
+  ! Set parameters of the PBEsol+rVV10s,
+  ! Ref. A. V. Terentjev et al, Phys. Rev. B 98, 214108 (2018)
+  real(dp),parameter:: vv_C = 0.0093
+  real(dp),parameter:: vv_b = 10.0
+  real(dp),parameter:: vv_beta = 0.00497
+
+  ! Mesh parameters for table of phi(k1,k2,r) and its Fourier transform
+  character(len=*),parameter:: kernelPrefactor='rho' ! Prefactor of the
+                                             ! nonlocal kernel for interpolation
+                                             ! ('kf'|'kf2'|'sqr_rho'|'rho')
+  integer, parameter:: nr = 2048             ! Radial mesh points (power of 2)
+  real(dp),parameter:: rcut = 100._dp        ! Radial cutoff: r>rcut => phi=0
+  real(dp),parameter:: rmin = 1.e-6_dp       ! Min. radius as denominator
+  integer, parameter:: nkf =  7              ! Number of Fermi wavevectors
+  integer, parameter:: nkg =  5              ! Num. of grad(n)/n wavevectors
+  integer, parameter:: nkfg = nkf*nkg        ! Num. of (kf,kg) mesh points
+  real(dp),parameter:: kfcut = 5.0_dp        ! Max. Fermi wavevec.
+  real(dp),parameter:: kgcut = 5.0_dp        ! Max. grad(n)/n
+  real(dp),parameter:: dkfmaxdkfmin = 20.0_dp ! Last/first kf mesh interval
+  real(dp),parameter:: dkgmaxdkgmin = 2.0_dp  ! Last/first kg mesh interval
+  real(dp),parameter:: ktol = 1.e-12_dp      ! 'out of range' tolerance
+  real(dp),parameter:: amin = 1.e-12_dp      ! tiny denominator to avoid /0
+  real(dp),parameter:: rhoMin = 1.e-10_dp    ! Min. density for nonzero Vxc
+
+  ! Parameters for cutoff function, used in radial Fourier transforms of phi
+  integer, parameter:: ncut1 =  8      ! cutoff(x)=(1-x**ncut1)**ncut2
+  integer, parameter:: ncut2 =  4
+
+  ! Parameters for saturate function, used to enforce that q<qcut
+  integer, parameter:: nsat = 15 ! h(x,xc)=1/(1/x**nsat+1/xc**nsat)**(1/nsat)
+                                 ! nsat=15 approximately corresponds to 
+                                 ! nsat=12 for eq.(5) of Roman-Soler PRL
+
+  ! Private module variables and arrays
+  logical, save:: kcut_set=.false.         ! Has kcut been set?
+  logical, save:: kmesh_set=.false.        ! Has (kf,kg) mesh been set?
+  logical, save:: phi_table_set=.false.    ! Has phi_table been set?
+  integer, save:: nk                       ! Num. of radial mesh k-points
+  real(dp),save:: dr                       ! r-mesh interval
+  real(dp),save:: dk                       ! k-mesh interval
+  real(dp),save:: kcut                     ! Planewave cutoff: k>kcut => phi=0
+  real(dp),save:: kmax                     ! Max. k vector in Fourier transforms
+  real(dp),save:: kfmesh(nkf)              ! Mesh points for Fermi wavevector
+  real(dp),save:: kgmesh(nkg)              ! Mesh points for grad(n)/n
+  real(dp),pointer,save:: &
+                  phir(:,:,:)=>null(),    &! Table of phi(r)
+                  phik(:,:,:)=>null(),    &! Table of phi(k)
+                  d2phidr2(:,:,:)=>null(),&! Table of d^2(phi)/dr^2
+                  d2phidk2(:,:,:)=>null()  ! Table of d^2(phi)/dk^2
+
+CONTAINS
+
+! -----------------------------------------------------------------------------
+
+real(dp) function cutoff( x )
+
+! Defines a smooth cutoff function that falls from y(0)=1 to y(1)=0
+
+  implicit none
+  real(dp),intent(in):: x
+
+  if (x<=0._dp) then
+    cutoff = 1
+  else if (x>=1._dp) then
+    cutoff = 0
+  else
+    cutoff = (1-x**ncut1)**ncut2
+  end if
+
+end function cutoff
+
+!-----------------------------------------------------------------------------
+
+subroutine iofk( kf, kg, ikf, ikg )
+
+! Finds indexes ikf and ikg such that kfmesh(ikf) <= kf < kfmesh(ikf+1)
+! and kgmesh(ikg) <= kg < kgmesh(ikg+1) in logarithmic meshes of the form
+! k(ik) = b*(exp((ik-1)*a)-1)
+
+  implicit none
+  real(dp), intent(in) :: kf, kg
+  integer,  intent(out):: ikf, ikg
+
+  real(dp),save:: akf, akg, bkf, bkg
+  logical, save:: first_call = .true.
+  integer:: ik
+
+  ! Mesh data initializations
+  if (first_call) then
+    ! Find logarithmic-mesh a & b parameters: x(j)=x(1)+b*(exp(a*(j-1))-1)
+    akf = log( (kfmesh(nkf)-kfmesh(nkf-1)) / (kfmesh(2)-kfmesh(1)) ) / (nkf-2)
+    akg = log( (kgmesh(nkg)-kgmesh(nkg-1)) / (kgmesh(2)-kgmesh(1)) ) / (nkg-2)
+    akf = max( akf, amin )
+    akg = max( akg, amin )
+    bkf = (kfmesh(nkf) - kfmesh(1)) / (exp(akf*(nkf-1)) - 1)
+    bkg = (kgmesh(nkg) - kgmesh(1)) / (exp(akg*(nkg-1)) - 1)
+    ! Check that meshes are indeed logarithmic
+    do ik = 1,nkf
+      if (abs(kfmesh(ik)-kfmesh(1)-bkf*(exp(akf*(ik-1))-1))>ktol) &
+         call die('pbesvvs_vdw_iofk ERROR: kfmesh not logarithmic')
+    end do
+    do ik = 1,nkg
+      if (abs(kgmesh(ik)-kgmesh(1)-bkg*(exp(akg*(ik-1))-1))>ktol) &
+         call die('pbesvvs_vdw_iofk ERROR: kgmesh not logarithmic')
+    end do
+    first_call = .false.
+  end if
+
+  ! Check that kf and kg are within interpolation range
+  if (kf<kfmesh(1)-ktol .or. kf>kfmesh(nkf)+ktol .or. &
+      kg<kgmesh(1)-ktol .or. kg>kgmesh(nkg)+ktol) then
+     call die('pbesvvs_vdw_iofk ERROR: (kf,kg) out of range')
+  endif
+
+  ! Find interpolation mesh intervals
+  ikf = 1 + int( log( 1 + (kf-kfmesh(1))/bkf ) / akf )
+  ikg = 1 + int( log( 1 + (kg-kgmesh(1))/bkg ) / akg )
+  ikf = max( 1, ikf )
+  ikg = max( 1, ikg )
+  ikf = min( nkf-1, ikf )
+  ikg = min( nkg-1, ikg )
+
+end subroutine iofk
+
+! -----------------------------------------------------------------------------
+
+subroutine kofn( n, gn, kf, kg, dkfdn, dkgdn, dkfdgn, dkgdgn )
+
+! Finds Fermi and gradient wavevectors from density and gradient
+
+  implicit none
+  real(dp), intent(in) :: n          ! Electron density
+  real(dp), intent(in) :: gn(3)      ! Density gradient
+  real(dp), intent(out):: kf         ! Local Fermi wavevector
+  real(dp), intent(out):: kg         ! |grad(n)|/n
+  real(dp), intent(out):: dkfdn      ! dkf/dn
+  real(dp), intent(out):: dkgdn      ! dkg/dn
+  real(dp), intent(out):: dkfdgn(3)  ! dkf/dgn
+  real(dp), intent(out):: dkgdgn(3)  ! dkg/dgn
+
+  real(dp):: gn2, pi
+
+! Trap exception for zero density
+  if (n <= 1.e-30_dp) then
+    kf = 0
+    kg = 0
+    dkfdn = 0
+    dkfdgn = 0
+    dkgdn = 0
+    dkgdgn = 0
+    return
+  end if
+
+! Find kf and kg
+  pi = acos(-1._dp)
+  kf = (3*pi**2*n)**(1._dp/3)
+  gn2 = sum(gn**2)
+  kg = sqrt(gn2)/n
+
+! Find derivatives
+  dkfdn = kf/n/3
+  dkfdgn = 0
+  dkgdn = -kg/n
+  dkgdgn = kg*gn/gn2
+
+end subroutine kofn
+
+!-----------------------------------------------------------------------------
+
+subroutine pofk( kf, kg, p, dpdkf, dpdkg )
+
+! Finds the values and derivatives, at (kf,kg), of the bicubic polynomials 
+! p_i(kf,kg) such that
+!    y(kf,kg) = Sum_i p_i(kf,kg) * y_i
+! is the bicubic spline interpolation at (kf,kg) of (any) function y with
+! values y_i at mesh points (kfmesh,kgmesh)_i
+
+  implicit none
+  real(dp),intent(in) :: kf, kg ! point at which the polynomials are required
+  real(dp),intent(out):: p(nkfg)      ! polynomial values at (kf,kg)
+  real(dp),intent(out):: dpdkf(nkfg)  ! dp/dkf at (kf,kg)
+  real(dp),intent(out):: dpdkg(nkfg)  ! dp/dkg at (kf,kg)
+
+  integer :: ikf, ikg, ikfg, ikf0, ikg0
+  real(dp):: a, b, dk, dkf0dkf, dkg0dkg, kf0, kg0
+  real(dp):: pkf0(nkf), dpkfdkf0(nkf), pkg0(nkg), dpkgdkg0(nkg)
+  logical, save :: first_call=.true.
+  real(dp),save :: pkf(nkf,nkf), d2pkfdkf2(nkf,nkf)
+  real(dp),save :: pkg(nkg,nkg), d2pkgdkg2(nkg,nkg)
+
+! Set up spline polynomial basis
+  if (first_call) then
+    do ikf = 1,nkf
+      pkf(:,ikf) = 0    ! ikf'th polynomial pkf(:,ikf) is one at kfmesh(ikf) 
+      pkf(ikf,ikf) = 1  ! and zero at all other points
+      call spline( kfmesh, pkf(:,ikf), nkf, 1.e30_dp, 1.e30_dp, &
+                   d2pkfdkf2(:,ikf) )
+!      call spline( kfmesh, pkf(:,ikf), nkf, 0._dp, 0._dp, d2pkfdkf2(:,ikf) )
+    end do
+    do ikg = 1,nkg
+      pkg(:,ikg) = 0
+      pkg(ikg,ikg) = 1
+      call spline( kgmesh, pkg(:,ikg), nkg, 1.e30_dp, 1.e30_dp, &
+                   d2pkgdkg2(:,ikg) )
+!      call spline( kgmesh, pkg(:,ikg), nkg, 0._dp, 0._dp, d2pkgdkg2(:,ikg) )
+    end do
+
+!   DEBUG
+    open(22,file='pkf.dat')
+    do ikf = 1,nkf
+      write(22,'(20e15.6)') kfmesh(ikf), pkf(:,ikf), d2pkfdkf2(:,ikf)
+    end do
+    close(22)
+    open(22,file='pkg.dat')
+    do ikg = 1,nkg
+      write(22,'(20e15.6)') kgmesh(ikg), pkg(:,ikg), d2pkgdkg2(:,ikg)
+    end do
+    close(22)
+!   END DEBUG
+
+    first_call = .false.
+  end if
+
+! 'Saturate' (kf,kg) values to bring them to interpolation range
+  call saturate( kf, kfcut, kf0, dkf0dkf )
+  call saturate( kg, kgcut, kg0, dkg0dkg )
+
+! Find interval of k mesh in which (kf0,kg0) point is included
+  call iofk( kf0, kg0, ikf0, ikg0 )
+
+! Evaluate pkf polynomials of spline basis at kf0
+! The splint code is in-lined here because it is a hot point
+  dk = kfmesh(ikf0+1) - kfmesh(ikf0)
+  a = (kfmesh(ikf0+1) - kf0) / dk   ! dadkf0 = -1/dk
+  b = (kf0 - kfmesh(ikf0)) / dk     ! dbdkf0 = +1/dk
+  do ikf = 1,nkf
+    pkf0(ikf) = a*pkf(ikf0,ikf) + b*pkf(ikf0+1,ikf) &
+              + ((a**3-a)*d2pkfdkf2(ikf0,ikf) +     &
+                 (b**3-b)*d2pkfdkf2(ikf0+1,ikf)) * dk**2/6
+    dpkfdkf0(ikf) = - (pkf(ikf0,ikf) - pkf(ikf0+1,ikf)) / dk &
+                    - ((3*a**2-1)*d2pkfdkf2(ikf0,ikf) -     &
+                       (3*b**2-1)*d2pkfdkf2(ikf0+1,ikf)) * dk/6
+  end do
+
+! Evaluate pkg polynomials of spline basis at kg0
+  dk = kgmesh(ikg0+1) - kgmesh(ikg0)
+  a = (kgmesh(ikg0+1) - kg0) / dk   ! dadkg0 = -1/dk
+  b = (kg0 - kgmesh(ikg0)) / dk     ! dbdkg0 = +1/dk
+  do ikg = 1,nkg
+    pkg0(ikg) = a*pkg(ikg0,ikg) + b*pkg(ikg0+1,ikg) &
+              + ((a**3-a)*d2pkgdkg2(ikg0,ikg) +     &
+                 (b**3-b)*d2pkgdkg2(ikg0+1,ikg)) * dk**2/6
+    dpkgdkg0(ikg) = - (pkg(ikg0,ikg) - pkg(ikg0+1,ikg)) / dk &
+                    - ((3*a**2-1)*d2pkgdkg2(ikg0,ikg) -     &
+                       (3*b**2-1)*d2pkgdkg2(ikg0+1,ikg)) * dk/6
+  end do
+
+! Evaluate pkf*pkg polynomials at (kf0,kg0)
+  ikfg = 0
+  do ikg = 1,nkg
+    do ikf = 1,nkf
+      ikfg = ikfg+1
+      p(ikfg) = pkf0(ikf) * pkg0(ikg)
+      dpdkf(ikfg) = dpkfdkf0(ikf)*dkf0dkf * pkg0(ikg)
+      dpdkg(ikfg) = pkf0(ikf) * dpkgdkg0(ikg)*dkg0dkg
+    end do
+  end do
+
+end subroutine pofk
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate( x, xc, y, dydx )
+
+  ! Defines a function y(x,xc) = 1/(1/x^nsat+1/xc^nsat)^(1/nsat), where nsat
+  ! is an integer set in the module header. This function is approximately
+  ! equal to x for x<xc and it saturates to xc when x->infinity
+
+  implicit none
+  real(dp),intent(in) :: x     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: y     ! Function value
+  real(dp),intent(out):: dydx  ! Derivative dy/dx
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('pbesvvs_vdwxc_saturate ERROR: xc<=0')
+  else if (x<0._dp) then
+    call die('pbesvvs_vdwxc_saturate ERROR: x<0')
+  else if (x==0._dp) then
+    y = 0
+    dydx = 1
+  else
+    z = 1/x**nsat + 1/xc**nsat
+    y = 1/z**(1._dp/nsat)
+    dydx = y/z/x**(nsat+1)
+  end if
+
+end subroutine saturate
+
+!-----------------------------------------------------------------------------
+
+subroutine saturate_inverse( y, xc, x, dxdy )
+
+! Finds the inverse of the function defined in saturate subroutine:
+!   y=1/(1/x^n+1/xc^n)^(1/n)  =>  x=1/(1/y^n-1/xc^n)^(1/n)
+
+  implicit none
+  real(dp),intent(in) :: y     ! Independent variable
+  real(dp),intent(in) :: xc    ! Saturation value
+  real(dp),intent(out):: x     ! Inverse function value
+  real(dp),intent(out):: dxdy  ! Derivative dx/dy
+
+  real(dp):: z
+
+  if (xc<=0._dp) then
+    call die('pbesvvs_vdwxc_saturate_inverse ERROR: xc<=0')
+  else if (y<0._dp .or. y>=xc) then
+    call die('pbesvvs_vdwxc_saturate_inverse ERROR: y out of range')
+  else if (y==0._dp) then
+    x = 0
+    dxdy = 1
+  else
+    z = 1/y**nsat - 1/xc**nsat
+    x = 1/z**(1._dp/nsat)
+    dxdy = x/z/y**(nsat+1)
+  end if
+
+end subroutine saturate_inverse
+
+!-----------------------------------------------------------------------------
+
+subroutine set_kmesh()
+
+! Sets mesh of q values
+
+  implicit none
+  integer :: mkf, mkg
+
+  if (.not.kmesh_set) then
+    call set_mesh( nkf, xmax=kfcut, dxndx1=dkfmaxdkfmin )
+    call get_mesh( nkf, mkf, kfmesh )
+    call set_mesh( nkg, xmax=kgcut, dxndx1=dkgmaxdkgmin )
+    call get_mesh( nkg, mkg, kgmesh )
+    kmesh_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(/,a,/,(10f8.4))') 'pbesvvs_vdw_set_kmesh: kfmesh =', kfmesh
+    write(udebug,'(/,a,/,(10f8.4))') 'pbesvvs_vdw_set_kmesh: kgmesh =', kgmesh
+#endif /* DEBUG_XC */
+  end if
+
+end subroutine set_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine set_phi_table()
+
+! Finds and stores in memory the interpolation table (mesh points and 
+! function values) for the kernel phi(k1,k2,k).
+
+  implicit none
+  character(len=*),parameter:: myName = 'pbesvvs_vdwxc/set_phi_table '
+  integer :: ik, ikf1, ikf2, ikg1, ikg2, ik1, ik2, ir
+  real(dp):: dkdk0, dphidk0, dphidkmax, dphidr0, dphidrmax, &
+             k, kf1, kf2, kg1, kg2, pi, r(0:nr)
+
+! Check that table was not set yet
+  if (phi_table_set) return
+
+! Check that kf, kg, r, and k meshes have been set
+  if (.not.kmesh_set) call set_kmesh()
+  if (.not.kcut_set) call die('pbesvvs_vdw_set_phi_table ERROR: kcut not set')
+  forall(ir=0:nr) r(ir) = ir*dr
+  pi = acos(-1.0_dp)
+
+! Allocate arrays
+  call re_alloc( phir,     0,nr, 1,nkfg, 1,nkfg, myName//'phir' )
+  call re_alloc( phik,     0,nr, 1,nkfg, 1,nkfg, myName//'phik' )
+  call re_alloc( d2phidr2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidr2' )
+  call re_alloc( d2phidk2, 0,nr, 1,nkfg, 1,nkfg, myName//'d2phidk2' )
+
+! Loop on (k1,k2) mesh points
+  do ikg2 = 1,nkg                    ! loop on kg2
+    do ikf2 = 1,nkf                  ! loop on kf2
+      ik2 = ikf2 + nkf*(ikg2-1)      ! combined (ikf2,ikg2) index
+      do ikg1 = 1,nkg                ! loop on kg1
+        do ikf1 = 1,nkf              ! loop on kf1
+          ik1 = ikf1 + nkf*(ikg1-1)  ! combined (ikf1,ikg1) index
+          if (ik1>ik2) cycle         ! since we will symmetrize at the end
+
+          ! Saturated (kf,kg) values
+          kf1 = kfmesh(ikf1)
+          kf2 = kfmesh(ikf2)
+          kg1 = kgmesh(ikg1)
+          kg2 = kgmesh(ikg2)
+
+          ! Find original (unsaturated) kf anf kg values
+          ! call saturate_inverse( kfmesh(ikf1), kfcut, kf1, dkdk0 )
+          ! call saturate_inverse( kfmesh(ikf2), kfcut, kf2, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg1), kgcut, kg1, dkdk0 )
+          ! call saturate_inverse( kgmesh(ikg2), kgcut, kg2, dkdk0 )
+
+          ! Find kernel as a function of r12
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = pbesvvs_vdw_phi_val( kf1, kf2, kg1, kg2, r(ir) )
+          end do
+
+          ! Kill kernel smoothly at r=rcut
+          do ir = 0,nr
+            phir(ir,ik1,ik2) = phir(ir,ik1,ik2) * cutoff( r(ir)/rcut )
+          end do
+
+          ! Find kernel in reciprocal space
+          call radfft( 0, nr, rcut, phir(:,ik1,ik2), phik(:,ik1,ik2) )
+          phik(:,ik1,ik2) = phik(:,ik1,ik2) * (2*pi)**1.5_dp
+
+          ! Filter out above kcut
+          phik(nk:nr,ik1,ik2) = 0
+
+          ! Soft filter below kcut
+          do ik = 1,nk
+            k = ik * dk
+            phik(ik,ik1,ik2) = phik(ik,ik1,ik2) * cutoff(k/kcut)
+          end do
+
+          ! Find filtered kernel in real space
+          call radfft( 0, nr, kmax, phik(:,ik1,ik2), phir(:,ik1,ik2) )
+          phir(:,ik1,ik2) = phir(:,ik1,ik2) / (2*pi)**1.5_dp
+
+          ! Set up spline interpolation tables
+          dphidr0 = 0      ! since phi(k1,k2,r) is even in r
+          dphidk0 = 0      ! and therefore phi(k1,k2,k) is also even in k
+          dphidrmax = 0    ! since phi->0 for r->infty
+          dphidkmax = 0    ! and also when k->infty
+          call spline( dr, phir(:,ik1,ik2), nr+1, dphidr0, dphidrmax, &
+                       d2phidr2(:,ik1,ik2) )
+          call spline( dk, phik(:,ik1,ik2), nk+1, dphidk0, dphidkmax, &
+                       d2phidk2(:,ik1,ik2) )
+
+          ! Fill symmetric elements
+          phir(:,ik2,ik1) = phir(:,ik1,ik2)
+          phik(:,ik2,ik1) = phik(:,ik1,ik2)
+          d2phidr2(:,ik2,ik1) = d2phidr2(:,ik1,ik2)
+          d2phidk2(:,ik2,ik1) = d2phidk2(:,ik1,ik2)
+
+!#ifdef DEBUG_XC
+!          if (.false. .and. ik1==ik2) then
+!            print*, 'pbesvvs_vdw_set_kcut: ik1,ik2=', ik1, ik2
+!            call window( 0._dp, 5._dp, -1._dp, 4._dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 1._dp )
+!            call plot( nr+1, r, phi, phir(:,ik1,ik2) )
+!            call window( 0._dp, 10._dp, -0.05_dp, 0.15_dp, 0 )
+!            call axes( 0._dp, 1._dp, 0._dp, 0.05_dp )
+!            call plot( nr+1, r, r**2*phi, r**2*phir(:,ik1,ik2))
+!            call show()
+!          end if
+!#endif /* DEBUG_XC */
+
+        end do ! ikf1
+      end do ! ikg1
+    end do ! ikf2
+  end do ! ikg2
+
+!#ifdef DEBUG_XC
+!  open(17,file='vv_phi.table')
+!  write(17,'(3a6,2a12,/,(3i6,2f15.9))') &
+!    'ik1', 'ik2', 'ir', 'phi', 'd2phi/dk2', &
+!!    (((ik1,ik2,ir,phir(ir,ik1,ik2),d2phidr2(ir,ik1,ik2), &
+!    (((ik1,ik2,ir,phik(ir,ik1,ik2),d2phidk2(ir,ik1,ik2), &
+!       ir=0,100),ik2=1,nkfg),ik1=1,nkfg)
+!  close(17)
+!#endif /* DEBUG_XC */
+
+! Mark table as set
+  phi_table_set = .true.
+
+end subroutine set_phi_table
+
+! -----------------------------------------------------------------------------
+
+real(dp) function pbesvvs_vdw_phi_val( kf1, kf2, kg1, kg2, r12 )
+
+! vdW energy kernel of Vydrov-vanVoorhis, eq.(2) of JCP 133, 244103 (2010)
+! This subroutine calculates the 'raw' kernel directly, without interpolarions
+! In practice, it returns n1*n2*phi(k1,k2,r12), which is smooth for kf->0
+! Input:
+!   real(dp):: kf1, kf2 ! Fermi wavevectors at points 1 and 2
+!   real(dp):: kg1, kg2 ! |grad(n)|/n at points 1 and 2
+!   real(dp):: r12      ! distance between points 1 and 2
+
+! Arguments
+  implicit none
+  real(dp),intent(in) :: kf1, kf2, kg1, kg2, r12
+
+! Internal variables
+  real(dp):: g1, g2, kappa1, kappa2, kf1m, kf2m, n1, n2, &
+             phi, pi, w01, w02, wg1, wg2, wp1, wp2, &
+             ss1, ss2, vv_C1, vv_C2, vv_Ct1, vv_Ct2
+
+!! New parameters of the PBEsol+rVV10s,
+!! Ref. A.V. Terentjev et al, Phys. Rev. B 98, 214108 (2018) 
+  vv_C1 = 0.5_dp
+  vv_C2 = 300.0_dp 
+
+! Avoid dividing by zero when kf=0 
+  kf1m = max(kf1,ktol)
+  kf2m = max(kf2,ktol)
+
+! Find kernel
+  pi = acos(-1.0_dp)
+  n1 = kf1m**3/(3*pi**2)          ! electron density at point 1
+  n2 = kf2m**3/(3*pi**2)          ! electron density at point 2
+!! New formulas of the PBEsol+rVV10s,
+!! Ref. A.V. Terentjev et al, Phys. Rev. B 98, 214108 (2018) 
+  ss1 = 0.1616204568582_dp*kg1/(n1**(1/3))
+  ss2 = 0.1616204568582_dp*kg2/(n2**(1/3))
+  vv_Ct1 = vv_C+vv_C1/(1.0_dp + vv_C2*(ss1-0.5_dp)**2)
+  vv_Ct2 = vv_C+vv_C1/(1.0_dp + vv_C2*(ss2-0.5_dp)**2)
+!! 
+  wp1 = sqrt(4*pi*n1)             ! local plasma frequency at point 1
+  wp2 = sqrt(4*pi*n2)             ! local plasma frequency at point 2
+!!  
+  wg1 = sqrt(vv_Ct1*kg1**4)         ! local band gap at point 1
+  wg2 = sqrt(vv_Ct2*kg2**4)         ! local band gap at point 2
+!!  
+  kappa1 = vv_b*kf1m**2/wp1       ! local VV kappa variable (eq.(9)) at point 1
+  kappa2 = vv_b*kf2m**2/wp2       ! kappa variable at point 2
+  w01 = sqrt(wg1**2+wp1**2/3)     ! local w0 frequency (eq.(5)) at point 1
+  w02 = sqrt(wg2**2+wp2**2/3)     ! local w0 frequency at point 2
+  g1 = w01*r12**2 + kappa1        ! local g variable (eq.(3)) at point 1
+  g2 = w02*r12**2 + kappa2        ! local g variable at point 2
+  phi = -1.5_dp/g1/g2/(g1+g2)     ! VV kernel phi (eq.(2))
+
+  if (kernelPrefactor=='rho') then
+    ! Return whole integrand of eq.(1)
+    pbesvvs_vdw_phi_val = n1*n2*phi
+  else if (kernelPrefactor=='kf') then
+    ! Return kf1*kf2*phi
+    pbesvvs_vdw_phi_val = kf1*kf2*phi
+  else if (kernelPrefactor=='kf2') then
+    ! Return (kf1*kf2)**2*phi
+    pbesvvs_vdw_phi_val = (kf1*kf2)**2*phi
+  else if (kernelPrefactor=='sqr_rho') then
+    ! Find and return sqrt(n1*n2)*phi
+    pbesvvs_vdw_phi_val = sqrt(n1*n2)*phi
+  else
+    call die('pbesvvs_vdw_phi_val ERROR: unknown kernelPrefactor')
+  end if
+
+end function pbesvvs_vdw_phi_val
+
+! -----------------------------------------------------------------------------
+
+real(dp) function pbesvvs_vdw_beta()
+
+! Returns parameter beta of VV functional
+
+  implicit none
+  pbesvvs_vdw_beta = vv_beta
+
+end function pbesvvs_vdw_beta
+
+!-----------------------------------------------------------------------------
+
+subroutine pbesvvs_vdw_get_kmesh( mkf, mkg, kf, kg )
+
+! Returns size and values of (kf,kg) mesh
+
+  implicit none
+  integer,          intent(out) :: mkf    ! Number of kf mesh points
+  integer,          intent(out) :: mkg    ! Number of kg mesh points
+  real(dp),optional,intent(out) :: kf(:)  ! Values of kf mesh points
+  real(dp),optional,intent(out) :: kg(:)  ! Values of kg mesh points
+  integer:: nmax
+  if (.not.kmesh_set) call set_kmesh()
+  mkf = nkf
+  mkg = nkg
+  if (present(kf)) then
+    nmax = max( nkf, size(kf) )
+    kf(1:nmax) = kfmesh(1:nmax)
+  end if
+  if (present(kg)) then
+    nmax = max( nkg, size(kg) )
+    kg(1:nmax) = kgmesh(1:nmax)
+  end if
+end subroutine pbesvvs_vdw_get_kmesh
+
+! -----------------------------------------------------------------------------
+
+subroutine pbesvvs_vdw_phi( k, phi, dphidk )
+
+! Finds by interpolation phi(k1,k2,k) (Fourier transform of phi(k1,k2,r)),
+! with k1=(kf1,kg1), k2=(kf2,kg2) for all mesh values of kf (Fermi 
+! wavevector) and kg (grad(n)/n). If the interpolation table does not exist,
+! it is calculated in the first call to pbesvvs_vdw_phi. It requires a previous 
+! call to pbesvvs_vdw_set_kc to set k mesh.
+
+  implicit none
+  real(dp),intent(in) :: k            ! Modulus of actual k vector
+  real(dp),intent(out):: phi(:,:)     ! phi(k1,k2,k) at given k
+                                      ! for all k1,k2 in (kf,kg) mesh
+  real(dp),intent(out):: dphidk(:,:)  ! dphi(k1,k2,k)/dk at given k
+
+  integer :: ik, ik1, ik2
+  real(dp):: a, a2, a3, b, b2, b3
+
+! Set interpolation table
+  if (.not.phi_table_set) call set_phi_table()
+
+! Check argument sizes
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    call die('pbesvvs_vdw_phi: ERROR: size(phi) too small')
+
+! Find phi values at point k
+  if (k >= kcut) then
+    phi(:,:) = 0
+    dphidk(:,:) = 0
+  else
+    ! Expand interpolation inline since this is the hottest point in VdW
+    ik = int(k/dk)
+    a = ((ik+1)*dk-k)/dk
+    b = 1 - a
+    a2 = (3*a**2 -1) * dk / 6
+    b2 = (3*b**2 -1) * dk / 6
+    a3 = (a**3 - a) * dk**2 / 6
+    b3 = (b**3 - b) * dk**2 / 6
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+!        call splint( dk, phik(:,ik1,ik2), d2phidk2(:,ik1,ik2), &
+!                     nk+1, k, phi(ik1,ik2), dphidk(ik1,ik2), pr )
+        phi(ik1,ik2) = a*phik(ik,ik1,ik2) + b*phik(ik+1,ik1,ik2) &
+                + a3*d2phidk2(ik,ik1,ik2) + b3*d2phidk2(ik+1,ik1,ik2)
+        dphidk(ik1,ik2) = (-phik(ik,ik1,ik2) &
+                               +phik(ik+1,ik1,ik2) )/dk &
+                - a2*d2phidk2(ik,ik1,ik2) + b2*d2phidk2(ik+1,ik1,ik2)
+        phi(ik2,ik1) = phi(ik1,ik2)
+        dphidk(ik2,ik1) = dphidk(ik1,ik2)
+      end do
+    end do
+  end if
+
+end subroutine pbesvvs_vdw_phi
+
+!-----------------------------------------------------------------------------
+
+subroutine pbesvvs_vdw_phiofr( r, phi )
+
+! Finds phi(k1,k2,r) with k1=(kf1,kg1), k2=(kf2,kg2) for mesh values of
+! kf's (Fermi wavevectors) and kg's (grad(n)/n)
+
+  implicit none
+  real(dp),intent(in) :: r
+  real(dp),intent(out):: phi(:,:)
+
+  integer :: ikf1, ikf2, ikg1, ikg2, ik1, ik2
+  real(dp):: dphidr
+
+  if (size(phi,1)<nkfg .or. size(phi,2)<nkfg) &
+    stop 'vv_phiofr: ERROR: size(phi) too small'
+  if (.not.phi_table_set) call set_phi_table()
+
+  if (r >= rcut) then
+    phi(:,:) = 0
+  else
+    do ik2 = 1,nkfg
+      do ik1 = 1,ik2
+        call splint( dr, phir(:,ik1,ik2), d2phidr2(:,ik1,ik2), &
+                     nr+1, r, phi(ik1,ik2), dphidr )
+        phi(ik2,ik1) = phi(ik1,ik2)
+      end do ! ik1
+    end do ! ik2
+  end if ! (r>=rcut)
+
+end subroutine pbesvvs_vdw_phiofr
+
+!-----------------------------------------------------------------------------
+
+subroutine pbesvvs_vdw_set_kcut( kc )
+
+! Sets the reciprocal planewave cutoff kc of the integration grid, and finds 
+! the interpolation table to be used by vdw_phi to obtain the vdW kernel phi
+! at the reciprocal mesh wavevectors.
+
+  implicit none
+  real(dp),intent(in):: kc  ! Planewave cutoff: k>kcut => phi=0
+
+  real(dp):: pi
+
+  ! Set kcut and radial mesh parameters, if not already set
+  if (.not. kc==kcut) then
+    kcut = kc
+    pi = acos(-1._dp)
+    dr = rcut / nr
+    dk = pi / rcut
+    kmax = pi / dr
+    nk = int(kcut/dk) + 1
+    if (nk>nr) stop 'pbesvvs_vdw_set_kcut: ERROR: nk>nr'
+    kcut_set = .true.
+#ifdef DEBUG_XC
+    write(udebug,'(a,5f8.3)') 'pbesvvs_vdw_set_kcut: kfcut,kgcut,rcut,kcut,kmax=', &
+      kfcut, kgcut, rcut, kc, kmax
+#endif /* DEBUG_XC */
+  end if
+
+  ! Set (kf,kg) mesh and phi table
+  call set_kmesh()
+  call set_phi_table()
+
+end subroutine pbesvvs_vdw_set_kcut
+
+!-----------------------------------------------------------------------------
+
+subroutine pbesvvs_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho )
+
+! Finds the value and derivatives of theta_i(rho,grad_rho) of eq.(8)
+! of Roman-Soler PRL 2009. In practice, theta_i=p_i rather than rho*p_i,
+! beacuse p_i is used here to expand rho1*rho2*phi rather than only phi.
+
+  implicit none
+  integer, intent(in) :: nspin                 ! Number of spin components
+  real(dp),intent(in) :: rhos(nspin)           ! Electron spin density
+  real(dp),intent(in) :: grhos(3,nspin)        ! Spin density gradient
+  real(dp),intent(out):: theta(nkfg)           ! Exp. polynomials at (kf,kg)
+  real(dp),intent(out):: dtdrho(nkfg,nspin)    ! dtheta(ik)/drhos
+  real(dp),intent(out):: dtdgrho(3,nkfg,nspin) ! dtheta(ik)/dgrhos(ix)
+
+  integer :: is, ix, ns
+  real(dp):: rho, grho(3), dpdkf(nkfg), dpdkg(nkfg), dkfdrho, dkfdgrho(3), &
+             dkgdrho, dkgdgrho(3), kf, kg, p(nkfg)
+
+  ! Sum spin components of electron density
+  ns = min(nspin,2)     ! num. of diagonal spin components (if noncollinear)
+  rho = sum(rhos(1:ns))            ! local electron density
+  grho = sum(grhos(:,1:ns),dim=2)  ! local density gradient
+
+  ! If density is between threshold, return zero
+  if (rho<rhoMin) then
+    theta = 0
+    dtdrho = 0
+    dtdgrho = 0
+    return
+  end if
+
+  ! Find local Fermi and gradient wavevectors, and their derivatives
+  call kofn( rho, grho, kf, kg, dkfdrho, dkgdrho, dkfdgrho, dkgdgrho )
+
+  ! Find expansion polynomials of integrand kernel of nonlocal vdW energy
+  call pofk( kf, kg, p, dpdkf, dpdkg )
+
+  ! Find theta functions and their derivatives with respect to rho and grho
+  if (kernelPrefactor=='rho') then
+    ! This is the right code if pbesvvs_vdw_phi_val returns rho1*rho2*phi
+    theta(1:nkfg) = p(1:nkfg)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = dpdkf(1:nkfg)*dkfdrho + dpdkg(1:nkfg)*dkgdrho
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix)
+      end do
+    end do
+  else if (kernelPrefactor=='kf') then
+    ! This is the right code if pbesvvs_vdw_phi_val returns kf1*kf2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf &
+                        + p(1:nkfg)*(1/kf-rho/kf**2*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf
+      end do
+    end do
+  else if (kernelPrefactor=='kf2') then
+    ! This is the right code if pbesvvs_vdw_phi_val returns (kf1*kf2)**2*phi
+    kf = max(kf,ktol)
+    theta(1:nkfg) = p(1:nkfg)*rho/kf**2
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = (dpdkf(1:nkfg)*dkfdrho &
+                          +dpdkg(1:nkfg)*dkgdrho)*rho/kf**2 &
+                        + p(1:nkfg)*(1/kf**2-2*rho/kf**3*dkfdrho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = (dpdkf(1:nkfg)*dkfdgrho(ix) &
+                                +dpdkg(1:nkfg)*dkgdgrho(ix))*rho/kf**2
+      end do
+    end do
+  else if (kernelPrefactor=='sqr_rho') then
+    ! This is the right code if pbesvvs_vdw_phi_val returns sqrt(rho1*rho2)*phi
+    rho = max(rho,1.e-12_dp)  ! to avoid division by zero
+    theta(1:nkfg) = p(1:nkfg) * sqrt(rho)
+    do is = 1,ns
+      dtdrho(1:nkfg,is) = ( dpdkf(1:nkfg)*dkfdrho &
+                          + dpdkg(1:nkfg)*dkgdrho ) * sqrt(rho) &
+                        + p(1:nkfg) * 0.5/sqrt(rho)
+      do ix = 1,3
+        dtdgrho(ix,1:nkfg,is) = dpdkf(1:nkfg)*dkfdgrho(ix) * sqrt(rho) &
+                              + dpdkg(1:nkfg)*dkgdgrho(ix) * sqrt(rho)
+      end do
+    end do
+  else
+    call die('pbesvvs_vdw_theta ERROR: unknown kernelPrefactor')
+  end if ! (kernelPrefactor=='rho')
+
+end subroutine pbesvvs_vdw_theta
+
+END MODULE m_pbesvvs_vdwxc
+
+
 
=== modified file 'Src/atom.F'
--- Src/atom.F	2019-06-20 10:59:42 +0000
+++ Src/atom.F	2019-09-06 14:47:07 +0000
@@ -2038,6 +2038,8 @@
         if (icorr .eq. "wc") ps_string ="GGA WC"
         if (icorr .eq. "bl") ps_string ="GGA BLYP"
         if (icorr .eq. "ps") ps_string ="GGA PBEsol"
+!
+        if (icorr .eq. "sg") ps_string ="GGA SG4"        
 
         if (icorr .eq. "jo") ps_string ="GGA PBEJsJrLO"
         if (icorr .eq. "jh") ps_string ="GGA PBEJsJrHEG"
@@ -2131,8 +2133,17 @@
      .          write(6,'(a,1x,2a)')
      .          'xc_check: WARNING: Pseudopotential generated with',
      .           trim(ps_string), " functional"
-
-
+!
+          elseif((XCauth(nf).eq.'SG4').and.(XCfunc(nf).eq.'GGA')) 
+     .        then
+
+            write(6,'(a)')
+     .       'xc_check: GGA SG4'
+            if (icorr.ne.'sg'.and.nXCfunc.eq.1) 
+     .          write(6,'(a,1x,2a)')
+     .          'xc_check: WARNING: Pseudopotential generated with',
+     .           trim(ps_string), " functional"
+            
           elseif((XCauth(nf).eq.'PW91').and.(XCfunc(nf).eq.'GGA')) then
 
             write(6,'(a)')
 
=== modified file 'Tests/Makefile'
--- Tests/Makefile	2019-01-30 15:22:16 +0000
+++ Tests/Makefile	2019-09-06 14:47:07 +0000
@@ -42,13 +42,14 @@
 
 label = work
 
-tests = 32_h2o ag anneal-cont ar2_vdw batio3 benzene bessel bessel-rich \
+tests = 32_h2o ag ag_sg4 anneal-cont ar2_vdw batio3 benzene bessel bessel-rich \
 	born born_spin carbon_nanoscroll ch4 chargeconf-h2o \
 	constant_volume dipole_correction fe fe_broyden \
 	fe_clust_noncollinear fe_clust_noncollinear-gga fe_cohp \
 	fen fe_noncol_kp fire_benzene floating \
 	force_2 force_constants gate_G_charge gate_G_hartree ge111 \
-	ge_fatbands_so	graphite_c6 graphite_c6_full graphite_vdw_df \
+	ge_fatbands_so	graphite_c6 graphite_c6_full graphite_c6_sg4 graphite_vdw_df \
+        graphite_vdw_pbesvvs graphite_vdw_pbevv graphite_vdw_sg4vv \
 	h2_bessel h2o h2o_2 h2o_am05 h2o_bands h2o_bands_nc \
 	h2o_bands_polarized h2o_basis h2o_coop h2o_dipole h2o_dipole2 \
 	h2o_dos h2o_filteret_basis h2o_findp_bug h2o_netcdf h2o_op_broyden \
 
=== modified file 'version.info'
--- version.info	2019-09-04 13:09:45 +0000
+++ version.info	2019-09-06 14:47:07 +0000
@@ -1,1 +1,1 @@
-trunk-781
+trunk-782