! Copyright (C) 2001-2009 Quantum ESPRESSO group
! Copyright (C) 2019 Brian Kolb, Timo Thonhauser
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
! ----------------------------------------------------------------------


MODULE vdW_DF

!! This module calculates the non-local correlation contribution to the
!! energy and potential according to:
!
!! M. Dion, H. Rydberg, E. Schroeder, D.C. Langreth, and
!! B.I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
!
!! henceforth referred to as DION. Further information about the
!! functional and its corresponding potential can be found in:
!
!! T. Thonhauser, V.R. Cooper, S. Li, A. Puzder, P. Hyldgaard,
!! and D.C. Langreth, Phys. Rev. B 76, 125112 (2007).
!
!! The proper spin extension of vdW-DF, i.e. svdW-DF, is derived in:
!
!! T. Thonhauser, S. Zuluaga, C.A. Arter, K. Berland, E. Schroder,
!! and P. Hyldgaard, Phys. Rev. Lett. 115, 136402 (2015).
!
!! henceforth referred to as THONHAUSER.
!
!! Two review articles show many of the vdW-DF applications:
!
!! D.C. Langreth et al., J. Phys.: Condens. Matter 21, 084203 (2009).
!
!! K. Berland et al., Rep. Prog. Phys. 78, 066501 (2015).
!
!! The method implemented is based on the method of G. Roman-Perez and
!! J.M. Soler described in:
!
!! G. Roman-Perez and J.M. Soler, Phys. Rev. Lett. 103, 096102 (2009).
!
!! henceforth referred to as SOLER.
!
!
! xc_vdW_DF and xc_vdW_DF_spin are the driver routines for vdW-DF
! calculations and are called from Modules/funct.f90. The routines here
! set up the parallel run (if any) and carry out the calls necessary to
! calculate the non-local correlation contributions to the energy and
! potential.


USE kinds,             ONLY : dp
USE constants,         ONLY : pi, fpi, e2
USE mp,                ONLY : mp_sum, mp_barrier, mp_get, mp_size, mp_rank, mp_bcast
USE mp_images,         ONLY : intra_image_comm
USE mp_bands,          ONLY : intra_bgrp_comm
USE io_global,         ONLY : stdout, ionode
USE fft_base,          ONLY : dfftp
USE fft_interfaces,    ONLY : fwfft, invfft
USE control_flags,     ONLY : iverbosity, gamma_only

! ----------------------------------------------------------------------
! No implicit variables

IMPLICIT NONE


! ----------------------------------------------------------------------
! By default everything is private

PRIVATE


! ----------------------------------------------------------------------
! Save all objects in this module

SAVE


! ----------------------------------------------------------------------
! Public functions

PUBLIC  :: xc_vdW_DF, xc_vdW_DF_spin, vdW_DF_stress,                   &
           vdW_DF_energy, vdW_DF_potential,                            &
           generate_kernel, interpolate_kernel,                        &
           initialize_spline_interpolation, spline_interpolation, pw,pw_spin


! ----------------------------------------------------------------------
! Public variables

PUBLIC  :: inlc, vdW_DF_analysis, Nr_points, r_max, q_min, q_cut, Nqs, q_mesh


! ----------------------------------------------------------------------
! General variables

INTEGER                  :: inlc            = 1
! The non-local correlation

INTEGER                  :: vdW_DF_analysis = 0
! vdW-DF analysis tool as described in PRB 97, 085115 (2018)

REAL(DP), PARAMETER      :: epsr            = 1.0D-12
! A small number to cut off densities

INTEGER                  :: idx
! Indexing variable


! ----------------------------------------------------------------------
! Kernel specific parameters and variables

INTEGER, PARAMETER       :: Nr_points = 1024
! The number of radial points (also the number of k points) used in the
! formation of the kernel functions for each pair of q values.
! Increasing this value will help in case you get a run-time error
! saying that you are trying to use a k value that is larger than the
! largest tabulated k point since the largest k point will be 2*pi/r_max
! * Nr_points. Memory usage of the vdW_DF piece of PWSCF will increase
! roughly linearly with this variable.

REAL(DP), PARAMETER      :: r_max     = 100.0D0
! The value of the maximum radius to use for the real-space kernel
! functions for each pair of q values. The larger this value is the
! smaller the smallest k value will be since the smallest k point value
! is 2*pi/r_max. Be careful though, since this will also decrease the
! maximum k point value and the vdW_DF code will crash if it encounters
! a g-vector with a magnitude greater than 2*pi/r_max *Nr_points.

REAL(DP), PARAMETER      :: dr = r_max/Nr_points, dk = 2.0D0*pi/r_max
! Real space and k-space spacing of grid points.

REAL(DP), PARAMETER      :: q_min = 1.0D-5, q_cut = 5.0D0
! The maximum and minimum values of q. During a vdW run, values of q0
! found larger than q_cut will be saturated (SOLER equation 5) to q_cut.

INTEGER,  PARAMETER                 :: Nqs    = 20
REAL(DP), PARAMETER, DIMENSION(Nqs) :: q_mesh = (/  &
   q_min              , 0.0449420825586261D0, 0.0975593700991365D0, 0.159162633466142D0, &
   0.231286496836006D0, 0.315727667369529D0 , 0.414589693721418D0 , 0.530335368404141D0, &
   0.665848079422965D0, 0.824503639537924D0 , 1.010254382520950D0 , 1.227727621364570D0, &
   1.482340921174910D0, 1.780437058359530D0 , 2.129442028133640D0 , 2.538050036534580D0, &
   3.016440085356680D0, 3.576529545442460D0 , 4.232271035198720D0 , q_cut /)

! The above two parameters define the q mesh to be used in the vdW_DF
! code. These are perhaps the most important to have set correctly.
! Increasing the number of q points will DRAMATICALLY increase the
! memory usage of the vdW_DF code because the memory consumption depends
! quadratically on the number of q points in the mesh. Increasing the
! number of q points may increase accuracy of the vdW_DF code, although,
! in testing it was found to have little effect. The largest value of
! the q mesh is q_cut. All values of q0 (DION equation 11) larger than
! this value during a run will be saturated to this value using equation
! 5 of SOLER. In testing, increasing the value of q_cut was found to
! have little impact on the results, although it is possible that in
! some systems it may be more important. Always make sure that the
! variable Nqs is consistent with the number of q points that are
! actually in the variable q_mesh. Also, do not set any q value to 0.
! This will cause an infinity in the Fourier transform.

INTEGER,  PARAMETER      :: Nintegration_points = 256
! Number of integration points for real-space kernel generation (see
! DION equation 14). This is how many a's and b's there will be.

REAL(DP), PARAMETER      :: a_min = 0.0D0, a_max = 64.0D0
! Min/max values for the a and b integration in DION equation 14.

REAL(DP) :: kernel( 0:Nr_points, Nqs, Nqs ), d2phi_dk2( 0:Nr_points, Nqs, Nqs )
! Matrices holding the Fourier transformed kernel function  and its
! second derivative for each pair of q values. The ordering is
! kernel(k_point, q1_value, q2_value).

REAL(DP) :: W_ab( Nintegration_points, Nintegration_points )
! Defined in DION equation 16.

REAL(DP) :: a_points( Nintegration_points ), a_points2( Nintegration_points )
! The values of the "a" points (DION equation 14) and their squares.

CONTAINS








  ! ####################################################################
  !                           |             |
  !                           |  functions  |
  !                           |_____________|
  !
  ! Functions to be used in get_q0_on_grid, get_q0_on_grid_spin, and
  ! phi_value.

  FUNCTION Fs(s)

     IMPLICIT NONE
     REAL(DP) :: s, Fs, Z_ab = 0.0D0

     IF ( inlc == 1 .OR. inlc == 3 ) THEN
        Z_ab = -0.8491D0
     ELSE IF ( inlc == 2 .OR. inlc == 4 .OR. inlc == 5 ) THEN
        Z_ab = -1.887D0
     END IF

     Fs = 1.0D0 - Z_ab * s * s / 9.0D0

  END FUNCTION Fs




  FUNCTION dFs_ds(s)

     IMPLICIT NONE
     REAL(DP)             :: s, dFs_ds, Z_ab = 0.0D0
     REAL(DP), PARAMETER  :: prefac = -2.0D0/9.0D0

     IF ( inlc == 1 .OR. inlc == 3 ) THEN
        Z_ab = -0.8491D0
     ELSE IF ( inlc == 2 .OR. inlc == 4 .OR. inlc == 5 ) THEN
        Z_ab = -1.887D0
     END IF

     dFs_ds =  prefac * s * Z_ab

  END FUNCTION dFs_ds




  FUNCTION kF(rho)

     IMPLICIT NONE
     REAL(DP)             :: rho, kF
     REAL(DP), PARAMETER  :: ex = 1.0D0/3.0D0

     kF = ( 3.0D0 * pi * pi * rho )**ex

  END FUNCTION kF




  FUNCTION dkF_drho(rho)

     IMPLICIT NONE
     REAL(DP)             :: rho, dkF_drho
     REAL(DP), PARAMETER  :: prefac = 1.0D0/3.0D0

     dkF_drho = prefac * kF(rho) / rho

  END FUNCTION dkF_drho




  FUNCTION ds_drho(rho, s)

     IMPLICIT NONE
     REAL(DP) :: rho, s, ds_drho

     ds_drho = -s * ( dkF_drho(rho) / kF(rho) + 1.0D0 / rho )

  END FUNCTION ds_drho




  FUNCTION ds_dgradrho(rho)

     IMPLICIT NONE
     REAL(DP) :: rho, ds_dgradrho

     ds_dgradrho = 0.5D0 / (kF(rho) * rho)

  END FUNCTION ds_dgradrho




  FUNCTION dqx_drho(rho, s)

     IMPLICIT NONE
     REAL(DP) :: rho, s, dqx_drho

     dqx_drho = dkF_drho(rho) * Fs(s) + kF(rho) * dFs_ds(s) * ds_drho(rho, s)

  END FUNCTION dqx_drho




  FUNCTION h_function(y)

     IMPLICIT NONE
     REAL(DP)             :: y, y2, y4, h_function
     REAL(DP), PARAMETER  :: g1 = fpi/9.0D0                                     ! vdW-DF1/2
     REAL(DP), PARAMETER  :: a3 = 0.94950D0, g3 = 1.12D0, g32 = g3*g3           ! vdW-DF3-opt1
     REAL(DP), PARAMETER  :: a4 = 0.28248D0, g4 = 1.29D0, g42 = g4*g4           ! vdW-DF3-opt2
     REAL(DP), PARAMETER  :: a5 = 2.01059D0, b5 = 8.17471D0, g5 = 1.84981D0, &  ! vdW-DF-C6
                             AA = ( b5 + a5*(a5/2.0D0-g5) ) / ( 1.0D0+g5-a5 )   !


     y2 = y*y

     IF ( inlc == 1 .OR. inlc == 2 ) THEN

        h_function = 1.0D0 - EXP( -g1*y2 )
     
     ELSE IF ( inlc == 3 ) THEN

        y4 = y2*y2
        h_function = 1.0D0 - 1.0D0 / ( 1.0D0 + g3*y2 + g32*y4 + a3*y4*y4 )

     ELSE IF ( inlc == 4 ) THEN 

        y4 = y2*y2
        h_function = 1.0D0 - 1.0D0 / ( 1.0D0 + g4*y2 + g42*y4 + a4*y4*y4 )

     ELSE IF ( inlc == 5 ) THEN
     
        y4 = y2*y2
        h_function = 1.0D0 - ( 1.0D0 + ( (a5-g5)*y2 + AA*y4 ) / ( 1.0D0+AA*y2 ) ) * EXP( -a5*y2 ) 

     END IF

  END FUNCTION








  ! ####################################################################
  !                           |             |
  !                           |  XC_VDW_DF  |
  !                           |_____________|

  SUBROUTINE xc_vdW_DF (rho_valence, rho_core, etxc, vtxc, v)
  !
  !! Driver routine for vdW-DF calculations, called from Modules/funct.f90.
  !! The routine here sets up the parallel run (if any) and carry out the 
  !! calls necessary to calculate the non-local correlation contributions
  !! to the energy and potential.
  
  USE gvect,                 ONLY : ngm, g
  USE cell_base,             ONLY : omega, tpiba

  IMPLICIT NONE

  ! --------------------------------------------------------------------
  ! Local variables
  !                                               _
  REAL(DP), INTENT(IN)    :: rho_valence(:,:)    !
  REAL(DP), INTENT(IN)    :: rho_core(:)         !  PWSCF input variables
  REAL(DP), INTENT(INOUT) :: etxc, vtxc, v(:,:)  !_

  INTEGER :: i_grid, theta_i, i_proc             ! Indexing variables over grid points,
                                                 ! theta functions, and processors.

  REAL(DP) :: grid_cell_volume                   ! The volume of the unit cell per G-grid point.

  REAL(DP), ALLOCATABLE ::  q0(:)                ! The saturated value of q (equations 11 and 12
                                                 ! of DION). This saturation is that of
                                                 ! equation 5 in SOLER.

  REAL(DP), ALLOCATABLE :: grad_rho(:,:)         ! The gradient of the charge density. The
                                                 ! format is as follows:
                                                 ! grad_rho(cartesian_component,grid_point).

  REAL(DP), ALLOCATABLE :: potential(:)          ! The vdW contribution to the potential.

  REAL(DP), ALLOCATABLE :: dq0_drho(:)           ! The derivative of the saturated q0
                                                 ! (equation 5 of SOLER) with respect
                                                 ! to the charge density (see
                                                 ! get_q0_on_grid subroutine for details).

  REAL(DP), ALLOCATABLE :: dq0_dgradrho(:)       ! The derivative of the saturated q0
                                                 ! (equation 5 of SOLER) with respect
                                                 ! to the gradient of the charge density
                                                 ! (again, see get_q0_on_grid subroutine).

  COMPLEX(DP), ALLOCATABLE :: thetas(:,:)        ! These are the functions of equation 8 of
                                                 ! SOLER. They will be forward Fourier transformed
                                                 ! in place to get theta(k) and worked on in
                                                 ! place to get the u_alpha(r) of equation 11
                                                 ! in SOLER. They are formatted as follows:
                                                 ! thetas(grid_point, theta_i).

  REAL(DP) :: Ec_nl                              ! The non-local vdW contribution to the energy.

  REAL(DP), ALLOCATABLE :: total_rho(:)          ! This is the sum of the valence and core
                                                 ! charge. This just holds the piece assigned
                                                 ! to this processor.

  LOGICAL, SAVE :: first_iteration = .TRUE.      ! Whether this is the first time this
                                                 ! routine has been called.




  ! --------------------------------------------------------------------
  ! Write out the vdW-DF information and initialize the calculation.

  IF ( first_iteration ) THEN
     IF ( inlc > 5 ) CALL errore( 'xc_vdW_DF', 'inlc not implemented', 1 )
     CALL generate_kernel
     IF ( ionode ) CALL vdW_info
     first_iteration = .FALSE.
  END IF


  ! --------------------------------------------------------------------
  ! Allocate arrays. nnr is a PWSCF variable that holds the number of
  ! points assigned to a given processor.

  allocate( total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr),                &
            potential(dfftp%nnr), thetas(dfftp%nnr, Nqs),               &
            q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )


  ! --------------------------------------------------------------------
  ! Add together the valence and core charge densities to get the total
  ! charge density. Note that rho_core is not the true core density and
  ! it is only non-zero for pseudopotentials with non-local core
  ! corrections.

  total_rho = rho_valence(:,1) + rho_core(:)


  ! --------------------------------------------------------------------
  ! Here we calculate the gradient in reciprocal space using FFT.

  CALL fft_gradient_r2r (dfftp, total_rho, g, grad_rho)


  ! --------------------------------------------------------------------
  ! Find the value of q0 for all assigned grid points. q is defined in
  ! equations 11 and 12 of DION and q0 is the saturated version of q
  ! defined in equation 5 of SOLER. This routine also returns the
  ! derivatives of the q0s with respect to the charge-density and the
  ! gradient of the charge-density. These are needed for the potential
  ! calculated below. This routine also calculates the thetas.

  CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)


  ! --------------------------------------------------------------------
  ! Carry out the integration in equation 8 of SOLER. This also turns
  ! the theta arrays into the precursor to the u_i(k) array which is
  ! inverse fourier transformed to get the u_i(r) functions of SOLER
  ! equation 11. Add the energy we find to the output variable etxc.

  CALL vdW_DF_energy (thetas, Ec_nl)
  etxc = etxc + Ec_nl

  IF ( iverbosity > 0 ) THEN
     CALL mp_sum(Ec_nl, intra_bgrp_comm)
     IF ( ionode ) THEN
        WRITE(stdout,'(/ / A)')       "     -----------------------------------------------"
        WRITE(stdout,'(A, F15.8, A)') "     Non-local corr. energy    =  ", Ec_nl, " Ry"
        WRITE(stdout,'(A /)')         "     -----------------------------------------------"
     END IF
  END IF


  ! --------------------------------------------------------------------
  ! Here we calculate the potential. This is calculated via equation 10
  ! of SOLER, using the u_i(r) calculated from quations 11 and 12 of
  ! SOLER. Each processor allocates the array to be the size of the full
  ! grid because, as can be seen in SOLER equation 10, processors need
  ! to access grid points outside their allocated regions. Begin by
  ! FFTing the u_i(k) to get the u_i(r) of SOLER equation 11.

  DO theta_i = 1, Nqs
     CALL invfft('Rho', thetas(:,theta_i), dfftp)
  END DO

  CALL vdW_DF_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, thetas, potential)
  v(:,1) = v(:,1) + e2 * potential(:)


  ! --------------------------------------------------------------------
  ! The integral of rho(r)*potential(r) for the vtxc output variable.

  grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)

  DO i_grid = 1, dfftp%nnr
     vtxc = vtxc + e2 * grid_cell_volume * rho_valence(i_grid,1) * potential(i_grid)
  END DO

  DEALLOCATE ( total_rho, grad_rho, potential, thetas, q0, dq0_drho, dq0_dgradrho )

  END SUBROUTINE xc_vdW_DF








  ! ####################################################################
  !                          |                  |
  !                          |  XC_VDW_DF_spin  |
  !                          |__________________|
  !
  !
  SUBROUTINE xc_vdW_DF_spin (rho_valence, rho_core, etxc, vtxc, v)
  
  !! This subroutine is as similar to \(\texttt{xc_vdW_DF}\) as possible,
  !! but handles the collinear \(\text{nspin}=2\) case.
  
  USE gvect,                 ONLY : ngm, g
  USE cell_base,             ONLY : omega, tpiba

  IMPLICIT NONE

  ! --------------------------------------------------------------------
  ! Local variables
  !                                              _
  REAL(DP), INTENT(IN) :: rho_valence(:,:)      !
  REAL(DP), INTENT(IN) :: rho_core(:)           ! PWSCF input variables.
  REAL(DP), INTENT(INOUT) :: etxc, vtxc, v(:,:) !_


  INTEGER :: i_grid, theta_i, i_proc            ! Indexing variables over grid points,
                                                ! theta functions, and processors, and a
                                                ! generic index.

  REAL(DP) :: grid_cell_volume                  ! The volume of the unit cell per G-grid point.

  REAL(DP), ALLOCATABLE :: q0(:)                ! The saturated value of q (equations 11 and 12
                                                ! of DION). This saturation is that of
                                                ! equation 5 in SOLER.

  REAL(DP), ALLOCATABLE :: grad_rho(:,:)        ! The gradient of the charge density. The
                                                ! format is as follows:
                                                ! grad_rho(cartesian_component,grid_point).
  REAL(DP), ALLOCATABLE :: grad_rho_up(:,:)     ! The gradient of the up charge density.
                                                ! Same format as grad_rho.
  REAL(DP), ALLOCATABLE :: grad_rho_down(:,:)   ! The gradient of the down charge density.
                                                ! Same format as grad_rho.

  REAL(DP), ALLOCATABLE :: potential_up(:)      ! The vdW contribution to the potential.
  REAL(DP), ALLOCATABLE :: potential_down(:)    ! The vdW contribution to the potential.

  REAL(DP), ALLOCATABLE :: dq0_drho_up(:)       ! The derivative of the saturated q0
  REAL(DP), ALLOCATABLE :: dq0_drho_down(:)     ! (equation 5 of SOLER) with respect
                                                ! to the charge density (see
                                                ! get_q0_on_grid subroutine for details).

  REAL(DP), ALLOCATABLE :: dq0_dgradrho_up(:)   ! The derivative of the saturated q0
  REAL(DP), ALLOCATABLE :: dq0_dgradrho_down(:) ! (equation 5 of SOLER) with respect
                                                ! to the gradient of the charge density
                                                ! (again, see get_q0_on_grid subroutine).

  COMPLEX(DP), ALLOCATABLE :: thetas(:,:)       ! These are the functions of equation 8 of
                                                ! SOLER. They will be forward Fourier transformed
                                                ! in place to get theta(k) and worked on in
                                                ! place to get the u_alpha(r) of equation 11
                                                ! in SOLER. They are formatted as follows:
                                                ! thetas(grid_point, theta_i).

  REAL(DP) :: Ec_nl                             ! The non-local vdW contribution to the energy.

  REAL(DP), ALLOCATABLE :: total_rho(:)         ! This is the sum of the valence (up and down)
                                                ! and core charge. This just holds the piece
                                                ! assigned to this processor.
  REAL(DP), ALLOCATABLE :: rho_up(:)            ! This is the just the up valence charge.
                                                ! This just holds the piece assigned
                                                ! to this processor.
  REAL(DP), ALLOCATABLE :: rho_down(:)          ! This is the just the down valence charge.
                                                ! This just holds the piece assigned
                                                ! to this processor.

  LOGICAL, SAVE :: first_iteration = .TRUE.     ! Whether this is the first time this
                                                ! routine has been called.




  ! --------------------------------------------------------------------
  ! Write out the vdW-DF information and initialize the calculation.

  IF ( first_iteration ) THEN
     IF ( inlc > 5 ) CALL errore( 'xc_vdW_DF_spin', 'inlc not implemented', 1 )
     CALL generate_kernel
     IF ( ionode ) CALL vdW_info
     first_iteration = .FALSE.
  END IF


  ! --------------------------------------------------------------------
  ! Allocate arrays. nnr is a PWSCF variable that holds the number of
  ! points assigned to a given processor.

  ALLOCATE( total_rho(dfftp%nnr), rho_up(dfftp%nnr), rho_down(dfftp%nnr),         &
     grad_rho(3,dfftp%nnr), grad_rho_up(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr), &
     potential_up(dfftp%nnr), potential_down(dfftp%nnr), thetas(dfftp%nnr, Nqs),  &
     q0(dfftp%nnr), dq0_drho_up(dfftp%nnr), dq0_dgradrho_up(dfftp%nnr),           &
     dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )


  ! --------------------------------------------------------------------
  ! Add together the valence and core charge densities to get the total
  ! charge density. Note that rho_core is not the true core density and
  ! it is only non-zero for pseudopotentials with non-local core
  ! corrections.

  rho_up    = ( rho_valence(:,1) + rho_valence(:,2) + rho_core(:) )*0.5D0
  rho_down  = ( rho_valence(:,1) - rho_valence(:,2) + rho_core(:) )*0.5D0
  total_rho = rho_up + rho_down

#if defined (__SPIN_BALANCED)
  rho_up   = total_rho*0.5D0
  rho_down = rho_up
  WRITE(stdout,'(/,/,"     Performing spin-balanced Ecnl calculation!")')
#endif


  ! --------------------------------------------------------------------
  ! Here we calculate the gradient in reciprocal space using FFT.

  CALL fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
  CALL fft_gradient_r2r (dfftp, rho_up,    g, grad_rho_up)
  CALL fft_gradient_r2r (dfftp, rho_down,  g, grad_rho_down)


  ! --------------------------------------------------------------------
  ! Find the value of q0 for all assigned grid points. q is defined in
  ! equations 11 and 12 of DION and q0 is the saturated version of q
  ! defined in equation 5 of SOLER. In the spin case, q0 is defined by
  ! equation 8 (and text above that equation) of THONHAUSER. This
  ! routine also returns the derivatives of the q0s with respect to the
  ! charge-density and the gradient of the charge-density. These are
  ! needed for the potential calculated below.

  CALL get_q0_on_grid_spin (total_rho, rho_up, rho_down, grad_rho, &
       grad_rho_up, grad_rho_down, q0, dq0_drho_up, dq0_drho_down, &
       dq0_dgradrho_up, dq0_dgradrho_down, thetas)


  ! --------------------------------------------------------------------
  ! Carry out the integration in equation 8 of SOLER. This also turns
  ! the thetas array into the precursor to the u_i(k) array which is
  ! inverse fourier transformed to get the u_i(r) functions of SOLER
  ! equation 11. Add the energy we find to the output variable etxc.

  CALL vdW_DF_energy(thetas, Ec_nl)
  etxc = etxc + Ec_nl

  IF ( iverbosity > 0 ) THEN
     CALL mp_sum(Ec_nl, intra_bgrp_comm)
     IF (ionode) THEN
        WRITE(stdout,'(/ / A)')       "     -----------------------------------------------"
        WRITE(stdout,'(A, F15.8, A)') "     Non-local corr. energy    =  ", Ec_nl, " Ry"
        WRITE(stdout,'(A /)')         "     -----------------------------------------------"
     END IF
  END IF


  ! --------------------------------------------------------------------
  ! Here we calculate the potential. This is calculated via equation 10
  ! of SOLER, using the u_i(r) calculated from quations 11 and 12 of
  ! SOLER. Each processor allocates the array to be the size of the full
  ! grid because, as can be seen in SOLER equation 10, processors need
  ! to access grid points outside their allocated regions. Begin by
  ! FFTing the u_i(k) to get the u_i(r) of SOLER equation 11.

  DO theta_i = 1, Nqs
     CALL invfft('Rho', thetas(:,theta_i), dfftp)
  END DO

  CALL vdW_DF_potential (q0, dq0_drho_up  , dq0_dgradrho_up  , grad_rho_up  , thetas, potential_up  )
  CALL vdW_DF_potential (q0, dq0_drho_down, dq0_dgradrho_down, grad_rho_down, thetas, potential_down)

  v(:,1) = v(:,1) + e2 * potential_up  (:)
  v(:,2) = v(:,2) + e2 * potential_down(:)


  ! --------------------------------------------------------------------
  ! The integral of rho(r)*potential(r) for the vtxc output variable

  grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)

  DO i_grid = 1, dfftp%nnr
     vtxc = vtxc + e2 * grid_cell_volume * (rho_valence(i_grid,1) +   &
            rho_valence(i_grid,2)) * 0.5_dp * potential_up  (i_grid)  &
                 + e2 * grid_cell_volume * (rho_valence(i_grid,1) -   &
            rho_valence(i_grid,2)) * 0.5_dp * potential_down(i_grid)
  END DO

  DEALLOCATE( total_rho, rho_up, rho_down, grad_rho, grad_rho_up, grad_rho_down, &
              potential_up, potential_down, thetas,                              &
              q0, dq0_drho_up, dq0_dgradrho_up, dq0_drho_down, dq0_dgradrho_down )

  END SUBROUTINE xc_vdW_DF_spin








  ! ####################################################################
  !                       |                  |
  !                       |  GET_Q0_ON_GRID  |
  !                       |__________________|
  !
  !
  SUBROUTINE get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
  !! This routine first calculates the q value defined in (DION equations
  !! 11 and 12), then saturates it according to (SOLER equation 5). More
  !! specifically it calculates the following:
  !
  !! * \(\text{q0(ir)}\) = saturated value of q;
  !! * \(\text{dq0_drho(ir)} = \text{total_rho} * d\text{q0} /d\text{rho}\)
  !! * \(\text{dq0_dgradrho} = \text{total_rho} / |\text{grad_rho}| * d\text{q0}
  !!   / d |\text{grad_rho}|\)
  
  
  IMPLICIT NONE

  REAL(DP), INTENT(IN)      :: total_rho(:), grad_rho(:,:)         ! Input variables needed.

  REAL(DP), INTENT(OUT)     :: q0(:), dq0_drho(:), dq0_dgradrho(:) ! Output variables that have been allocated
                                                                   ! outside this routine but will be set here.
  COMPLEX(DP), INTENT(INOUT):: thetas(:,:)                         ! The thetas from SOLER.

  INTEGER, PARAMETER        :: m_cut = 12                          ! How many terms to include in the sum
                                                                   ! of SOLER equation 5.

  REAL(DP)                  :: rho                                 ! Local variable for the density.
  REAL(DP)                  :: r_s                                 ! Wigner-Seitz radius.
  REAL(DP)                  :: s                                   ! Reduced gradient.
  REAL(DP)                  :: q
  REAL(DP)                  :: ec
  REAL(DP)                  :: dq0_dq                              ! The derivative of the saturated
                                                                   ! q0 with respect to q.

  INTEGER                   :: i_grid                              ! Indexing variable.




  ! --------------------------------------------------------------------
  ! Initialize q0-related arrays.

  q0(:)           = q_cut
  dq0_drho(:)     = 0.0D0
  dq0_dgradrho(:) = 0.0D0


  DO i_grid = 1, dfftp%nnr

     rho = total_rho(i_grid)


     ! -----------------------------------------------------------------
     ! This prevents numerical problems. If the charge density is
     ! negative (an unphysical situation), we simply treat it as very
     ! small. In that case, q0 will be very large and will be saturated.
     ! For a saturated q0 the derivative dq0_dq will be 0 so we set q0 =
     ! q_cut and dq0_drho = dq0_dgradrho = 0 and go on to the next
     ! point.

     IF ( rho < epsr ) CYCLE


     ! -----------------------------------------------------------------
     ! Calculate some intermediate values needed to find q.

     r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)

     s   = SQRT( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 ) / &
           (2.0D0 * kF(rho) * rho )


     ! -----------------------------------------------------------------
     ! This is the q value defined in equations 11 and 12 of DION.
     ! Use pw() from flib/functionals.f90 to get qc = kf/eps_x * eps_c.
     !
     CALL pw(r_s, 1, ec, dq0_drho(i_grid))
     q = -4.0D0*pi/3.0D0 * ec + kF(rho) * Fs(s)


     ! -----------------------------------------------------------------
     ! Bring q into its proper bounds.

     CALL saturate_q ( q, q_cut, q0(i_grid), dq0_dq )
     IF (q0(i_grid) < q_min) q0(i_grid) = q_min


     ! -----------------------------------------------------------------
     ! Here we find derivatives. These are actually the density times
     ! the derivative of q0 with respect to rho and grad_rho. The
     ! density factor comes in since we are really differentiating
     ! theta = (rho)*P(q0) with respect to density (or its gradient)
     ! which will be
     !
     !    dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho]
     !
     ! and
     !
     !    dtheta_dgrad_rho = dP_dq0 * [rho * dq0_dq * dq_dgrad_rho]
     !
     ! The parts in square brackets are what is calculated here. The
     ! dP_dq0 term will be interpolated later.

     dq0_drho(i_grid)     = dq0_dq * rho * ( -4.0D0*pi/3.0D0 * &
                            (dq0_drho(i_grid) - ec)/rho + dqx_drho(rho, s) )
     dq0_dgradrho(i_grid) = dq0_dq * rho * kF(rho) * dFs_ds(s) * ds_dgradrho(rho)

  END DO


  ! --------------------------------------------------------------------
  ! Here we calculate the theta functions of SOLER equation 8. These are
  ! defined as
  !
  !    rho * P_i(q0(rho, grad_rho))
  !
  ! where P_i is a polynomial that interpolates a Kroneker delta
  ! function at the point q_i (taken from the q_mesh) and q0 is the
  ! saturated version of q. q is defined in equations 11 and 12 of DION
  ! and the saturation proceedure is defined in equation 5 of SOLER.
  ! This is the biggest memory consumer in the method since the thetas
  ! array is (total # of FFT points)*Nqs complex numbers. In a parallel
  ! run, each processor will hold the values of all the theta functions
  ! on just the points assigned to it. thetas are stored in reciprocal
  ! space as theta_i(k) because this is the way they are used later for
  ! the convolution (equation 8 of SOLER). Start by interpolating the
  ! P_i polynomials defined in equation 3 in SOLER for the particular q0
  ! values we have.

  CALL spline_interpolation (q_mesh, q0, thetas)

  DO i_grid = 1, dfftp%nnr
     thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
  END DO

  DO idx = 1, Nqs
     CALL fwfft ('Rho', thetas(:,idx), dfftp)
  END DO

  END SUBROUTINE get_q0_on_grid


  !-----------------------------------------------------------------------
  SUBROUTINE pw( rs, iflag, ec, vc )
    !-----------------------------------------------------------------------
    ! --A provisional copy of the pw routine in XC lib to avoid external calls--
    !! * iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
    !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994)
    !
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: rs
    INTEGER, INTENT(IN)  :: iflag
    REAL(DP), INTENT(OUT) :: ec, vc
    !
    REAL(DP), PARAMETER :: a=0.031091d0, b1=7.5957d0, b2=3.5876d0, c0=a, &
                           c1=0.046644d0, c2=0.00664d0, c3=0.01043d0, d0=0.4335d0, &
                           d1=1.4408d0
    REAL(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog
    REAL(DP) :: a1(2), b3(2), b4(2)
    DATA a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, &
         b4 / 0.49294d0, 0.13354d0 /
    ! high- and low-density formulae implemented but not used in PW case
    ! (reason: inconsistencies in PBE/PW91 functionals).
    IF ( rs < 1d0 .AND. iflag == 2 ) THEN
       ! high density formula
       lnrs = LOG(rs)
       ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs
       vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * &
                 lnrs - (2.d0 * c3 + c2) / 3.d0 * rs
    ELSEIF ( rs > 100.d0 .AND. iflag == 2 ) THEN
       ! low density formula
       ec = - d0 / rs + d1 / rs**1.5d0
       vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0
    ELSE
       ! interpolation formula
       rs12 = SQRT(rs)
       rs32 = rs * rs12
       rs2  = rs**2
       om   = 2.d0*a*( b1*rs12 + b2*rs + b3(iflag) * rs32 + b4(iflag)*rs2 )
       dom  = 2.d0*a*( 0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3(iflag) * &
              rs32 + 2.d0 * b4(iflag) * rs2 )
       olog = LOG( 1.d0 + 1.0d0 / om )
       !
       ec = - 2.d0 * a * (1.d0 + a1(iflag) * rs) * olog
       vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1(iflag) * rs) &
                * olog - 2.d0 / 3.d0 * a * (1.d0 + a1(iflag) * rs) * dom / &
                (om * (om + 1.d0) )
    ENDIF
    !
    RETURN
    !
  END SUBROUTINE pw
  !---------------------

  ! ####################################################################
  !                       |                       |
  !                       |  GET_Q0_ON_GRID_spin  |
  !                       |_______________________|

  SUBROUTINE get_q0_on_grid_spin (total_rho, rho_up, rho_down, grad_rho, &
             grad_rho_up, grad_rho_down, q0, dq0_drho_up, dq0_drho_down, &
             dq0_dgradrho_up, dq0_dgradrho_down, thetas)
  
  !! Find the value of \(\text{q0}\) for all assigned grid points. q is defined in
  !! equations 11 and 12 of DION and q0 is the saturated version of q
  !! defined in equation 5 of SOLER. In the spin case, q0 is defined by
  !! equation 8 (and text above that equation) of THONHAUSER. This
  !! routine also returns the derivatives of the \(\text{q0}\)s with respect to the
  !! charge-density and the gradient of the charge-density. These are
  !! needed for the potential.
  
  IMPLICIT NONE

  REAL(DP),  INTENT(IN)      :: total_rho(:), grad_rho(:,:)              ! Input variables.
  REAL(DP),  INTENT(IN)      :: rho_up(:), grad_rho_up(:,:)              ! Input variables.
  REAL(DP),  INTENT(IN)      :: rho_down(:), grad_rho_down(:,:)          ! Input variables.

  REAL(DP),  INTENT(OUT)     :: q0(:), dq0_drho_up(:), dq0_drho_down(:)  ! Output variables.
  REAL(DP),  INTENT(OUT)     :: dq0_dgradrho_up(:), dq0_dgradrho_down(:) ! Output variables.
  COMPLEX(DP), INTENT(INOUT) :: thetas(:,:)                              ! The thetas from SOLER.

  REAL(DP)                   :: rho, up, down                            ! Local copy of densities.
  REAL(DP)                   :: zeta                                     ! Spin polarization.
  REAL(DP)                   :: r_s                                      ! Wigner-Seitz radius.
  REAL(DP)                   :: q, qc, qx, qx_up, qx_down                ! q for exchange and correlation.
  REAL(DP)                   :: q0x_up, q0x_down                         ! Saturated q values.
  REAL(DP)                   :: fac
  REAL(DP)                   :: ec, vc(2)
  REAL(DP)                   :: dq0_dq, dq0x_up_dq, dq0x_down_dq         ! Derivative of q0 w.r.t q.
  REAL(DP)                   :: dqc_drho_up, dqc_drho_down               ! Intermediate values.
  REAL(DP)                   :: dqx_drho_up, dqx_drho_down               ! Intermediate values.
  REAL(DP)                   :: s_up, s_down                             ! Reduced gradients.
  INTEGER                    :: i_grid                                   ! Indexing variable.
  LOGICAL                    :: calc_qx_up, calc_qx_down


  fac = 2.0D0**(-1.0D0/3.0D0)


  ! --------------------------------------------------------------------
  ! Initialize q0-related arrays.

  q0(:)                = q_cut
  dq0_drho_up(:)       = 0.0D0
  dq0_drho_down(:)     = 0.0D0
  dq0_dgradrho_up(:)   = 0.0D0
  dq0_dgradrho_down(:) = 0.0D0


  DO i_grid = 1, dfftp%nnr

     rho  = total_rho(i_grid)
     up   = rho_up(i_grid)
     down = rho_down(i_grid)


     ! -----------------------------------------------------------------
     ! This prevents numerical problems. If the charge density is
     ! negative (an unphysical situation), we simply treat it as very
     ! small. In that case, q0 will be very large and will be saturated.
     ! For a saturated q0 the derivative dq0_dq will be 0 so we set q0 =
     ! q_cut and dq0_drho = dq0_dgradrho = 0 and go on to the next
     ! point.

     IF ( rho < epsr ) CYCLE

     calc_qx_up   = .TRUE.
     calc_qx_down = .TRUE.

     IF ( up   < epsr/2.0D0 ) calc_qx_up   = .FALSE.
     IF ( down < epsr/2.0D0 ) calc_qx_down = .FALSE.


     ! -----------------------------------------------------------------
     ! The spin case is numerically even more tricky and we have to
     ! saturate each spin channel separately. Note that we are
     ! saturating at a higher value here, so that very large q values
     ! get saturated to exactly q_cut in the second, overall saturation.

     q0x_up        = 0.0D0
     q0x_down      = 0.0D0
     dqx_drho_up   = 0.0D0
     dqx_drho_down = 0.0D0


     IF ( calc_qx_up ) THEN
        s_up    = SQRT( grad_rho_up(1,i_grid)**2 + grad_rho_up(2,i_grid)**2 + &
                  grad_rho_up(3,i_grid)**2 ) / (2.0D0 * kF(up) * up)
        qx_up   = kF(2.0D0*up) * Fs(fac*s_up)
        CALL saturate_q (qx_up, 4.0D0*q_cut, q0x_up, dq0x_up_dq)
     END IF

     IF ( calc_qx_down ) THEN
        s_down  = SQRT( grad_rho_down(1,i_grid)**2 + grad_rho_down(2,i_grid)**2 + &
                  grad_rho_down(3,i_grid)**2) / (2.0D0 * kF(down) * down)
        qx_down = kF(2.0D0*down) * Fs(fac*s_down)
        CALL saturate_q (qx_down, 4.0D0*q_cut, q0x_down, dq0x_down_dq)
     END IF


     ! -----------------------------------------------------------------
     ! This is the q value defined in equations 11 and 12 of DION and
     ! equation 8 of THONHAUSER (also see text above that equation).

     r_s  = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
     zeta = (up - down) / rho
     IF ( ABS(zeta) > 1.0D0 ) zeta = SIGN(1.0D0, zeta)
     call pw_spin( r_s, zeta, ec, vc(1), vc(2) )
     dqc_drho_up   = vc(1)
     dqc_drho_down = vc(2)

     qx = ( up * q0x_up + down * q0x_down ) / rho
     qc = -4.0D0*pi/3.0D0 * ec
     q  = qx + qc


     ! -----------------------------------------------------------------
     ! Bring q into its proper bounds.

     CALL saturate_q (q, q_cut, q0(i_grid), dq0_dq)
     IF (q0(i_grid) < q_min) q0(i_grid) = q_min


     ! -----------------------------------------------------------------
     ! Here we find derivatives. These are actually the density times
     ! the derivative of q0 with respect to rho and grad_rho. The
     ! density factor comes in since we are really differentiating
     ! theta = (rho)*P(q0) with respect to density (or its gradient)
     ! which will be
     !
     !    dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho]
     !
     ! and
     !
     !    dtheta_dgrad_rho = dP_dq0 * [rho * dq0_dq * dq_dgrad_rho]
     !
     ! The parts in square brackets are what is calculated here. The
     ! dP_dq0 term will be interpolated later.

     IF ( calc_qx_up ) THEN
        dqx_drho_up   = 2.0D0*dq0x_up_dq*up*dqx_drho(2.0D0*up, fac*s_up) + q0x_up*down/rho
        dq0_dgradrho_up (i_grid) = 2.0D0 * dq0_dq * dq0x_up_dq * up * kF(2.0D0*up) * &
                        dFs_ds(fac*s_up) * ds_dgradrho(2.0D0*up)
     END IF

     IF ( calc_qx_down ) THEN
        dqx_drho_down = 2.0D0*dq0x_down_dq*down*dqx_drho(2.0D0*down, fac*s_down) + q0x_down*up/rho
        dq0_dgradrho_down(i_grid) = 2.0D0 * dq0_dq * dq0x_down_dq * down * kF(2.0D0*down) * &
                        dFs_ds(fac*s_down) * ds_dgradrho(2.0D0*down)
     END IF

     IF ( calc_qx_down ) dqx_drho_up   = dqx_drho_up   - q0x_down*down/rho
     IF ( calc_qx_up )   dqx_drho_down = dqx_drho_down - q0x_up  *up  /rho

     dqc_drho_up   = -4.0D0*pi/3.0D0 * (dqc_drho_up   - ec)
     dqc_drho_down = -4.0D0*pi/3.0D0 * (dqc_drho_down - ec)

     dq0_drho_up  (i_grid) = dq0_dq * (dqc_drho_up   + dqx_drho_up  )
     dq0_drho_down(i_grid) = dq0_dq * (dqc_drho_down + dqx_drho_down)

  END DO


  ! --------------------------------------------------------------------
  ! Here we calculate the theta functions of SOLER equation 8. These are
  ! defined as
  !
  !    rho * P_i(q0(rho, grad_rho))
  !
  ! where P_i is a polynomial that interpolates a Kroneker delta
  ! function at the point q_i (taken from the q_mesh) and q0 is the
  ! saturated version of q. q is defined in equations 11 and 12 of DION
  ! and the saturation proceedure is defined in equation 5 of SOLER.
  ! This is the biggest memory consumer in the method since the thetas
  ! array is (total # of FFT points)*Nqs complex numbers. In a parallel
  ! run, each processor will hold the values of all the theta functions
  ! on just the points assigned to it. thetas are stored in reciprocal
  ! space as theta_i(k) because this is the way they are used later for
  ! the convolution (equation 8 of SOLER). Start by interpolating the
  ! P_i polynomials defined in equation 3 in SOLER for the particular q0
  ! values we have.

  CALL spline_interpolation (q_mesh, q0, thetas)

  DO i_grid = 1, dfftp%nnr
     thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
  END DO

  DO idx = 1, Nqs
     CALL fwfft ('Rho', thetas(:,idx), dfftp)
  END Do

  END SUBROUTINE get_q0_on_grid_spin

  
  !-----------------------------------------------------------------------
  SUBROUTINE pw_spin( rs, zeta, ec, vc_up, vc_dw )
    !-----------------------------------------------------------------------
    !--A provisional copy of the pw routine in XC lib to avoid external calls--
    !! J.P. Perdew and Y. Wang, PRB 45, 13244 (1992).
    IMPLICIT NONE
    REAL(DP), INTENT(IN) :: rs
    !! Wigner-Seitz radius
    REAL(DP), INTENT(IN) :: zeta
    !! zeta = (rho_up - rho_dw)/rho_tot
    REAL(DP), INTENT(OUT) :: ec, vc_up, vc_dw
    !
    REAL(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz
    REAL(DP) :: om, dom, olog, epwc, vpwc
    REAL(DP) :: omp, domp, ologp, epwcp, vpwcp
    REAL(DP) :: oma, doma, ologa, alpha, vpwca
    !
    ! xc parameters, unpolarised
    REAL(DP), PARAMETER :: a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = &
             3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, &
             c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0
    ! xc parameters, polarised
    REAL(DP), PARAMETER :: ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p &
                 = 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p =     &
                0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p &
                = 1.7697d0
    ! xc PARAMETERs, antiferro
    REAL(DP), PARAMETER :: aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = &
                 3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a =       &
                 0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a &
                 = 0.3969d0
    REAL(DP), PARAMETER :: fz0 = 1.709921d0
    !
    !     if (rs < 0.5d0) then
    ! high density formula (not implemented)
    !
    !     elseif (rs > 100.d0) then
    ! low density formula (not implemented)
    !
    !     else
    ! interpolation formula
    !
    zeta2 = zeta * zeta
    zeta3 = zeta2 * zeta
    zeta4 = zeta3 * zeta
    rs12 = SQRT(rs)
    rs32 = rs * rs12
    rs2 = rs**2
    ! unpolarised
    om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2)
    dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 * rs32 &
         + 2.d0 * b4 * rs2)
    olog = LOG(1.d0 + 1.0d0 / om)
    epwc = - 2.d0 * a * (1.d0 + a1 * rs) * olog
    vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs) * olog - 2.d0 / &
           3.d0 * a * (1.d0 + a1 * rs) * dom / (om * (om + 1.d0) )
    ! polarized
    omp  = 2.d0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2)
    domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs + 1.5d0 * b3p * &
           rs32 + 2.d0 * b4p * rs2)
    ologp = LOG(1.d0 + 1.0d0 / omp)
    epwcp = - 2.d0 * ap * (1.d0 + a1p * rs) * ologp
    vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs) * ologp - &
            2.d0 / 3.d0 * ap * (1.d0 + a1p * rs) * domp / (omp * (omp + 1.d0))
    ! antiferro
    oma = 2.d0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2)
    doma = 2.d0 * aa * ( 0.5d0 * b1a * rs12 + b2a * rs + 1.5d0 * b3a * &
           rs32 + 2.d0 * b4a * rs2 )
    ologa = LOG( 1.d0 + 1.0d0/oma )
    alpha = 2.d0 * aa * (1.d0 + a1a*rs) * ologa
    vpwca = + 2.d0 * aa * (1.d0 + 2.d0/3.d0 * a1a * rs) * ologa + &
            2.d0 / 3.d0 * aa * (1.d0 + a1a*rs) * doma / (oma * (oma + 1.d0))
    !
    fz = ( (1.d0 + zeta)**(4.d0 / 3.d0) + (1.d0 - zeta)**(4.d0 / &
            3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0)
    dfz = ( (1.d0 + zeta)**(1.d0 / 3.d0) - (1.d0 - zeta)**(1.d0 / &
            3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) )
    !
    ec = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) &
                * fz * zeta4
    vc_up = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
                   * fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
                   zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
                   * (1.d0 - zeta)
    vc_dw = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
                   * fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
                   zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
                   * (1.d0 + zeta)
    RETURN
    !
  END SUBROUTINE pw_spin

  

  ! ####################################################################
  !                            |              |
  !                            |  saturate_q  |
  !                            |______________|

  SUBROUTINE saturate_q (q, q_cutoff, q0, dq0_dq)
  
  !! Here, we calculate \(\text{q0}\) by saturating \(\text{q}\) according to
  !! equation 5 of SOLER. Also, we find the derivative \(\text{dq0/dq}\) needed
  !! for the derivatives \(\text{dq0/drho}\) and \(\text{dq0/dgradrh0}\).
  
  IMPLICIT NONE

  REAL(DP),  INTENT(IN)      :: q             ! Input q.
  REAL(DP),  INTENT(IN)      :: q_cutoff      ! Cutoff q.
  REAL(DP),  INTENT(OUT)     :: q0            ! Output saturated q.
  REAL(DP),  INTENT(OUT)     :: dq0_dq        ! Derivative of dq0/dq.

  REAL(DP)                   :: e_exp         ! Exponent.
  INTEGER,   PARAMETER       :: m_cut = 12    ! How many terms to include in
                                              ! the sum of SOLER equation 5.
  ! --------------------------------------------------------------------

  e_exp  = 0.0D0
  dq0_dq = 0.0D0

  DO idx = 1, m_cut
     e_exp  = e_exp + (q/q_cutoff)**idx/idx
     dq0_dq = dq0_dq + (q/q_cutoff)**(idx-1)
  END Do

  q0     = q_cutoff*(1.0D0 - EXP(-e_exp))
  dq0_dq = dq0_dq * EXP(-e_exp)

  END SUBROUTINE saturate_q








  ! ####################################################################
  !                          |               |
  !                          | vdW_DF_energy |
  !                          |_______________|
  !
  !
  SUBROUTINE vdW_DF_energy (thetas, vdW_xc_energy)

  !! This routine carries out the integration of equation 8 of SOLER.  It
  !! returns the non-local exchange-correlation energy and the \(\text{u_alpha(k)}\)
  !! arrays used to find the \(\text{u_alpha(r)}\) arrays via equations 11 and 12 in
  !! SOLER.
  
  USE gvect,           ONLY : gg, ngm, igtongl, gl, ngl, gstart
  USE cell_base,       ONLY : tpiba, omega

  IMPLICIT NONE

  COMPLEX(DP), INTENT(INOUT) :: thetas(:,:)            ! On input this variable holds the theta
                                                       ! functions (equation 8, SOLER) in the format
                                                       ! thetas(grid_point, theta_i). On output
                                                       ! this array holds u_alpha(k) =
                                                       ! Sum_j[theta_beta(k)phi_alpha_beta(k)].

  REAL(DP), INTENT(OUT) :: vdW_xc_energy               ! The non-local correlation energy.

  REAL(DP), ALLOCATABLE :: kernel_of_k(:,:)            ! This array will hold the interpolated kernel
                                                       ! values for each pair of q values in the q_mesh.

  REAL(DP)    :: g                                     ! The magnitude of the current g vector.
  INTEGER     :: last_g                                ! The shell number of the last g vector.


  INTEGER     :: g_i, q1_i, q2_i, i_grid               ! Index variables.

  COMPLEX(DP) :: theta(Nqs), thetam(Nqs), theta_g(Nqs) ! Temporary storage arrays used since we
                                                       ! are overwriting the thetas array here.
  REAL(DP)    :: G0_term, G_multiplier

  COMPLEX(DP), ALLOCATABLE :: u_vdw(:,:)               ! Temporary array holding u_alpha(k).




  ALLOCATE ( u_vdW(dfftp%nnr,Nqs), kernel_of_k(Nqs, Nqs) )
  vdW_xc_energy = 0.0D0
  u_vdW(:,:)    = CMPLX(0.0_DP, 0.0_DP, kind=dp)


  ! --------------------------------------------------------------------
  ! Loop over PWSCF's array of magnitude-sorted g-vector shells. For
  ! each shell, interpolate the kernel at this magnitude of g, then find
  ! all points on the shell and carry out the integration over those
  ! points. The PWSCF variables used here are ngm = number of g-vectors
  ! on this processor, nl = an array that gives the indices into the FFT
  ! grid for a particular g vector, igtongl = an array that gives the
  ! index of which shell a particular g vector is in, gl = an array that
  ! gives the magnitude of the g vectors for each shell. In essence, we
  ! are forming the reciprocal-space u(k) functions of SOLER equation
  ! 11. These are kept in thetas array. Here we should use gstart,ngm
  ! but all the cases are handled by conditionals inside the loop

  G_multiplier = 1.0D0
  IF ( gamma_only ) G_multiplier = 2.0D0

  last_g = -1

  DO g_i = 1, ngm

     IF ( igtongl(g_i) .NE. last_g) THEN
        g = SQRT(gl(igtongl(g_i))) * tpiba
        CALL interpolate_kernel(g, kernel_of_k)
        last_g = igtongl(g_i)
     END IF

     theta = thetas(dfftp%nl(g_i),:)

     DO q2_i = 1, Nqs
        DO q1_i = 1, Nqs
           u_vdW(dfftp%nl(g_i),q2_i)  = u_vdW(dfftp%nl(g_i),q2_i) + kernel_of_k(q2_i,q1_i)*theta(q1_i)
        END DO
        vdW_xc_energy = vdW_xc_energy + G_multiplier * (u_vdW(dfftp%nl(g_i),q2_i)*CONJG(theta(q2_i)))
     END DO

     IF (g_i < gstart ) vdW_xc_energy = vdW_xc_energy / G_multiplier

  END DO

  IF ( gamma_only ) u_vdW(dfftp%nlm(:),:) = CONJG( u_vdW(dfftp%nl(:),:) )


  ! --------------------------------------------------------------------
  ! Apply scaling factors. The e2 comes from PWSCF's choice of units.
  ! This should be 0.5 * e2 * vdW_xc_energy * (2pi)^3/omega * (omega)^2,
  ! with the (2pi)^3/omega being the volume element for the integral
  ! (the volume of the reciprocal unit cell) and the 2 factors of omega
  ! being used to cancel the factor of 1/omega PWSCF puts on forward
  ! FFTs of the 2 theta factors. 1 omega cancels and the (2pi)^3
  ! cancels because there should be a factor of 1/(2pi)^3 on the radial
  ! Fourier transform of phi that was left out to cancel with this
  ! factor.

  vdW_xc_energy = 0.5D0 * e2 * omega * vdW_xc_energy

  thetas(:,:) = u_vdW(:,:)
  DEALLOCATE ( u_vdW, kernel_of_k )

  END SUBROUTINE vdW_DF_energy








  ! ####################################################################
  !                        |                   |
  !                        |  vdW_DF_potential |
  !                        |___________________|
  !
  !
  SUBROUTINE vdW_DF_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, u_vdW, potential)

  !! This routine finds the non-local correlation contribution to the potential (i.e.
  !! the derivative of the non-local piece of the energy with respect to density)
  !! given in SOLER equation 10. The \(\text{u_alpha(k)}\) functions were found 
  !! while calculating the energy. They are passed in as the matrix \(\text{u_vdW}\).
  !! Most of the required derivatives were calculated in the \(\texttt{get_q0_on_grid}\)
  !! routine, but the derivative of the interpolation polynomials, \(P_\alpha(q)\), 
  !! (SOLER equation 3) with respect to q is interpolated here, along with the polynomials
  !! themselves.
  
  USE gvect,               ONLY : g
  USE cell_base,           ONLY : alat, tpiba

  IMPLICIT NONE

  REAL(DP), INTENT(IN) ::  q0(:), grad_rho(:,:)       ! Input arrays holding the value of q0 for
                                                      ! all points assigned to this processor and
                                                      ! the gradient of the charge density for
                                                      ! points assigned to this processor.

  REAL(DP), INTENT(IN) :: dq0_drho(:), dq0_dgradrho(:)! The derivative of q0 with respect to the
                                                      ! charge density and gradient of the charge
                                                      ! density (almost). See comments in the
                                                      ! get_q0_on_grid subroutine above.

  COMPLEX(DP), INTENT(IN) :: u_vdW(:,:)               ! The functions u_alpha(r) obtained by
                                                      ! inverse transforming the functions
                                                      ! u_alph(k). See equations 11 and 12 in SOLER.

  REAL(DP), INTENT(INOUT) :: potential(:)             ! The non-local correlation potential for
                                                      ! points on the grid over the whole cell (not
                                                      ! just those assigned to this processor).

  REAL(DP), ALLOCATABLE, SAVE :: d2y_dx2(:,:)         ! Second derivatives of P_alpha polynomials
                                                      ! for interpolation.

  INTEGER :: i_grid, P_i,icar                         ! Index variables.

  INTEGER :: q_low, q_hi, q                           ! Variables to find the bin in the q_mesh that
                                                      ! a particular q0 belongs to (for interpolation).
  REAL(DP) :: dq, a, b, c, d, e, f                    ! Intermediate variables used in the
                                                      ! interpolation of the polynomials.

  REAL(DP) :: y(Nqs), dP_dq0, P                       ! The y values for a given polynomial (all 0
                                                      ! exept for element i of P_i) The derivative
                                                      ! of P at a given q0 and the value of P at a
                                                      ! given q0. Both of these are interpolated
                                                      ! below.

  REAL(DP) :: gradient2                               ! Squared gradient.

  REAL(DP)   , ALLOCATABLE ::h_prefactor(:)
  COMPLEX(DP), ALLOCATABLE ::h(:)




  ALLOCATE ( h_prefactor(dfftp%nnr), h(dfftp%nnr) )

  potential     = 0.0D0
  h_prefactor   = 0.0D0


  ! --------------------------------------------------------------------
  ! Get the second derivatives of the P_i functions for interpolation.
  ! We have already calculated this once but it is very fast and it's
  ! just as easy to calculate it again.

  IF (.NOT. ALLOCATED( d2y_dx2) ) THEN

     ALLOCATE( d2y_dx2(Nqs, Nqs) )
     CALL initialize_spline_interpolation ( q_mesh, d2y_dx2(:,:) )

  end if


  DO i_grid = 1, dfftp%nnr

     q_low = 1
     q_hi = Nqs


     ! -----------------------------------------------------------------
     ! Figure out which bin our value of q0 is in in the q_mesh.

     DO WHILE ( (q_hi - q_low) > 1)

        q = INT((q_hi + q_low)/2)

        IF (q_mesh(q) > q0(i_grid)) THEN
           q_hi = q
        ELSE
           q_low = q
        END IF

     END DO

     IF ( q_hi == q_low ) CALL errore('vdW_DF_potential','qhi == qlow',1)

     dq = q_mesh(q_hi) - q_mesh(q_low)

     a = (q_mesh(q_hi) - q0(i_grid))/dq
     b = (q0(i_grid) - q_mesh(q_low))/dq
     c = (a**3 - a)*dq**2/6.0D0
     d = (b**3 - b)*dq**2/6.0D0
     e = (3.0D0*a**2 - 1.0D0)*dq/6.0D0
     f = (3.0D0*b**2 - 1.0D0)*dq/6.0D0

     DO P_i = 1, Nqs
        y = 0.0D0
        y(P_i) = 1.0D0

        P      = a*y(q_low) + b*y(q_hi)  + c*d2y_dx2(P_i,q_low) + d*d2y_dx2(P_i,q_hi)
        dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(P_i,q_low) + f*d2y_dx2(P_i,q_hi)


        ! --------------------------------------------------------------
        ! The first term in equation 10 of SOLER.
        potential(i_grid) = potential(i_grid) + u_vdW(i_grid,P_i)* (P + dP_dq0 * dq0_drho(i_grid))
        IF (q0(i_grid) .NE. q_mesh(Nqs)) THEN
           h_prefactor(i_grid) = h_prefactor(i_grid) + u_vdW(i_grid,P_i)*dP_dq0*dq0_dgradrho(i_grid)
        END IF

     END DO

  END DO

  DO icar = 1,3

     h(:) = CMPLX( h_prefactor(:) * grad_rho(icar,:), 0.0_DP, kind=dp )

     DO i_grid = 1, dfftp%nnr
        gradient2 = grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2
        IF ( gradient2 > 0.0D0 ) h(i_grid) = h(i_grid) / SQRT( gradient2 )
     END DO

     CALL fwfft ('Rho', h, dfftp)
     h(dfftp%nl(:)) = CMPLX(0.0_DP,1.0_DP, kind=dp) * tpiba * g(icar,:) * h(dfftp%nl(:))
     IF (gamma_only) h(dfftp%nlm(:)) = CONJG(h(dfftp%nl(:)))
     CALL invfft ('Rho', h, dfftp)
     potential(:) = potential(:) - REAL(h(:))

  END Do

  DEALLOCATE ( h_prefactor, h )

  END SUBROUTINE vdW_DF_potential








  ! ####################################################################
  !                       |                        |
  !                       |  SPLINE_INTERPOLATION  |
  !                       |________________________|
  !
  !
  SUBROUTINE spline_interpolation (x, evaluation_points, values)
  
  !! This routine is modeled after an algorithm from NUMERICAL_RECIPES It
  !! was adapted for Fortran, of course and for the problem at hand, in
  !! that it finds the bin a particular x value is in and then loops over
  !! all the \(P_i\) functions so we only have to find the bin once.

  IMPLICIT NONE

  REAL(DP), INTENT(IN) :: x(:), evaluation_points(:)     ! Input variables. The x values used to
                                                         ! form the interpolation (q_mesh in this
                                                         ! case) and the values of q0 for which we
                                                         ! are interpolating the function.

  COMPLEX(DP), INTENT(INOUT) :: values(:,:)              ! An output array (allocated outside this
                                                         ! routine) that stores the interpolated
                                                         ! values of the P_i (SOLER equation 3)
                                                         ! polynomials. The format is
                                                         ! values(grid_point, P_i).

  INTEGER :: Ngrid_points, Nx                            ! Total number of grid points to evaluate
                                                         ! and input x points.

  REAL(DP), ALLOCATABLE, SAVE :: d2y_dx2(:,:)            ! The second derivatives required to do
                                                         ! the interpolation.

  INTEGER :: i_grid, lower_bound, upper_bound, P_i       ! Some indexing variables.

  REAL(DP), ALLOCATABLE :: y(:)                          ! Temporary variables needed for the
  REAL(DP) :: a, b, c, d, dx                             ! interpolation.




  Nx = size(x)
  Ngrid_points = size(evaluation_points)


  ! --------------------------------------------------------------------
  ! Allocate the temporary array.

  ALLOCATE( y(Nx) )


  ! --------------------------------------------------------------------
  ! If this is the first time this routine has been called we need to
  ! get the second derivatives (d2y_dx2) required to perform the
  ! interpolations. So we allocate the array and call
  ! initialize_spline_interpolation to get d2y_dx2.

  IF (.NOT. ALLOCATED(d2y_dx2) ) THEN

     ALLOCATE( d2y_dx2(Nx,Nx) )
     CALL initialize_spline_interpolation(x, d2y_dx2)

  END IF

  DO i_grid=1, Ngrid_points

     lower_bound = 1
     upper_bound = Nx

     DO WHILE ( (upper_bound - lower_bound) > 1 )

        idx = (upper_bound+lower_bound) / 2

        IF ( evaluation_points(i_grid) > x(idx) ) THEN
           lower_bound = idx
        ELSE
           upper_bound = idx
        END IF

     END DO

     dx = x(upper_bound)-x(lower_bound)

     a = (x(upper_bound) - evaluation_points(i_grid))/dx
     b = (evaluation_points(i_grid) - x(lower_bound))/dx
     c = ((a**3-a)*dx**2)/6.0D0
     d = ((b**3-b)*dx**2)/6.0D0

     DO P_i = 1, Nx

        y = 0
        y(P_i) = 1

        values(i_grid, P_i) = a*y(lower_bound) + b*y(upper_bound) &
             + (c*d2y_dx2(P_i,lower_bound) + d*d2y_dx2(P_i, upper_bound))
     END DO

  END DO

  DEALLOCATE( y )

  END SUBROUTINE spline_interpolation








  ! ####################################################################
  !                  |                                   |
  !                  |  INITIALIZE_SPLINE_INTERPOLATION  |
  !                  |___________________________________|
  !
  
  SUBROUTINE initialize_spline_interpolation (x, d2y_dx2)
  
  !! This routine is modeled after an algorithm from NUMERICAL RECIPES It
  !! was adapted for Fortran and for the problem at hand.

  IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: x(:)               ! The input abscissa values.
  REAL(DP), INTENT(INOUT) :: d2y_dx2(:,:)       ! The output array (allocated outside this routine)
                                                ! that holds the second derivatives required for
                                                ! interpolating the function.

  INTEGER :: Nx, P_i                            ! The total number of x points and some indexing
                                                ! variables.

  REAL(DP), ALLOCATABLE :: temp_array(:), y(:)  ! Some temporary arrays required. y is the array
                                                ! that holds the funcion values (all either 0 or
                                                ! 1 here).

  REAL(DP) :: temp1, temp2                      ! Some temporary variables required.




  Nx = SIZE(x)

  ALLOCATE( temp_array(Nx), y(Nx) )

  DO P_i=1, Nx

     ! -----------------------------------------------------------------
     ! In the Soler method, the polynomicals that are interpolated are Kroneker
     ! delta funcions at a particular q point. So, we set all y values to 0
     ! except the one corresponding to the particular function P_i.

     y = 0.0D0
     y(P_i) = 1.0D0

     d2y_dx2(P_i,1) = 0.0D0
     temp_array(1) = 0.0D0

     DO idx = 2, Nx-1

        temp1 = (x(idx)-x(idx-1))/(x(idx+1)-x(idx-1))
        temp2 = temp1 * d2y_dx2(P_i,idx-1) + 2.0D0
        d2y_dx2(P_i,idx) = (temp1-1.0D0)/temp2

        temp_array(idx) = (y(idx+1)-y(idx))/(x(idx+1)-x(idx)) &
             - (y(idx)-y(idx-1))/(x(idx)-x(idx-1))
        temp_array(idx) = (6.0D0*temp_array(idx)/(x(idx+1)-x(idx-1)) &
             - temp1*temp_array(idx-1))/temp2

     END DO

     d2y_dx2(P_i,Nx) = 0.0D0

     DO idx=Nx-1, 1, -1

        d2y_dx2(P_i,idx) = d2y_dx2(P_i,idx) * d2y_dx2(P_i,idx+1) + temp_array(idx)

     END DO

  END DO

  DEALLOCATE( temp_array, y )

  END SUBROUTINE initialize_spline_interpolation








  ! ####################################################################
  !                          |                    |
  !                          | INTERPOLATE_KERNEL |
  !                          |____________________|
  !
  
  SUBROUTINE interpolate_kernel (k, kernel_of_k)

  !! This routine is modeled after an algorithm from NUMERICAL RECIPES
  !! Adapted for Fortran and the problem at hand. This function is used
  !! to find the Phi-alpha-beta needed for equations 8 and 11 of SOLER.
  
  IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: k                 ! Input value, the magnitude of the g-vector
                                               ! for the current point.

  REAL(DP), INTENT(INOUT) :: kernel_of_k(:,:)  ! An output array (allocated outside this routine)
                                               ! that holds the interpolated value of the kernel
                                               ! for each pair of q points (i.e. the phi_alpha_beta
                                               ! of the Soler method.

  INTEGER  :: q1_i, q2_i, k_i                  ! Indexing variables.

  REAL(DP) :: A, B, C, D                       ! Intermediate values for the interpolation.




  ! --------------------------------------------------------------------
  ! Check to make sure that the kernel table we have is capable of
  ! dealing with this value of k. If k is larger than
  ! Nr_points*2*pi/r_max then we can't perform the interpolation. In
  ! that case, a kernel file should be generated with a larger number of
  ! radial points.

  IF ( k >= Nr_points*dk ) THEN

     WRITE(*,'(A,F10.5,A,F10.5)') "k =  ", k, "     k_max =  ", Nr_points*dk
     CALL errore('interpolate kernel', 'k value requested is out of range',1)

  END IF

  kernel_of_k = 0.0D0


  ! --------------------------------------------------------------------
  ! This integer division figures out which bin k is in since the kernel
  ! is set on a uniform grid.

  k_i = INT(k/dk)


  ! --------------------------------------------------------------------
  ! Test to see if we are trying to interpolate a k that is one of the
  ! actual function points we have. The value is just the value of the
  ! function in that case.

  IF ( MOD(k,dk) == 0 ) THEN

     DO q1_i = 1, Nqs
        DO q2_i = 1, q1_i

           kernel_of_k(q1_i, q2_i) = kernel(k_i,q1_i, q2_i)
           kernel_of_k(q2_i, q1_i) = kernel(k_i,q2_i, q1_i)

        END DO
     END DO

     RETURN

  END IF


  ! --------------------------------------------------------------------
  ! If we are not on a function point then we carry out the
  ! interpolation.

  A = (dk*(k_i+1.0D0) - k)/dk
  B = (k - dk*k_i)/dk
  C = (A**3-A)*dk**2/6.0D0
  D = (B**3-B)*dk**2/6.0D0

  DO q1_i = 1, Nqs
     DO q2_i = 1, q1_i

        kernel_of_k(q1_i, q2_i) = A*kernel(k_i, q1_i, q2_i) + B*kernel(k_i+1, q1_i, q2_i) &
                        + (C*d2phi_dk2(k_i, q1_i, q2_i) + D*d2phi_dk2(k_i+1, q1_i, q2_i))

        kernel_of_k(q2_i, q1_i) = kernel_of_k(q1_i, q2_i)

     END DO
  END DO

  END SUBROUTINE interpolate_kernel








  ! ####################################################################
  !                         |                 |
  !                         |  VDW_DF_STRESS  |
  !                         |_________________|

  SUBROUTINE vdW_DF_stress (rho_valence, rho_core, nspin, sigma) ! PH adjusted for wrapper spin/nospin
  
  !! vdW-DF stress calculation.
  
  use gvect,           ONLY : ngm, g
  USE cell_base,       ONLY : tpiba

  IMPLICIT NONE

  REAL(DP), INTENT(IN)     :: rho_valence(:,:)       !
  REAL(dp), INTENT(IN)     :: rho_core(:)            ! Input variables.
  INTEGER,  INTENT(IN)     :: nspin                  !
  REAL(dp), INTENT(INOUT)  :: sigma(3,3)             !

  REAL(DP), ALLOCATABLE    :: grad_rho(:,:)          !
  REAL(DP), ALLOCATABLE    :: total_rho(:)           ! Rho values.

  real(dp), allocatable    :: rho_up(:)              !
  real(dp), allocatable    :: rho_down(:)            !

  real(dp), allocatable :: grad_rho_up(:,:)          ! The gradient of the up charge density.
                                                     ! Same format as grad_rho
  real(dp), allocatable :: grad_rho_down(:,:)        ! The gradient of the down charge density.
                                                     ! Same format as grad_rho



  REAL(DP), ALLOCATABLE    :: q0(:)                  !
  REAL(DP), ALLOCATABLE    :: dq0_drho(:)            ! q-values.
  REAL(DP), ALLOCATABLE    :: dq0_dgradrho(:)        !
  real(dp), allocatable    :: dq0_drho_up(:)         ! The derivative of the saturated q0
  real(dp), allocatable    :: dq0_drho_down(:)       ! with respect to the spin charge density 
  real(dp), allocatable    :: dq0_dgradrho_up(:)     ! The derivative of the saturated q0
  real(dp), allocatable    :: dq0_dgradrho_down(:)   ! with respect to the gradient of the spin charge sensity

  COMPLEX(DP), ALLOCATABLE :: thetas(:,:)            ! Thetas.
  INTEGER                  :: i_proc, theta_i, l, m

  REAL(DP)                 :: sigma_grad(3,3)
  REAL(DP)                 :: sigma_ker(3,3)




  ! --------------------------------------------------------------------
  ! Tests
 
  IF ( inlc > 5 ) CALL errore( 'xc_vdW_DF', 'inlc not implemented', 1 )

#if defined (__SPIN_BALANCED)
  IF ( nspin==2 ) THEN
     WRITE(stdout,'(/,/ "     Performing spin-balanced Ecnl stress calculation!")')
  ELSE IF ( nspin > 2 ) THEN
     CALL errore ('vdW_DF_stress', 'noncollinear vdW stress not implemented', 1)
  END IF
#else
  IF ( nspin>2 ) THEN
     CALL errore ('vdW_DF_stress', 'vdW stress not implemented for nspin > 2', 1)
  END IF
#endif

  sigma(:,:)      = 0.0_DP
  sigma_grad(:,:) = 0.0_DP
  sigma_ker(:,:)  = 0.0_DP


  ! --------------------------------------------------------------------
  ! Allocations

  ALLOCATE( total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr), thetas(dfftp%nnr, Nqs), q0(dfftp%nnr) )
  ALLOCATE ( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
#if defined (__SPIN_BALANCED)
#else
  IF (nspin==2) THEN
     allocate( rho_up(dfftp%nnr), rho_down(dfftp%nnr) )
     allocate( grad_rho_up(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr) )
     allocate( dq0_drho_up (dfftp%nnr), dq0_dgradrho_up  (dfftp%nnr) )
     allocate( dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )
  ENDIF
#endif

  ! --------------------------------------------------------------------
  ! Charge

  total_rho = rho_valence(:,1) + rho_core(:)
#if defined (__SPIN_BALANCED)
#else
  IF (nspin==2) THEN
    rho_up    = ( rho_valence(:,1) + rho_valence(:,2) + rho_core(:) )*0.5D0
    rho_down  = ( rho_valence(:,1) - rho_valence(:,2) + rho_core(:) )*0.5D0
  ENDIF
#endif



  ! --------------------------------------------------------------------
  ! Here we calculate the gradient in reciprocal space using FFT.

  CALL fft_gradient_r2r (dfftp, total_rho,  g, grad_rho)
#if defined (__SPIN_BALANCED)
#else
  IF (nspin==2) THEN
     call fft_gradient_r2r (dfftp, rho_up,    g, grad_rho_up)
     call fft_gradient_r2r (dfftp, rho_down,  g, grad_rho_down)
  ENDIF
#endif

  ! --------------------------------------------------------------------
  ! Get q0.

#if defined (__SPIN_BALANCED)
  CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
#else
  IF (nspin == 1) THEN
     CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
  ELSEIF (nspin==2) THEN
     CALL get_q0_on_grid_spin ( total_rho, rho_up, rho_down, grad_rho, grad_rho_up, grad_rho_down, &
          q0, dq0_drho_up, dq0_drho_down, dq0_dgradrho_up, dq0_dgradrho_down, thetas)
  ENDIF
#endif

  ! --------------------------------------------------------------------
  ! Stress

#if defined (__SPIN_BALANCED)
     CALL vdW_DF_stress_gradient (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, &
                                  thetas, sigma_grad)
#else
  IF (nspin == 1) THEN
     CALL vdW_DF_stress_gradient (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, &
                                  thetas, sigma_grad)
  ELSEIF (nspin == 2) THEN
     CALL vdW_DF_stress_gradient_spin (total_rho, grad_rho_up, grad_rho_down, q0, &
                                       dq0_dgradrho_up, dq0_dgradrho_down, &
                                       thetas, sigma_grad)
  ENDIF
#endif


  CALL vdW_DF_stress_kernel   (total_rho, q0, thetas, sigma_ker)

  sigma = - (sigma_grad + sigma_ker)

  DO l = 1, 3
     DO m = 1, l - 1
        sigma (m, l) = sigma (l, m)
     END DO
  END DO

  DEALLOCATE( total_rho, grad_rho, thetas, q0 )
  DEALLOCATE( dq0_drho, dq0_dgradrho )
#if defined (__SPIN_BALANCED)
#else
  IF (nspin == 2) THEN
     deallocate( rho_up, rho_down )
     deallocate( grad_rho_up, grad_rho_down )
     deallocate( dq0_drho_up, dq0_drho_down, dq0_dgradrho_up, dq0_dgradrho_down )
  ENDIF
#endif

  END SUBROUTINE vdW_DF_stress ! PH adjusted for wrapper spin/nospin


! -------------------------------------------------------------------------
 ! Begin Spin vdW-DF_strees_gradient implemented Per Hyldgaard 2019, GPL. No Waranties
 ! Adapted from the original nspin = 1 code (subroutine below) by Thonhauser and coauthors
! -------------------------------------------------------------------------

  ! ####################################################################
  !                     |                               |
  !                     |  VDW_DF_STRESS_GRADIENT_SPIN  |
  !                     |_______________________________|

  SUBROUTINE vdW_DF_stress_gradient_spin (total_rho, grad_rho_up, grad_rho_down, q0, &
                                          dq0_dgradrho_up, dq0_dgradrho_down, &
                                          thetas, sigma)

  USE gvect,                 ONLY : ngm, g, gg, gstart
  USE cell_base,             ONLY : omega, tpiba, alat, at, tpiba2

  implicit none

  real(dp), intent(IN)     :: total_rho(:)           !
  real(dp), intent(IN)     :: grad_rho_up (:, :)     ! Input variables.
  real(dp), intent(IN)     :: grad_rho_down(:, :)    !
  real(dp), intent(inout)  :: sigma(:,:)             !
  real(dp), intent(IN)     :: q0(:)                  !
  real(dp), intent(IN)     :: dq0_dgradrho_up(:)     !
  real(dp), intent(IN)     :: dq0_dgradrho_down(:)   !
  complex(dp), intent(IN)  :: thetas(:,:)            !

  complex(dp), allocatable :: u_vdW(:,:)             !

  real(dp), allocatable    :: d2y_dx2(:,:)           !
  real(dp) :: y(Nqs), dP_dq0, P, a, b, c, d, e, f    ! Interpolation.
  real(dp) :: dq                                     !

  integer  :: q_low, q_hi, q, q1_i, q2_i , g_i       ! Loop and q-points.

  integer  :: l, m
  real(dp) :: prefactor_up, prefactor_down           ! Final summation of sigma.
  real(dp) :: grad2_up, grad2_down                   ! Magnitude of density gradient.


  integer  :: i_proc, theta_i, i_grid, q_i, &        !
              ix, iy, iz                             ! Iterators.

  character(LEN=1) :: intvar

  allocate( d2y_dx2(Nqs, Nqs) )
  allocate( u_vdW(dfftp%nnr, Nqs) )

  sigma(:,:) = 0.0_DP
  prefactor_up  = 0.0_DP
  prefactor_down  = 0.0_DP

  ! --------------------------------------------------------------------
  ! Get u in k-space.

  call thetas_to_uk(thetas, u_vdW)

  ! --------------------------------------------------------------------
  ! Get u in real space.

  do theta_i = 1, Nqs
     CALL invfft('Rho', u_vdW(:,theta_i), dfftp)
  end do


  ! --------------------------------------------------------------------
  ! Get the second derivatives for interpolating the P_i.

  call initialize_spline_interpolation(q_mesh, d2y_dx2(:,:))


  ! --------------------------------------------------------------------
  ! Do the real space integration to obtain the stress component.

  do i_grid = 1, dfftp%nnr

     if ( total_rho(i_grid) < epsr ) cycle

     q_low = 1
     q_hi  = Nqs
     grad2_up = sqrt( grad_rho_up(1,i_grid)**2 &
                + grad_rho_up(2,i_grid)**2 + grad_rho_up(3,i_grid)**2 )
     grad2_down = sqrt( grad_rho_down(1,i_grid)**2 &
                  + grad_rho_down(2,i_grid)**2 + grad_rho_down(3,i_grid)**2 )

     if ( grad2_up == 0.0_dp .or. grad2_down == 0.0_dp) cycle

     ! -----------------------------------------------------------------
     ! Figure out which bin our value of q0 is in the q_mesh.

     do while ( (q_hi - q_low) > 1)

        q = int((q_hi + q_low)/2)

        if (q_mesh(q) > q0(i_grid)) then
            q_hi = q
        else
            q_low = q
        end if

     end do

     if (q_hi == q_low) call errore('vdW_DF_stress_gradient_spin','qhi == qlow',1)

     dq = q_mesh(q_hi) - q_mesh(q_low)

     a = (q_mesh(q_hi) - q0(i_grid))/dq
     b = (q0(i_grid) - q_mesh(q_low))/dq
     c = (a**3 - a)*dq**2/6.0D0
     d = (b**3 - b)*dq**2/6.0D0
     e = (3.0D0*a**2 - 1.0D0)*dq/6.0D0
     f = (3.0D0*b**2 - 1.0D0)*dq/6.0D0

    do q_i = 1, Nqs

        y(:)   = 0.0D0
        y(q_i) = 1.0D0

        dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(q_i,q_low) + f*d2y_dx2(q_i,q_hi)

        prefactor_up = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho_up(i_grid) / grad2_up
        prefactor_down = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho_down(i_grid) / grad2_down

        do l = 1, 3
            do m = 1, l

                sigma (l, m) = sigma (l, m) -  e2 * prefactor_up * &
                               (grad_rho_up(l,i_grid) * grad_rho_up(m,i_grid))
                sigma (l, m) = sigma (l, m) -  e2 * prefactor_down * &
                               (grad_rho_down(l,i_grid) * grad_rho_down(m,i_grid))
            end do
        end do

     end do

  end do

  call mp_sum(  sigma, intra_bgrp_comm )

  call dscal (9, 1.d0 / (dfftp%nr1 * dfftp%nr2 * dfftp%nr3), sigma, 1)

  deallocate( d2y_dx2, u_vdW )

  END SUBROUTINE vdW_DF_stress_gradient_spin


  ! ####################################################################
  !                     |                          |
  !                     |  VDW_DF_STRESS_GRADIENT  |
  !                     |__________________________|

  SUBROUTINE vdW_DF_stress_gradient (total_rho, grad_rho, q0, &
             dq0_drho, dq0_dgradrho, thetas, sigma)

  USE gvect,                 ONLY : ngm, g, gstart
  USE cell_base,             ONLY : omega, tpiba, alat, at, tpiba2

  IMPLICIT NONE

  REAL(DP), INTENT(IN)     :: total_rho(:)           !
  REAL(DP), INTENT(IN)     :: grad_rho(:, :)         !
  REAL(DP), INTENT(INOUT)  :: sigma(:,:)             !
  REAL(DP), INTENT(IN)     :: q0(:)                  ! Input variables.
  REAL(DP), INTENT(IN)     :: dq0_drho(:)            !
  REAL(DP), INTENT(IN)     :: dq0_dgradrho(:)        !
  COMPLEX(DP), INTENT(IN)  :: thetas(:,:)            !

  COMPLEX(DP), ALLOCATABLE :: u_vdW(:,:)             !

  REAL(DP), ALLOCATABLE    :: d2y_dx2(:,:)           !
  REAL(DP) :: y(Nqs), dP_dq0, P, a, b, c, d, e, f    ! Interpolation.
  REAL(DP) :: dq                                     !

  INTEGER  :: q_low, q_hi, q, q1_i, q2_i , g_i       ! Loop and q-points.

  INTEGER  :: l, m
  REAL(DP) :: prefactor                              ! Final summation of sigma
  REAL(DP) :: grad2                                  ! magnitude of density gradient.

  INTEGER  :: i_proc, theta_i, i_grid, q_i, &        ! Iterators.
              ix, iy, iz                             !

  CHARACTER(LEN=1) :: intvar




  ALLOCATE( d2y_dx2(Nqs, Nqs) )
  ALLOCATE( u_vdW(dfftp%nnr, Nqs) )

  sigma(:,:) = 0.0_DP
  prefactor  = 0.0_DP


  ! --------------------------------------------------------------------
  ! Get u in k-space.

  CALL thetas_to_uk(thetas, u_vdW)


  ! --------------------------------------------------------------------
  ! Get u in real space.

  DO theta_i = 1, Nqs
     CALL invfft('Rho', u_vdW(:,theta_i), dfftp)
  END DO


  ! --------------------------------------------------------------------
  ! Get the second derivatives for interpolating the P_i.

  CALL initialize_spline_interpolation(q_mesh, d2y_dx2(:,:))


  ! --------------------------------------------------------------------
  ! Do the real space integration to obtain the stress component.

  DO i_grid = 1, dfftp%nnr

     IF ( total_rho(i_grid) < epsr ) CYCLE

     q_low = 1
     q_hi  = Nqs
     grad2 = sqrt( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 )

     IF ( grad2 == 0.0_dp ) CYCLE

     ! -----------------------------------------------------------------
     ! Figure out which bin our value of q0 is in the q_mesh.
     DO WHILE ( (q_hi - q_low) > 1)

        q = INT((q_hi + q_low)/2)

        IF (q_mesh(q) > q0(i_grid)) THEN
            q_hi = q
        ELSE
            q_low = q
        END IF

     END DO

     IF (q_hi == q_low) call errore('stress_vdW_gradient','qhi == qlow', 1)

     dq = q_mesh(q_hi) - q_mesh(q_low)
     a  = (q_mesh(q_hi) - q0(i_grid))/dq
     b  = (q0(i_grid) - q_mesh(q_low))/dq
     c  = (a**3 - a)*dq**2/6.0D0
     d  = (b**3 - b)*dq**2/6.0D0
     e  = (3.0D0*a**2 - 1.0D0)*dq/6.0D0
     f  = (3.0D0*b**2 - 1.0D0)*dq/6.0D0

     DO q_i = 1, Nqs

        y(:)   = 0.0D0
        y(q_i) = 1.0D0

        dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(q_i,q_low) + f*d2y_dx2(q_i,q_hi)

        prefactor = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho(i_grid) / grad2

        DO l = 1, 3
        DO m = 1, l

            sigma (l, m) = sigma (l, m) -  e2 * prefactor * &
                           (grad_rho(l,i_grid) * grad_rho(m,i_grid))
        END DO
        END DO

     END DO

  END DO

  CALL mp_sum(  sigma, intra_bgrp_comm )

  CALL dscal (9, 1.0D0 / (dfftp%nr1 * dfftp%nr2 * dfftp%nr3), sigma, 1)

  DEALLOCATE( d2y_dx2, u_vdW )

  END SUBROUTINE vdW_DF_stress_gradient








  ! ####################################################################
  !                      |                        |
  !                      |  VDW_DF_STRESS_KERNEL  |
  !                      |________________________|

  SUBROUTINE vdW_DF_stress_kernel (total_rho, q0, thetas, sigma)

  USE gvect,                 ONLY : ngm, g, gg, igtongl, gl, ngl, gstart
  USE cell_base,             ONLY : omega, tpiba, tpiba2

  IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: q0(:)
  REAL(DP), INTENT(IN)    :: total_rho(:)
  REAL(DP), INTENT(INOUT) :: sigma(3,3)
  COMPLEX(DP), INTENT(IN) :: thetas(:,:)

  REAL(DP), ALLOCATABLE   :: dkernel_of_dk(:,:)

  INTEGER                 :: l, m, q1_i, q2_i , g_i
  INTEGER                 :: last_g, theta_i
  REAL(DP)                :: g2, ngmod2, g_kernel, G_multiplier




  ALLOCATE( dkernel_of_dk(Nqs, Nqs) )

  sigma(:,:) = 0.0_DP

  ! --------------------------------------------------------------------
  ! Integration in g-space.

  last_g = -1

  G_multiplier = 1.0D0

  IF ( gamma_only ) G_multiplier = 2.0D0

  DO g_i = gstart, ngm

     g2 = gg (g_i) * tpiba2
     g_kernel = SQRT(g2)

     IF ( igtongl(g_i) .NE. last_g) THEN

        CALL interpolate_Dkernel_Dk(g_kernel, dkernel_of_dk)  ! Gets the derivatives.
        last_g = igtongl(g_i)

     END IF

     DO q2_i = 1, Nqs
     DO q1_i = 1, Nqs
        DO l = 1, 3
        DO m = 1, l

           sigma (l, m) = sigma (l, m) - G_multiplier * 0.5 * e2 * thetas(dfftp%nl(g_i),q1_i) * &
                          dkernel_of_dk(q1_i,q2_i)*conjg(thetas(dfftp%nl(g_i),q2_i))* &
                          (g (l, g_i) * g (m, g_i) * tpiba2) / g_kernel
        END DO
        END DO
     END DO
     END DO

     IF ( g_i < gstart ) sigma(:,:) = sigma(:,:) / G_multiplier

  END DO

  CALL mp_sum( sigma, intra_bgrp_comm )

  DEALLOCATE( dkernel_of_dk )

  END SUBROUTINE vdW_DF_stress_kernel








  ! ####################################################################
  !                        |                        |
  !                        | INTERPOLATE_DKERNEL_DK |
  !                        |________________________|

  SUBROUTINE interpolate_Dkernel_Dk (k, dkernel_of_dk)

  IMPLICIT NONE

  REAL(DP), INTENT(IN)    :: k                      ! Input value, the magnitude of the g-vector
                                                    ! for the current point.

  REAL(DP), INTENT(INOUT) :: dkernel_of_dk(Nqs,Nqs) ! An output array (allocated outside this
                                                    ! routine) that holds the interpolated value of
                                                    ! the kernel for each pair of q points (i.e. the
                                                    ! phi_alpha_beta of the Soler method.

  INTEGER :: q1_i, q2_i, k_i                        ! Indexing variables.

  REAL(DP) :: A, B, dAdk, dBdk, dCdk, dDdk          ! Intermediate values for the interpolation.




  IF ( k >= Nr_points*dk ) THEN

     WRITE(*,'(A,F10.5,A,F10.5)') "k =  ", k, "     k_max =  ", Nr_points*dk
     CALL errore('interpolate kernel', 'k value requested is out of range',1)

  END IF

  dkernel_of_dk = 0.0D0

  k_i  = INT(k/dk)

  A    = (dk*(k_i+1.0D0) - k)/dk
  B    = (k - dk*k_i)/dk

  dAdk = -1.0D0/dk
  dBdk = 1.0D0/dk
  dCdk = -((3*A**2 -1.0D0)/6.0D0)*dk
  dDdk = ((3*B**2 -1.0D0)/6.0D0)*dk

  DO q1_i = 1, Nqs
     DO q2_i = 1, q1_i

        dkernel_of_dk(q1_i, q2_i) = dAdk*kernel(k_i, q1_i, q2_i) + dBdk*kernel(k_i+1, q1_i, q2_i) &
                            + dCdk*d2phi_dk2(k_i, q1_i, q2_i) + dDdk*d2phi_dk2(k_i+1, q1_i, q2_i)

        dkernel_of_dk(q2_i, q1_i) = dkernel_of_dk(q1_i, q2_i)

     END DO
  END DO

  END SUBROUTINE interpolate_Dkernel_Dk








  ! ####################################################################
  !                          |              |
  !                          | thetas_to_uk |
  !                          |______________|


  SUBROUTINE thetas_to_uk (thetas, u_vdW)

  USE gvect,           ONLY : gg, ngm, igtongl, gl, ngl, gstart
  USE cell_base,       ONLY : tpiba, omega

  IMPLICIT NONE

  COMPLEX(DP), INTENT(IN)  :: thetas(:,:)       ! On input this variable holds the theta functions
                                                ! (equation 8, SOLER) in the format
                                                ! thetas(grid_point, theta_i).
  COMPLEX(DP), INTENT(OUT) :: u_vdW(:,:)        ! On output this array holds u_alpha(k) =
                                                ! Sum_j[theta_beta(k)phi_alpha_beta(k)].

  REAL(DP), ALLOCATABLE    :: kernel_of_k(:,:)  ! This array will hold the interpolated kernel
                                                ! values for each pair of q values in the q_mesh.

  REAL(DP) :: g
  INTEGER  :: last_g, g_i, q1_i, q2_i, i_grid   ! Index variables.

  COMPLEX(DP) :: theta(Nqs)                     ! Temporary storage vector used since we are
                                                ! overwriting the thetas array here.




  ALLOCATE( kernel_of_k(Nqs, Nqs) )

  u_vdW(:,:) = CMPLX(0.0_DP, 0.0_DP, kind=dp)

  last_g = -1

  DO g_i = 1, ngm

     IF ( igtongl(g_i) .ne. last_g) THEN

        g = SQRT(gl(igtongl(g_i))) * tpiba
        CALL interpolate_kernel(g, kernel_of_k)
        last_g = igtongl(g_i)

     END IF

     theta = thetas(dfftp%nl(g_i),:)

     DO q2_i = 1, Nqs
        DO q1_i = 1, Nqs
           u_vdW(dfftp%nl(g_i),q2_i) = u_vdW(dfftp%nl(g_i),q2_i) + kernel_of_k(q2_i,q1_i)*theta(q1_i)
        END DO
     END DO

  END Do

  IF ( gamma_only ) u_vdW(dfftp%nlm(:),:) = CONJG(u_vdW(dfftp%nl(:),:))

  DEALLOCATE( kernel_of_k )

  END SUBROUTINE thetas_to_uk








  ! ####################################################################
  !                           |                 |
  !                           | GENERATE_KERNEL |
  !                           |_________________|
  !
  ! The original definition of the kernel function is given in DION
  ! equations 14-16. The Soler method makes the kernel function a
  ! function of only 1 variable (r) by first putting it in the form
  ! phi(q1*r, q2*r). Then, the q-dependence is removed by expanding the
  ! function in a special way (see SOLER equation 3). This yields a
  ! separate function for each pair of q points that is a function of r
  ! alone. There are (Nqs^2+Nqs)/2 unique functions, where Nqs is the
  ! number of q points used. In the Soler method, the kernel is first
  ! made in the form phi(d1, d2) but this is not done here. It was found
  ! that, with q's chosen judiciously ahead of time, the kernel and the
  ! second derivatives required for interpolation could be tabulated
  ! ahead of time for faster use of the vdW-DF functional. Through
  ! testing we found no need to soften the kernel and correct for this
  ! later (see SOLER eqations 6-7).
  !
  ! The algorithm employed here is "embarrassingly parallel," meaning
  ! that it parallelizes very well up to (Nqs^2+Nqs)/2 processors,
  ! where, again, Nqs is the number of q points chosen. However,
  ! parallelization on this scale is unnecessary. In testing the code
  ! runs in under a minute on 16 Intel Xeon processors.
  !
  ! IMPORTANT NOTICE: Results are very sensitive to compilation details.
  ! In particular, the usage of FMA (Fused Multiply-and-Add)
  ! instructions used by modern CPUs such as AMD Interlagos (Bulldozer)
  ! and Intel Ivy Bridge may affect quite heavily some components of the
  ! kernel (communication by Ake Sandberg, Umea University). In practice
  ! this should not be a problem, since most affected elements are the
  ! less relevant ones.
  !
  ! Some of the algorithms here are somewhat modified versions of those
  ! found in:
  !
  !    Numerical Recipes in C; William H. Press, Brian P. Flannery, Saul
  !    A. Teukolsky, and William T. Vetterling. Cambridge University
  !    Press (1988).
  !
  ! hereafter referred to as NUMERICAL_RECIPES. The routines were
  ! translated to Fortran, of course and variable names are generally
  ! different.
  !
  ! For the calculation of the kernel we have benefited from access to
  ! earlier vdW-DF implementation into PWscf and ABINIT, written by Timo
  ! Thonhauser, Valentino Cooper, and David Langreth. These codes, in
  ! turn, benefited from earlier codes written by Maxime Dion and Henrik
  ! Rydberg.

  SUBROUTINE generate_kernel

  IMPLICIT NONE

  INTEGER  :: a_i, b_i, q1_i, q2_i, r_i
  ! Indexing variables.

  REAL(DP) :: weights( Nintegration_points )
  ! Array to hold dx values for the Gaussian-Legendre integration of the kernel.

  REAL(DP) :: sin_a( Nintegration_points ), cos_a( Nintegration_points )
  ! Sine and cosine values of the aforementioned points a.

  REAL(DP) :: d1, d2, d, integral
  ! Intermediate values.

  ! --------------------------------------------------------------------
  ! The following variables control the parallel environment.

  INTEGER :: my_start_q, my_end_q, Ntotal
  ! Starting and ending q value for each  processor, also the total
  ! number of calculations to do, i.e. (Nqs^2 + Nqs)/2.

  REAL(DP), ALLOCATABLE :: phi(:,:), phi_deriv(:,:)
  ! Arrays to store the kernel functions and their second derivatives.
  ! They are stored as phi(radial_point, idx).

  INTEGER, ALLOCATABLE  :: indices(:,:), proc_indices(:,:)
  ! Indices holds the values of q1 and q2 as partitioned out to the
  ! processors. It is an Ntotal x 2 array stored as indices(index of
  ! point number, q1:q2). Proc_indices holds the section of the indices
  ! array that is assigned to each processor. This is a Nproc x 2
  ! array, stored as proc_indices(processor_number,
  ! starting_index:ending_index)

  INTEGER :: Nper, Nextra, start_q, end_q
  ! Baseline number of jobs per processor, number of processors that
  ! get an extra job in case the number of jobs doesn't split evenly
  ! over the number of processors, starting index into the indices
  ! array, ending index into the indices array.

  INTEGER :: nproc, mpime
  ! Number or procs, rank of current processor.

  INTEGER :: proc_i, my_Nqs




  ! --------------------------------------------------------------------
  ! Start the timer.

  CALL start_clock ( 'vdW_kernel' )


  ! --------------------------------------------------------------------
  ! The total number of phi_alpha_beta functions that have to be
  ! calculated.

  Ntotal = (Nqs**2 + Nqs)/2
  ALLOCATE ( indices(Ntotal, 2) )


  ! --------------------------------------------------------------------
  ! This part fills in the indices array. It just loops through the q1
  ! and q2 values and stores them. Sections of this array will be
  ! assigned to each of the processors later.

  idx = 1

  DO q1_i = 1, Nqs
     DO q2_i = 1, q1_i
        indices(idx, 1) = q1_i
        indices(idx, 2) = q2_i
        idx = idx + 1
     END DO
  END DO


  ! --------------------------------------------------------------------
  ! Figure out the baseline number of functions to be calculated by each
  ! processor and how many processors get one extra job.

  nproc  = mp_size( intra_image_comm )
  mpime  = mp_rank( intra_image_comm )
  Nper   = Ntotal/nproc
  Nextra = MOD(Ntotal, nproc)

  ALLOCATE( proc_indices(nproc, 2) )

  start_q = 0
  end_q   = 0


  ! --------------------------------------------------------------------
  ! Loop over all the processors and figure out which section of the
  ! indices array each processor should do. All processors figure this
  ! out for every processor so there is no need to communicate results.

  DO proc_i = 1, nproc

     start_q = end_q + 1
     end_q   = start_q + (Nper - 1)
     IF (proc_i <= Nextra) end_q = end_q + 1

     ! This is to prevent trouble if number of processors exceeds Ntotal.
     IF ( proc_i > Ntotal ) THEN
        start_q    = Ntotal
        end_q      = Ntotal
     END IF

     IF ( proc_i == (mpime+1) ) THEN
        my_start_q = start_q
        my_end_q   = end_q
     END IF

     proc_indices(proc_i, 1) = start_q
     proc_indices(proc_i, 2) = end_q

  END DO


  ! --------------------------------------------------------------------
  ! Store how many jobs are assigned to me.

  my_Nqs    = my_end_q - my_start_q + 1
  ALLOCATE( phi( 0:Nr_points, my_Nqs ), phi_deriv( 0:Nr_points, my_Nqs ) )

  phi       = 0.0D0
  phi_deriv = 0.0D0
  kernel    = 0.0D0
  d2phi_dk2 = 0.0D0


  ! --------------------------------------------------------------------
  ! Find the integration points we are going to use in the
  ! Gaussian-Legendre integration.

  CALL prep_gaussian_quadrature( weights )


  ! --------------------------------------------------------------------
  ! Get a, a^2, sin(a), cos(a) and the weights for the Gaussian-Legendre
  ! integration.

  DO a_i=1, Nintegration_points
     a_points (a_i) = TAN( a_points(a_i) )
     a_points2(a_i) = a_points(a_i)**2
     weights(a_i)   = weights(a_i)*( 1 + a_points2(a_i) )
     cos_a(a_i)     = COS( a_points(a_i) )
     sin_a(a_i)     = SIN( a_points(a_i) )
  END DO


  ! --------------------------------------------------------------------
  ! Calculate the value of the W function defined in DION equation 16
  ! for each value of a and b.

  DO a_i = 1, Nintegration_points
  DO b_i = 1, Nintegration_points
     W_ab(a_i, b_i) = 2.0D0 * weights(a_i)*weights(b_i) * (           &
        (3.0D0-a_points2(a_i))*a_points(b_i) *sin_a(a_i)*cos_a(b_i) + &
        (3.0D0-a_points2(b_i))*a_points(a_i) *cos_a(a_i)*sin_a(b_i) + &
        (a_points2(a_i)+a_points2(b_i)-3.0D0)*sin_a(a_i)*sin_a(b_i) - &
        3.0D0*a_points(a_i)*a_points(b_i)*cos_a(a_i)*cos_a(b_i) )   / &
        (a_points(a_i)*a_points(b_i) )
  END DO
  END DO


  ! --------------------------------------------------------------------
  ! vdW-DF analysis tool as described in PRB 97, 085115 (2018).

  IF      ( vdW_DF_analysis == 1 ) THEN

     DO a_i = 1, Nintegration_points
     DO b_i = 1, Nintegration_points
        W_ab(a_i, b_i) = weights(a_i)*weights(b_i) *                  &
           a_points(a_i)*a_points(b_i)*sin_a(a_i)*sin_a(b_i)
     END DO
     END DO

  ELSE IF ( vdW_DF_analysis == 2 ) THEN

     DO a_i = 1, Nintegration_points
     DO b_i = 1, Nintegration_points
        W_ab(a_i, b_i) = W_ab(a_i, b_i) - weights(a_i)*weights(b_i) *  &
           a_points(a_i)*a_points(b_i)*sin_a(a_i)*sin_a(b_i)
     END DO
     END DO

  END IF


  ! --------------------------------------------------------------------
  ! Now, we loop over all the pairs (q1,q2) that are assigned to us and
  ! perform our calculations.

  DO idx = 1, my_Nqs

     ! -----------------------------------------------------------------
     ! First, get the value of phi(q1*r, q2*r) for each r and the
     ! particular values of q1 and q2 we are using.

     DO r_i = 1, Nr_points
        d1  = q_mesh( indices(idx+my_start_q-1, 1) ) * dr * r_i
        d2  = q_mesh( indices(idx+my_start_q-1, 2) ) * dr * r_i
        phi(r_i, idx) = phi_value(d1, d2)
     END DO


     ! -----------------------------------------------------------------
     ! Now, perform a radial FFT to turn our phi_alpha_beta(r) into
     ! phi_alpha_beta(k) needed for SOLER equation 8.

     CALL radial_fft( phi(:,idx) )


     ! -----------------------------------------------------------------
     ! Determine the spline interpolation coefficients for the Fourier
     ! transformed kernel function.

     CALL set_up_splines( phi(:, idx), phi_deriv(:, idx) )

  END DO


  ! --------------------------------------------------------------------
  ! Finally, we collect the results after letting everybody catch up.

  CALL mp_barrier( intra_image_comm )

  DO proc_i = 0, nproc-1

     IF ( proc_i >= Ntotal ) EXIT

     CALL mp_get ( phi      , phi      , mpime, 0, proc_i, 0, intra_image_comm )
     CALL mp_get ( phi_deriv, phi_deriv, mpime, 0, proc_i, 0, intra_image_comm )

     IF ( mpime == 0 ) THEN

        DO idx = proc_indices(proc_i+1,1), proc_indices(proc_i+1,2)
           q1_i = indices(idx, 1)
           q2_i = indices(idx, 2)
           kernel    (:, q1_i, q2_i) = phi       (:, idx - proc_indices(proc_i+1,1) + 1)
           d2phi_dk2 (:, q1_i, q2_i) = phi_deriv (:, idx - proc_indices(proc_i+1,1) + 1)
           kernel    (:, q2_i, q1_i) = kernel    (:, q1_i, q2_i)
           d2phi_dk2 (:, q2_i, q1_i) = d2phi_dk2 (:, q1_i, q2_i)
        END DO

     END IF

  END DO

  CALL mp_bcast ( kernel   , 0, intra_image_comm )
  CALL mp_bcast ( d2phi_dk2, 0, intra_image_comm )


  ! --------------------------------------------------------------------
  ! Keep the lines below for testing and combatibility with the old
  ! kernel file reading/writing method.
  !
  ! Writing the calculated kernel.
  !
  ! IF ( ionode ) THEN
  !    WRITE(stdout,'(/ / A)') "     vdW-DF kernel table calculated and written to file."
  !    OPEN(UNIT=21, FILE='kernel_table', STATUS='replace', FORM='formatted', ACTION='write')
  !    WRITE(21, '(2i5,f13.8)') Nqs, Nr_points
  !    WRITE(21, '(1p4e23.14)') r_max
  !    WRITE(21, '(1p4e23.14)') q_mesh
  !    DO q1_i = 1, Nqs
  !       DO q2_i = 1, q1_i
  !          WRITE(21, '(1p4e23.14)') kernel(:, q1_i, q2_i)
  !       END DO
  !    END DO
  !    DO q1_i = 1, Nqs
  !       DO q2_i = 1, q1_i
  !          WRITE(21, '(1p4e23.14)') d2phi_dk2(:, q1_i, q2_i)
  !       END DO
  !    END DO
  !    CLOSE (21)
  ! END IF
  !
  !
  ! Reading the kernel from an old kernel file.
  !
  ! IF (ionode) WRITE(stdout,'(/ / A)') "     vdW-DF kernel read from file."
  ! OPEN(UNIT=21, FILE='vdW_kernel_table', STATUS='old', FORM='formatted', ACTION='read')
  ! read(21, '(/ / / / / /)')
  ! DO q1_i = 1, Nqs
  !    DO q2_i = 1, q1_i
  !       READ(21, '(1p4e23.14)') kernel(:, q1_i, q2_i)
  !       kernel(:, q2_i, q1_i) = kernel(:, q1_i, q2_i)
  !    END DO
  ! END DO
  ! DO q1_i = 1, Nqs
  !    DO q2_i = 1, q1_i
  !       READ(21, '(1p4e23.14)')    d2phi_dk2(:, q1_i, q2_i)
  !       d2phi_dk2(:, q2_i, q1_i) = d2phi_dk2(:, q1_i, q2_i)
  !    END DO
  ! END DO
  ! CLOSE (21)


  DEALLOCATE( indices, proc_indices, phi, phi_deriv )


  ! --------------------------------------------------------------------
  ! Stop the timer.

  CALL stop_clock ( 'vdW_kernel' )

  END SUBROUTINE generate_kernel








  ! ####################################################################
  !                    |                            |
  !                    |  PREP_GAUSSIAN_QUADRATURE  |
  !                    |____________________________|
  !

  SUBROUTINE prep_gaussian_quadrature( weights )

  !! Routine to calculate the points and weights for the
  !! Gaussian-Legendre integration. This routine is modeled after the
  !! routine GAULEG from NUMERICAL RECIPES.
  
  REAL(DP), INTENT(INOUT) :: weights(:)
  ! The points and weights for the Gaussian-Legendre integration.

  INTEGER  :: Npoints
  ! The number of points we actually have to calculate. The rest will
  ! be obtained from symmetry.

  REAL(DP) :: poly_1, poly_2, poly_3
  ! Temporary storage for Legendre polynomials.

  INTEGER  :: i_point, i_poly
  ! Indexing variables.

  REAL(DP) :: root, dp_dx, last_root
  ! The value of the root of a given Legendre polynomial, the derivative
  ! of the polynomial at that root and the value of the root in the last
  ! iteration (to check for convergence of Newton's method).

  real(dp) :: midpoint, length
  ! The middle of the x-range and the length to that point.




  Npoints  = (Nintegration_points + 1)/2
  midpoint = 0.5D0 * ( ATAN(a_min) + ATAN(a_max) )
  length   = 0.5D0 * ( ATAN(a_max) - ATAN(a_min) )

  DO i_point = 1, Npoints
     ! -----------------------------------------------------------------
     ! Make an initial guess for the root.

     root = COS(DBLE(pi*(i_point - 0.25D0)/(Nintegration_points + 0.5D0)))

     DO
        ! --------------------------------------------------------------
        ! Use the recurrence relations to find the desired polynomial,
        ! evaluated at the approximate root. See NUMERICAL_RECIPES.

        poly_1 = 1.0D0
        poly_2 = 0.0D0

        DO i_poly = 1, Nintegration_points

           poly_3 = poly_2
           poly_2 = poly_1
           poly_1 = ((2.0D0 * i_poly - 1.0D0)*root*poly_2 - (i_poly-1.0D0)*poly_3)/i_poly

        END DO


        ! --------------------------------------------------------------
        ! Use the recurrence relations to find the desired polynomial.
        ! Find the derivative of the polynomial and use it in Newton's
        ! method to refine our guess for the root.

        dp_dx = Nintegration_points * (root*poly_1 - poly_2)/(root**2 - 1.0D0)

        last_root = root
        root      = last_root - poly_1/dp_dx


        ! --------------------------------------------------------------
        ! Check for convergence.

        IF (abs(root - last_root) <= 1.0D-14) EXIT

     END DO


     ! -----------------------------------------------------------------
     ! Fill in the array of evaluation points.

     a_points(i_point) = midpoint - length*root
     a_points(Nintegration_points + 1 - i_point) = midpoint + length*root


     ! -----------------------------------------------------------------
     ! Fill in the array of weights.

     weights(i_point) = 2.0D0 * length/((1.0D0 - root**2)*dp_dx**2)
     weights(Nintegration_points + 1 - i_point) = weights(i_point)

  END DO

  END SUBROUTINE prep_gaussian_quadrature








  ! ####################################################################
  !                            |             |
  !                            |  PHI_VALUE  |
  !                            |_____________|
  !
  
  REAL(DP) FUNCTION phi_value(d1, d2)

  !! This function returns the value of the kernel calculated via DION
  !! equation 14.
  
  REAL(DP), INTENT(IN) :: d1, d2
  ! The point at which to evaluate the kernel. d1 = q1*r and d2 = q2*r.

  REAL(DP) :: w, x, y, z, T
  ! Intermediate values.

  REAL(DP) :: nu(Nintegration_points), nu1(Nintegration_points)
  ! Defined in the discussio below equation 16 of DION.

  INTEGER  :: a_i, b_i
  ! Indexing variables.




  ! --------------------------------------------------------------------
  ! Loop over all integration points and calculate the value of the nu
  ! functions defined in the discussion below equation 16 in DION.

  DO a_i = 1, Nintegration_points
     nu(a_i)  = a_points2(a_i)/( 2.0D0 * h_function( a_points(a_i)/d1 ))
     nu1(a_i) = a_points2(a_i)/( 2.0D0 * h_function( a_points(a_i)/d2 ))
  END DO


  ! --------------------------------------------------------------------
  ! Carry out the integration of DION equation 13.

  phi_value = 0.0D0

  DO a_i = 1, Nintegration_points
     w = nu(a_i)
     y = nu1(a_i)
     DO b_i = 1, Nintegration_points
        x = nu(b_i)
        z = nu1(b_i)
        T = (1.0D0/(w+x) + 1.0D0/(y+z))*(1.0D0/((w+y)*(x+z)) + 1.0D0/((w+z)*(y+x)))
        phi_value = phi_value + T * W_ab(a_i, b_i)
     END DO
  END DO

  phi_value = 1.0D0/pi**2*phi_value

  END FUNCTION phi_value








  ! ####################################################################
  !                            |              |
  !                            |  RADIAL_FFT  |
  !                            |______________|
  !
  
  SUBROUTINE radial_fft(phi)
  
  !! This subroutine performs a radial Fourier transform on the
  !! real-space kernel functions.
  ! Basically, this is just:
  !            int(4*pi*r^2*phi*sin(k*r)/(k*r))dr
  ! integrated from 0 to r_max.
  ! That is, it is the kernel function phi integrated with the 0^th spherical
  ! Bessel function radially, with a 4*pi assumed from angular
  ! integration since we have spherical symmetry. The spherical symmetry
  ! comes in because the kernel function depends only on the magnitude
  ! of the vector between two points. The integration is done using the
  ! trapezoid rule.
  
  REAL(DP), INTENT(INOUT) :: phi(0:Nr_points)
  ! On input holds the real-space function phi_q1_q2(r).
  ! On output hold the reciprocal-space function phi_q1_q2(k).

  REAL(DP) :: phi_k(0:Nr_points)
  ! Temporary storage for phi_q1_q2(k).

  INTEGER  :: k_i, r_i
  ! Indexing variables.

  REAL(DP) :: r, k
  ! The real and reciprocal space points.




  phi_k = 0.0D0

  ! --------------------------------------------------------------------
  ! Handle the k=0 point separately.

  DO r_i = 1, Nr_points
     r        = r_i * dr
     phi_k(0) = phi_k(0) + phi(r_i)*r**2
  END DO


  ! --------------------------------------------------------------------
  ! Subtract half of the last value off because of the trapezoid rule.

  phi_k(0) = phi_k(0) - 0.5D0 * (Nr_points*dr)**2 * phi(Nr_points)


  ! --------------------------------------------------------------------
  ! Integration for the rest of the k-points.

  DO k_i = 1, Nr_points
     k = k_i * dk
     DO r_i = 1, Nr_points
        r          = r_i * dr
        phi_k(k_i) = phi_k(k_i) + phi(r_i) * r * SIN(k*r) / k
     END DO
     phi_k(k_i) = phi_k(k_i) - 0.5D0 * phi(Nr_points) * r * SIN(k*r) / k
  END DO


  ! --------------------------------------------------------------------
  ! Add in the 4*pi and the dr factor for the integration.

  phi = 4.0D0 * pi * phi_k * dr

  END SUBROUTINE radial_fft








  ! ####################################################################
  !                          |                  |
  !                          |  SET UP SPLINES  |
  !                          |__________________|
  !
  
  SUBROUTINE set_up_splines(phi, D2)
  
  !! This subroutine accepts a function (phi) and finds at each point the
  !! second derivative (D2) for use with spline interpolation. This
  !! function assumes we are using the expansion described in SOLER
  !! equation 3.
  ! That is, the derivatives are those needed to interpolate
  ! Kronecker delta functions at each of the q values. Other than some
  ! special modification to speed up the algorithm in our particular
  ! case, this algorithm is taken directly from NUMERICAL_RECIPES.

  REAL(DP), INTENT(IN)    :: phi(0:Nr_points)
  ! The k-space kernel function for a particular q1 and q2.

  REAL(DP), INTENT(INOUT) :: D2(0:Nr_points)
  ! The second derivatives to be used in the interpolation expansion
  ! (SOLER equation 3).

  REAL(DP), ALLOCATABLE   :: temp_array(:)         ! Temporary storage.
  REAL(DP)                :: temp_1, temp_2

  INTEGER  :: r_i
  ! Indexing variable.




  ALLOCATE( temp_array(0:Nr_points) )

  D2         = 0
  temp_array = 0

  DO r_i = 1, Nr_points - 1
     temp_1  = DBLE(r_i - (r_i - 1))/DBLE( (r_i + 1) - (r_i - 1) )
     temp_2  = temp_1 * D2(r_i-1) + 2.0D0
     D2(r_i) = (temp_1 - 1.0D0)/temp_2
     temp_array(r_i) = ( phi(r_i+1) - phi(r_i))/DBLE( dk*((r_i+1) - r_i) ) - &
          ( phi(r_i) - phi(r_i-1))/DBLE( dk*(r_i - (r_i-1)) )
     temp_array(r_i) = (6.0D0*temp_array(r_i)/DBLE( dk*((r_i+1) - (r_i-1)) )-&
          temp_1*temp_array(r_i-1))/temp_2
  END DO

  D2(Nr_points) = 0.0D0
  DO  r_i = Nr_points-1, 0, -1
     D2(r_i) = D2(r_i)*D2(r_i+1) + temp_array(r_i)
  END DO

  DEALLOCATE( temp_array )

  END SUBROUTINE set_up_splines








  ! ####################################################################
  !                          |            |
  !                          |  VDW_INFO  |
  !                          |____________|

  SUBROUTINE vdW_info

  IMPLICIT NONE




  WRITE(stdout,'(/)')
  WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"% You are using vdW-DF, which was implemented by the Thonhauser group. %")')
  WRITE(stdout,'(5x,"% Please cite the following two papers that made this development      %")')
  WRITE(stdout,'(5x,"% possible and the two reviews that describe the various versions:     %")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"%   T. Thonhauser et al., PRL 115, 136402 (2015).                      %")')
  WRITE(stdout,'(5x,"%   T. Thonhauser et al., PRB 76, 125112 (2007).                       %")')
  WRITE(stdout,'(5x,"%   K. Berland et al., Rep. Prog. Phys. 78, 066501 (2015).             %")')
  WRITE(stdout,'(5x,"%   D.C. Langreth et al., J. Phys.: Condens. Matter 21, 084203 (2009). %")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"% If you are calculating the stress with vdW-DF, please also cite:     %")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"%   R. Sabatini et al., J. Phys.: Condens. Matter 24, 424209 (2012).   %")')
  WRITE(stdout,'(5x,"%                                                                      %")')
  WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")')
  WRITE(stdout,'(/)')

  IF ( iverbosity > 0 ) THEN
     WRITE(stdout,'(5x,"Carrying out vdW-DF run using the following parameters:")')
     WRITE(stdout,'(5X,A,I3,A,I5,A,F8.3)' ) "Nqs    = ", Nqs, "  Npoints = ", Nr_points, &
                  "  r_max = ", r_max
     WRITE(stdout,'(5X,"q_mesh =",4F12.8)') (q_mesh(idx), idx=1, 4)
     WRITE(stdout,'(13X,4F12.8)') (q_mesh(idx), idx=5, Nqs)
  END IF

  END SUBROUTINE




END MODULE vdW_DF
