LCOV - code coverage report
Current view: top level - vgen - vdW.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 493 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 19 0.0 %

          Line data    Source code
       1             : #define POTENTIAL
       2             :       MODULE param
       3             :       USE m_constants
       4             :       IMPLICIT NONE
       5             : !
       6             : 
       7             :       REAL,PARAMETER   :: au2A =0.529177249
       8             :       REAL,PARAMETER   :: Ha2eV=hartree_to_ev_const
       9             :       REAL,PARAMETER   :: pi=3.14159265358979323846
      10             : !
      11             :       INTEGER, PARAMETER   :: jnlout=6
      12             :       INTEGER, PARAMETER   :: ftable=20,fdebug=21,fxcrysd=22,fwarn=23
      13             : !
      14             : ! parameters used to calculate the vdW-DF functional
      15             : !
      16             :       REAL,PARAMETER   ::  Zab_v1=-0.8491        ! vdW-DF1
      17             :       REAL,PARAMETER   ::  Zab_v2=-1.887         ! vdW-DF2
      18             : !
      19             : ! parameters used to calculate the LDA correlation energy
      20             : !
      21             :       REAL,PARAMETER   ::  A       = 0.031091    ! some fitted parameters for the
      22             :       REAL,PARAMETER   ::  alpha_1 = 0.2137      ! evaluation of LDA correlation
      23             :       REAL,PARAMETER   ::  beta_1  = 7.5957      ! energy according to:
      24             :       REAL,PARAMETER   ::  beta_2  = 3.5876      !
      25             :       REAL,PARAMETER   ::  beta_3  = 1.6382      ! J. P. Perdew and Yue Wang,
      26             :       REAL,PARAMETER   ::  beta_4  = 0.49294     ! Phys. Rev. B 45, 13244 (1992)
      27             :       REAL,PARAMETER   ::  beta_H  = 0.066725    ! for H function from PRL 77, 3865 (1996)
      28             :       REAL,PARAMETER   ::  gamma_H = 0.031091    ! for H function
      29             : !
      30             : ! parameters used to calculate the PBE exchange energy
      31             : !
      32             :       REAL,PARAMETER   ::  kappa_revPBE = 1.245   ! revPBE
      33             :       REAL,PARAMETER   ::  kappa_PBE    = 0.804   ! PBE from PRL 77, 3865 (1996)
      34             :       REAL,PARAMETER   ::  mu = 0.21951           ! PBE from PRL 77, 3865 (1996)
      35             : !
      36             :       END MODULE param
      37             : !
      38             : ! module containing the main non-local variables and subroutines
      39             : !
      40             :       MODULE nonlocal_data
      41             : 
      42             :       IMPLICIT NONE
      43             : !
      44             :       REAL             ::  Zab             ! This is the Zab really used. Zab_v1
      45             :                                                ! is the default. Can be switched to
      46             :                                                ! Zab_v2 by invoking the program with
      47             :                                                ! command line option 'vdw2'
      48             : !
      49             :       INTEGER              :: nx, ny, nz       ! number of grid points in x, y, z direction
      50             :       INTEGER              :: n_grid           ! total number of grid points
      51             :       INTEGER              :: n_k              ! number of k-points for which the kernel was
      52             :                                                ! tabulated
      53             : !
      54             :       REAL             :: r_max            ! maximum r for which the kernel has been generated
      55             : !
      56             :       REAL,ALLOCATABLE :: q_alpha(:)       ! q mesh for the interpolation
      57             :       REAL,ALLOCATABLE :: phi(:,:,:)       !
      58             :       REAL,ALLOCATABLE :: d2phi_dk2(:,:,:) !
      59             : !
      60             :       INTEGER                 :: n_gvectors    ! number of G vectors
      61             :       REAL                :: G_cut         ! radius of cutoff sphere
      62             :       REAL,   ALLOCATABLE :: G(:,:)        ! the G vectors
      63             :       INTEGER,ALLOCATABLE :: G_ind(:)      ! the index of the G vectors
      64             : !
      65             :       REAL             :: a1(3),a2(3),a3(3)  ! lattice vectors
      66             :       REAL             :: b1(3),b2(3),b3(3)  ! reciprocal lattice vectors
      67             : !
      68             :       INTEGER              :: n_alpha            ! number of q points
      69             :       REAL             :: q_cut              ! maximum q value
      70             :       INTEGER              :: m_c                ! maximum m for the saturation of q_0
      71             : !
      72             :       REAL             :: omega              ! volume of the unit cell
      73             :       REAL             :: tpibya             ! (2*pi/omega)
      74             : !
      75             :       REAL             :: lambda             ! parameter for the logarithmic q mesh
      76             :       REAL             :: dk                 ! spacing of the uniform radial k grid
      77             : !
      78             :       INTEGER              :: time1(8),time2(8)  ! arrays for measuring the execution time
      79             : !
      80             :       END MODULE nonlocal_data
      81             : !
      82             :       MODULE driver_fft
      83             : !
      84             :  CONTAINS
      85             : !
      86             : !==========================================================================================
      87             : !
      88             : !     this subroutine is an interface for the fft. It will fourier transform in place the
      89             : !     complex array fftin. idir = -1 means forward and idir = +1 means backward fourier
      90             : !     transform. At the moment it is designed to use fftw3.
      91             : 
      92           0 :       SUBROUTINE inplfft( fftin, idir )
      93             :       USE m_juDFT
      94             :       USE param,        ONLY: jnlout
      95             :       USE nonlocal_data,ONLY: nx,ny,nz,n_grid
      96             : 
      97             :       implicit none
      98             : 
      99             : 
     100             :       complex, intent(inout) :: fftin(n_grid)
     101             : 
     102             :       integer                    :: idir
     103             : 
     104             :       integer*8                  :: plan
     105             : #ifdef CPP_FFTW
     106             :       include 'fftw3.f'  ! some definitions for fftw3
     107             : 
     108             :       if ( idir == -1 ) then
     109             : 
     110             :          CALL dfftw_plan_dft_3d(plan,nz,ny,nx,fftin,fftin,FFTW_FORWARD,FFTW_ESTIMATE)
     111             :          CALL dfftw_execute_dft(plan,fftin,fftin)
     112             :          CALL dfftw_destroy_plan(plan)
     113             : 
     114             :          fftin = fftin / real(n_grid) ! rescale as fftw3 puts a factor n_grid on forward FFT
     115             : 
     116             :       elseif( idir == 1 ) then
     117             : 
     118             :          CALL dfftw_plan_dft_3d(plan,nz,ny,nx,fftin,fftin,FFTW_BACKWARD,FFTW_ESTIMATE)
     119             :          CALL dfftw_execute_dft(plan,fftin,fftin)
     120             :          CALL dfftw_destroy_plan(plan)
     121             : 
     122             :       else
     123             : 
     124             :         write(jnlout,*) 'ERROR during FFT: neither FORWARD &
     125             :                          nor BACKWARD FFT was chosen.'
     126             :         STOP 'Error in FFT'
     127             : 
     128             :       endif
     129             : #else
     130           0 :         CALL judft_error("VdW functionals only available if compiled with CPP_FFTW")
     131             : #endif
     132           0 :       END SUBROUTINE inplfft
     133             : !
     134             :       END MODULE driver_fft
     135             : !
     136             :       MODULE functionals
     137             : !
     138             :  CONTAINS
     139             : !
     140             : !================================================================================
     141             : !
     142           0 :       SUBROUTINE calc_PBE_correlation(n,grad_n,Ec_PBE,Ec_LDA,e_cLDA,e_cSL)
     143             : !
     144             :       USE param,        ONLY: Ha2eV,pi,        &
     145             :                               jnlout,             &
     146             :                               beta_H,gamma_H,     &                        ! Parameters for LDA_c
     147             :                               A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
     148             :       USE nonlocal_data,ONLY: n_grid,omega
     149             : !
     150             :       IMPLICIT NONE
     151             : !
     152             :       REAL,INTENT(in)  :: n(n_grid)          ! charge density
     153             :       REAL,INTENT(in)  :: grad_n(n_grid,3)   ! gradient of charge density
     154             :       REAL,INTENT(out) :: Ec_PBE             ! PBE-correlation energy
     155             :       REAL,INTENT(out) :: Ec_LDA             ! LDA-correlation energy
     156             :       REAL,INTENT(out) :: e_cLDA(n_grid)     ! LDA correlation energy density
     157             :       REAL,INTENT(out) :: e_cSL(n_grid)      ! SL correlation energy density in Ha
     158             : !
     159             :       REAL             :: r_s                ! intermediate values needed for the
     160             :       REAL             :: k_F                ! formulas given below.
     161             :       REAL             :: zeta               !
     162             :       REAL             :: phi_zeta           !
     163             :       REAL             :: grad_n_squ         !
     164             :       REAL             :: t                  !             -- "" --
     165             :       REAL             :: k_s                !
     166             :       REAL             :: AA                 !
     167             :       REAL             :: H                  !
     168             : !
     169             :       REAL             :: sqrt_r_s
     170             :       REAL             :: eLDA_c
     171             :       REAL             :: LDA_dummy_1
     172             :       REAL             :: LDA_dummy_2
     173             : !
     174             :       integer              :: i1
     175             : !
     176           0 :       e_cLDA(:)=0.0
     177           0 :       e_cSL(:) =0.0
     178             : !
     179           0 :       Ec_PBE=0.0
     180           0 :       Ec_LDA=0.0
     181             : !
     182             : !     Formula for PBE correlation ( Perdew, Burke and Ernzerhof PRL 77, 18 (1996) eqn. [3]):
     183             : !     Ec_PBE = int d^3 r n(r) (ec_LSDA (r_s, zeta) + H (r_s, zeta, t))
     184             : !
     185             : !     r_s = 3 / (4* pi*n)**(1/3) Seitz radius.
     186             : !     zeta is relative spin polarization will be set to 0 as vdW-DF only works without spin
     187             : !     t = | grad_n | / (2*phi(zeta)*k_s*n) is a dimensionless gradient
     188             : !     phi(zeta) = [ (1 + zeta)**(2/3) + (1 - zeta)**(2/3) ] / 2
     189             : !     k_s = sqrt( 4*k_F / (pi*a_0) )
     190             : !     a_0 = ( hbar / (m*e^2) ) = 1 in a. u.
     191             : !
     192             : !     H = (e**2/a_0) g * phi**3
     193             : !          * ln(1 + beta_H/gamma_H * t**2 * [1 + AA*t**2 / (1 + AA*t**2 + AA**2 * t**4 )] )
     194             : !     AA = beta_H/gamma_H * [ exp( -ec_LSDA / ( gamma_H * phi**3 * e**2 / a_0)) - 1 ]**(-1)
     195             : !
     196             : !
     197           0 :       zeta = 0.0     ! non-spin polarized case
     198             : !
     199             :       phi_zeta = ( (1.0 + zeta)**(2.0/3.0) + &
     200           0 :                    (1.0 - zeta)**(2.0/3.0))/2.0
     201             : !
     202           0 :       do i1=1,n_grid
     203             : !
     204           0 :          if(n(i1) < 1e-12) cycle
     205             : 
     206           0 :          r_s = (3.0/(4.0*pi*n(i1)))**(1.0/3.0)
     207             : 
     208             : !     we also need LDA correlation per particle so I repeat here the formulas from calc_q0
     209             : 
     210           0 :          sqrt_r_s = sqrt(r_s)
     211             : 
     212           0 :          LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
     213             :          LDA_dummy_2 =  2.0*A*(beta_1*sqrt_r_s     + beta_2*r_s + &
     214           0 :                                   beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
     215             : 
     216           0 :          eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
     217             : 
     218             : !     now start calculation of PBE correlation. first check wether density is very small (approx
     219             : !     zero) or negative which is also not correct
     220             : !
     221           0 :          k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
     222             : 
     223           0 :          k_s = sqrt( 4.0 * k_F / pi)
     224             : 
     225           0 :          grad_n_squ = grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2
     226             : 
     227           0 :          t = sqrt(grad_n_squ)/(2.0*phi_zeta*k_s*n(i1))
     228             : 
     229             :          AA = (beta_H/gamma_H)* &
     230           0 :               (1.0/(exp(-1.0*eLDA_c/(gamma_H*phi_zeta**3))- 1.0 ))
     231             : 
     232             :          H = gamma_H*phi_zeta**3 *         &
     233             :              log(1.0+                  &
     234             :                   (beta_H/gamma_H)*(t**2)* &
     235           0 :                   (1.0+AA*t**2)/(1.0+AA*t**2+AA**2 * t**4))
     236             : !
     237             : !     LDA and SL correlation energy density
     238             : !
     239           0 :          e_cLDA(i1) = eLDA_c*n(i1)
     240           0 :          e_cSL(i1)  =      H*n(i1)
     241             : !
     242             : #if 0
     243             :          Ec_PBE = Ec_PBE + n(i1)*eLDA_c + n(i1)*H
     244             :          Ec_LDA = Ec_LDA + n(i1)*eLDA_c
     245             : #else
     246           0 :          Ec_PBE = Ec_PBE + e_cLDA(i1) + e_cSL(i1)
     247           0 :          Ec_LDA = Ec_LDA + e_cLDA(i1)
     248             : #endif
     249             :       enddo
     250             : !
     251           0 :       Ec_PBE = Ec_PBE*omega/real(n_grid)
     252           0 :       Ec_LDA = Ec_LDA*omega/real(n_grid)
     253             : !
     254           0 :       END SUBROUTINE calc_PBE_correlation
     255             : !
     256             : !================================================================================
     257             : !
     258           0 :       SUBROUTINE calc_GGA_exchange(n,grad_n,Ex_PBE,Ex_revPBE,Ex_PW86,Ex_LDA)
     259             : !
     260             :       USE param,        ONLY: pi, &
     261             :                               mu,kappa_PBE,kappa_revPBE  ! Param. for PBE_ex
     262             :       USE nonlocal_data,ONLY: n_grid,omega
     263             : !
     264             :       IMPLICIT NONE
     265             : !
     266             :       REAL,INTENT(IN)  :: n(n_grid)          ! charge density
     267             :       REAL,INTENT(IN)  :: grad_n(n_grid,3)   ! gradient of charge density
     268             :       REAL,INTENT(OUT) :: Ex_PBE             !    PBE-exchange energy in Ha
     269             :       REAL,INTENT(OUT) :: Ex_revPBE          ! revPBE-exchange energy in Ha
     270             :       REAL,INTENT(OUT) :: Ex_PW86            ! refit PW86-exchange energy
     271             :                                                  ! in Ha
     272             :       REAL,INTENT(OUT) :: Ex_LDA             ! LDA part of the PBE exchange
     273             :                                                  ! energy Ex_PBE
     274             : !
     275             :       INTEGER              :: i1
     276             :       REAL             :: k_F                ! formulas given below.
     277             :       REAL             :: ss                 ! reduced density gradient
     278             :       REAL             :: sqrt_grad_n
     279             : !
     280             : ! ss is actually the reduced density gradient
     281             : !                      ss = |\grad n|/2(3\pi^2)^1/3n^4/3
     282             : !
     283           0 :       Ex_PBE   =0.0
     284           0 :       Ex_revPBE=0.0
     285           0 :       Ex_PW86  =0.0
     286           0 :       Ex_LDA   =0.0
     287             : !
     288             : ! here the charge density has already the correct indices
     289             : !
     290           0 :       do i1=1,n_grid
     291             : !
     292           0 :       if(n(i1) < 1e-12) cycle
     293             : !
     294           0 :         k_F = (3.0*pi*pi*n(i1))**(1.0/3.0)
     295             : !
     296           0 :         sqrt_grad_n = sqrt(grad_n(i1,1)**2 + grad_n(i1,2)**2 + grad_n(i1,3)**2)
     297             : !
     298           0 :         ss=sqrt_grad_n/(2.0*k_F*n(i1))
     299             : !
     300             : ! LDA part of the PBE exchange energy density
     301             : !
     302             :         Ex_LDA = Ex_LDA +                                            &
     303             :                  (-3.0/4.0)*(3.0/pi)**(1.0/3.0)*      &
     304           0 :                  (n(i1)**(4.0/3.0))
     305             : !
     306             : ! PBE exchange energy
     307             : !
     308             :         Ex_PBE = Ex_PBE +                                            &
     309             :                  (-3.0/4.0)*(3.0/pi)**(1.0/3.0)*      &
     310             :                  (n(i1)**(4.0/3.0))*                           &
     311             :                  (1.0+kappa_PBE-                                  &
     312           0 :                          kappa_PBE/(1.0+mu*ss**2.0/kappa_PBE))
     313             : !
     314             : ! revPBE exchange energy
     315             : !
     316             :         Ex_revPBE = Ex_revPBE +                                      &
     317             :                  (-3.0/4.0)*(3.0/pi)**(1.0/3.0)*      &
     318             :                  (n(i1)**(4.0/3.0))*                           &
     319             :                  (1.0+kappa_revPBE-                               &
     320           0 :                          kappa_revPBE/(1.0+mu*ss**2.0/kappa_revPBE))
     321             : !
     322             : ! use refit PW86
     323             : !
     324             :         Ex_PW86 = Ex_PW86 +                                          &
     325             :                  (-3.0/4.0)*(3.0/pi)**(1.0/3.0)*      &
     326             :                  (n(i1)**(4.0/3.0))*                           &
     327             :                  (1.0+15.0*0.1234*ss**2.0 +              &
     328           0 :                   17.33*ss**4.0+0.163*ss**6.0)**(1.0/15.0)
     329             :       enddo
     330             : !
     331           0 :       Ex_PBE    = Ex_PBE   *omega/real(n_grid)
     332           0 :       Ex_revPBE = Ex_revPBE*omega/real(n_grid)
     333           0 :       Ex_PW86   = Ex_PW86  *omega/real(n_grid)
     334           0 :       Ex_LDA    = Ex_LDA   *omega/real(n_grid)
     335             : !
     336           0 :       END SUBROUTINE calc_GGA_exchange
     337             : !
     338             : !==========================================================================================
     339             : !
     340           0 :       SUBROUTINE calc_ehartree(n)
     341             :       USE param,        ONLY: Ha2ev,pi,       &
     342             :                               jnlout
     343             :       USE nonlocal_data,ONLY: n_grid,n_gvectors, &
     344             :                               omega,             &
     345             :                               G,G_ind
     346             :       USE driver_fft
     347             : !
     348             :       IMPLICIT NONE
     349             : !
     350             :       REAL,INTENT(in)     :: n(n_grid)
     351             : !
     352             :       INTEGER                 :: i1,ind
     353             :       REAL                :: E_Hartree
     354             :       REAL                :: fac,size
     355           0 :       COMPLEX,ALLOCATABLE :: n_cmplx(:)
     356             : !
     357             : ! E_Hartree= 2*pi*omega \sum_{G/=0} n(G)^2/G^2
     358             : !
     359           0 :       allocate(n_cmplx(n_grid))
     360             : !
     361           0 :       n_cmplx = (0.0,0.0)
     362           0 :       n_cmplx(1:n_grid) = cmplx( n(1:n_grid), 0.0 )
     363             : !
     364           0 :       call inplfft( n_cmplx, -1 )
     365             : !
     366           0 :       E_Hartree=0.0
     367             : !
     368           0 :       do i1 = 1, n_gvectors
     369             : !
     370           0 :       ind = G_ind(i1)
     371           0 :       size=sum(G(i1,1:3)**2)
     372           0 :       if (size > 1e-12) then
     373           0 :          fac=1.0/size
     374             :       else
     375           0 :       print'(A,F10.3)','Total nr. of electrons: ',omega*real(n_cmplx(ind))
     376           0 :          fac=0.0
     377             :       endif
     378             : !
     379             :       E_Hartree=E_Hartree+ &
     380           0 :                 (real(n_cmplx(ind))**2+aimag(n_cmplx(ind))**2)*fac
     381             : !
     382             :       enddo
     383             : !
     384           0 :       deallocate(n_cmplx)
     385             : !
     386           0 :       E_Hartree = E_Hartree*omega*2.0*pi
     387             : !
     388           0 :       write(jnlout,'(/,A)') '-------------- Hartree energy  --------------'
     389           0 :       write(jnlout,'(A)')      'E_Hartree (Ha/Ry/eV):'
     390           0 :       write(jnlout,'(3F24.15)') E_Hartree,E_Hartree*2.0,E_Hartree*Ha2eV
     391             : !
     392           0 :       END SUBROUTINE calc_ehartree
     393             : !
     394             :       END MODULE functionals
     395             : !
     396             :       MODULE plot_functions
     397             : !
     398             :  CONTAINS
     399             : !
     400             : !==========================================================================================
     401             : !
     402           0 :       SUBROUTINE write_xcrysden_LDA(file_name,ene_dens)
     403             : !
     404             :       USE param,        ONLY: au2A,fxcrysd
     405             :       USE nonlocal_data,ONLY: nx,ny,nz,n_grid, &
     406             :                               a1,a2,a3
     407             : !
     408             :       IMPLICIT NONE
     409             : !
     410             :       CHARACTER(LEN=*),INTENT(IN) :: file_name
     411             :       REAL,        INTENT(IN) :: ene_dens(n_grid)
     412             :       INTEGER                     :: i1,ind,ix,iy,iz
     413           0 :       REAL,ALLOCATABLE        :: ene_temp(:,:,:)
     414             : !
     415           0 :       allocate(ene_temp(nx,ny,nz))
     416             : !
     417           0 :       open(fxcrysd,FILE=trim(file_name))
     418             : !
     419           0 :       write(fxcrysd,'(A)') 'CRYSTAL'
     420           0 :       write(fxcrysd,'(A)') 'PRIMVEC'
     421           0 :       write(fxcrysd,'(3F16.10)') a1(:)*au2A
     422           0 :       write(fxcrysd,'(3F16.10)') a2(:)*au2A
     423           0 :       write(fxcrysd,'(3F16.10)') a3(:)*au2A
     424             : !
     425           0 :       write(fxcrysd,'(A)')       'PRIMCOORD'
     426           0 :       write(fxcrysd,'(A,2x,I1)') 'XXX',1
     427           0 :       write(fxcrysd,'(A)')       'Atomic_positions'
     428             : !
     429           0 :       write(fxcrysd,'(A)') 'BEGIN_BLOCK_DATAGRID_3D'
     430           0 :       write(fxcrysd,'(A)') '3D_PWSCF'
     431           0 :       write(fxcrysd,'(A)') 'DATAGRID_3D_UNKNOWN'
     432           0 :       write(fxcrysd,'(3(3x,I5))') nx,ny,nz
     433           0 :       write(fxcrysd,'(3F16.10)') 0.0,0.0,0.0 ! origin of the system
     434           0 :       write(fxcrysd,'(3F16.10)') a1(:)*au2A
     435           0 :       write(fxcrysd,'(3F16.10)') a2(:)*au2A
     436           0 :       write(fxcrysd,'(3F16.10)') a3(:)*au2A
     437             : !
     438           0 :       do ix=1,nx
     439           0 :         do iy=1,ny
     440           0 :           do iz=1,nz
     441           0 :             ind = ix + nx*(iy - 1) + nx*ny*(iz - 1)
     442           0 :             ene_temp(ix,iy,iz)=ene_dens(ind)
     443             :             ! in the case of the semi-local corr. ene. dens.
     444             :             !   ene_temp(ix,iy,iz) can be positive
     445             :             !if(ene_temp(ix,iy,iz) > 0.0) then
     446             :             !    print'(A)','WARNING: ene_temp(ix,iy,iz) > 0.0!'
     447             :             !endif
     448             :           enddo
     449             :         enddo
     450             :       enddo
     451             : !
     452           0 :       do i1=1,nz
     453           0 :         write(fxcrysd,'(5(1x,E15.8))') ene_temp(1:nx,1:ny,i1)
     454             :       enddo
     455             : !
     456             : ! write the end of the XSF file
     457             : !
     458           0 :       write(fxcrysd,'(A)') 'END_DATAGRID_3D'
     459           0 :       write(fxcrysd,'(A)') 'END_BLOCK_DATAGRID_3D'
     460             : !
     461           0 :       close(fxcrysd)
     462             : !
     463           0 :       deallocate(ene_temp)
     464             : !
     465           0 :       END SUBROUTINE write_xcrysden_LDA
     466             : !
     467             : !==========================================================================================
     468             : !
     469           0 :       SUBROUTINE write_xcrysden_NL(file_name,ene_dens)
     470             : !
     471             :       USE param,        ONLY: au2A,fxcrysd
     472             :       USE nonlocal_data,ONLY: nx,ny,nz,       &
     473             :                               n_grid,n_alpha, &
     474             :                               a1,a2,a3
     475             : !
     476             :       IMPLICIT NONE
     477             : !
     478             :       CHARACTER(LEN=*),INTENT(IN) :: file_name
     479             :       REAL,        INTENT(IN) :: ene_dens(n_grid,n_alpha)
     480             :       INTEGER                     :: i1,ind,ix,iy,iz
     481           0 :       REAL,ALLOCATABLE        :: ene_temp(:,:,:)
     482             : !
     483           0 :       allocate(ene_temp(nx,ny,nz))
     484             : !
     485           0 :       open(fxcrysd,FILE=trim(file_name))
     486             : !
     487           0 :       write(fxcrysd,'(A)') 'CRYSTAL'
     488           0 :       write(fxcrysd,'(A)') 'PRIMVEC'
     489           0 :       write(fxcrysd,'(3F16.10)') a1(:)*au2A
     490           0 :       write(fxcrysd,'(3F16.10)') a2(:)*au2A
     491           0 :       write(fxcrysd,'(3F16.10)') a3(:)*au2A
     492             : !
     493           0 :       write(fxcrysd,'(A)')       'PRIMCOORD'
     494           0 :       write(fxcrysd,'(A,2x,I1)') 'XXX',1
     495           0 :       write(fxcrysd,'(A)')       'Atomic_positions'
     496             : !
     497           0 :       write(fxcrysd,'(A)') 'BEGIN_BLOCK_DATAGRID_3D'
     498           0 :       write(fxcrysd,'(A)') '3D_PWSCF'
     499           0 :       write(fxcrysd,'(A)') 'DATAGRID_3D_UNKNOWN'
     500           0 :       write(fxcrysd,'(3(3x,I5))') nx,ny,nz
     501           0 :       write(fxcrysd,'(3F16.10)') 0.0,0.0,0.0 ! origin of the system
     502           0 :       write(fxcrysd,'(3F16.10)') a1(:)*au2A
     503           0 :       write(fxcrysd,'(3F16.10)') a2(:)*au2A
     504           0 :       write(fxcrysd,'(3F16.10)') a3(:)*au2A
     505             : !
     506           0 :       do ix=1,nx
     507           0 :         do iy=1,ny
     508           0 :           do iz=1,nz
     509           0 :             ind = ix + nx*(iy - 1) + nx*ny*(iz - 1)
     510           0 :             ene_temp(ix,iy,iz)=sum(ene_dens(ind,1:n_alpha))
     511             :           enddo
     512             :         enddo
     513             :       enddo
     514             : !
     515           0 :       do i1=1,nz
     516           0 :         write(fxcrysd,'(5(1x,E15.8))') ene_temp(1:nx,1:ny,i1)
     517             :       enddo
     518             : !
     519             : ! write the end of the XSF file
     520             : !
     521           0 :       write(fxcrysd,'(A)') 'END_DATAGRID_3D'
     522           0 :       write(fxcrysd,'(A)') 'END_BLOCK_DATAGRID_3D'
     523             : !
     524           0 :       close(fxcrysd)
     525             : !
     526           0 :       deallocate(ene_temp)
     527             : !
     528           0 :       END SUBROUTINE write_xcrysden_NL
     529             : !
     530             :       END MODULE plot_functions
     531             : !
     532             :       MODULE nonlocal_funct
     533             :       PRIVATE
     534             :       PUBLIC :: soler
     535             : !
     536             :  CONTAINS
     537             : !
     538             : !========================================================================================
     539             : !
     540           0 :       SUBROUTINE soler(n, Ecnl, v_nl)
     541             :       USE param,        ONLY: Ha2eV,au2A,jnlout
     542             :       USE nonlocal_data,ONLY: n_grid,n_alpha, &
     543             :                               omega,          &
     544             :                               q_alpha,phi,    &
     545             :                               G,G_ind,        &
     546             :                               time1,time2
     547             :       USE functionals
     548             :       USE driver_fft
     549             :       USE plot_functions
     550             : !
     551             :       IMPLICIT NONE
     552             : !
     553             :       REAL,          INTENT(INOUT) :: n(n_grid)       ! charge density
     554             :       REAL,          INTENT(INOUT) :: v_nl(n_grid)
     555             :       REAL                         :: Ecnl            ! the non-local correlation
     556             : 
     557             :       CHARACTER(LEN=132)               :: option
     558             : !
     559           0 :       REAL,ALLOCATABLE    :: grad_n(:,:)     ! gradient of charge density
     560           0 :       REAL,ALLOCATABLE    :: e_cLDA(:)       ! LDA correlation energy density
     561           0 :       REAL,ALLOCATABLE    :: e_cSL(:)        ! semi local correlation energy density
     562           0 :       REAL,ALLOCATABLE    :: e_cNL(:,:)      ! non local correlation energy density
     563           0 :       REAL,ALLOCATABLE    :: q_0(:)          ! array holding the q0 on the grid
     564             : 
     565             : #ifdef POTENTIAL
     566             : 
     567           0 :       REAL,ALLOCATABLE    :: dq0_dn(:)
     568           0 :       REAL,ALLOCATABLE    :: dq0_dgrad_n(:)
     569             : 
     570             : #endif
     571             : 
     572           0 :       COMPLEX,ALLOCATABLE :: theta_alpha(:,:)! array holding theta_i
     573           0 :       COMPLEX,ALLOCATABLE :: u_a(:,:)        ! quantity needed for the non local
     574             :                                                  ! correlation energy density
     575             :       REAL                :: Ecnl_rsp        ! Ecnl calculated in real space
     576             :       REAL                :: Ec_LDA_rsp      ! Ec_LDA calculated in real space
     577             :       REAL                :: Ec_LDA          ! LDA correlation
     578             :       REAL                :: Ec_PBE          ! PBE correlation energy. will be replaced by
     579             :                                                  ! LDA correlation + non local correlation
     580             :       REAL                :: Ex_PBE          ! PBE exchange energy
     581             :       REAL                :: Ex_revPBE       ! revPBE echange energy
     582             :       REAL                :: Ex_PW86         ! refit PW86 exchange energy
     583             :       REAL                :: Ex_LDA          ! LDA part of the PBE exchange
     584             :                                                  ! energy Ex_PBE
     585             :       INTEGER                 :: i, j, k
     586             :       INTEGER                 :: ind ,mem_size
     587             :       REAL                :: sum_test
     588             : !
     589             : !     the actual program starts here. allocation for n and grad_n have to be reconsidered when
     590             : !     putting this as a subroutine into a different program.
     591             : !
     592             : !     read the pregenerated kernel from file 'vdW_kernel_table'. This will set variables n_alpha
     593             : !     n_k and dk. The arrays phi(n_alpha,n_alpha,n_k) and d2phi_dk2(n_alpha,n_alpha,n_k) needed for the
     594             : !     later interpolation of the kernel will be allocated and read.
     595             : !
     596           0 :       call read_kernel()
     597             : !
     598           0 :       allocate(theta_alpha(n_grid,n_alpha))
     599           0 :       allocate(u_a(n_grid,n_alpha))
     600           0 :       allocate(grad_n(n_grid,3),q_0(n_grid))
     601           0 :       allocate(e_cLDA(n_grid),e_cSL(n_grid),e_cNL(n_grid,n_alpha))
     602             : !
     603             :       mem_size=2*n_grid*n_alpha + 2*n_grid*n_alpha + n_grid*3 + n_grid + &
     604           0 :                n_grid + n_grid + n_grid*n_alpha
     605             :       write(jnlout,'(A,F12.3,A,/)')  &
     606           0 :                    'Memory required: ',real(mem_size)*8.0/(1024.0*1024.0),' MB'
     607             : !
     608             : !     setup the G vectors which we need for the calculation of the gradient, the non local correlation
     609             : !     energy and the non local contribution to the potential.
     610             : !
     611           0 :       call DATE_AND_TIME(values=time1)
     612             : !
     613           0 :       call setup_g_vectors()
     614             : !
     615           0 :       call DATE_AND_TIME(values=time2)
     616             : !
     617           0 :       write(jnlout,'(A)') 'time to setup G vectors:'
     618           0 :       call timing( time2 - time1 )
     619             : !
     620           0 :       write(jnlout,*)
     621             : !
     622             : !     calculate the gradient of the density numerically by FFT.
     623             : !
     624           0 :       call DATE_AND_TIME(values=time1)
     625             : !
     626           0 :       call calc_gradient( n, grad_n )
     627             : !
     628           0 :       call DATE_AND_TIME(values=time2)
     629             : !
     630           0 :       write(jnlout,'(A)') 'time to calculate gradient of n numerically:'
     631           0 :       call timing( time2 - time1 )
     632             : !
     633             : !     calculate q_0 for every grid point according to equations (11),(12) from Dion and (7)
     634             : !     from Soler. The latter cares for the saturation.
     635             : !
     636           0 :       call DATE_AND_TIME(values=time1)
     637             : !
     638           0 :       call calc_q0(n, grad_n, q_0)
     639             : !
     640           0 :       call DATE_AND_TIME(values=time2)
     641             : !
     642           0 :       write(jnlout,'(A)') 'time to calculate q0:'
     643           0 :       call timing( time2 - time1 )
     644             : !
     645             : !     calculate theta_i defined as theta_i = n(r_i) * p_alpha(q_i) where p_alpha are the
     646             : !     polynomials obtained by interpolating dirac delta. These polynomials will also be
     647             : !     derived in this subroutine
     648             : !
     649           0 :       call DATE_AND_TIME(values=time1)
     650             : !
     651           0 :       call calc_theta_i(n, q_alpha, q_0, theta_alpha)
     652             : !
     653           0 :       call DATE_AND_TIME(values=time2)
     654             : !
     655           0 :       write(jnlout,'(A)') 'time to setup Thetas:'
     656           0 :       call timing( time2 - time1 )
     657             : !
     658             : !     to calculate the non local energy density we need theta_alpha_i in real space. As they are
     659             : !     overwritten during the fourier transform we save them here.
     660             : !
     661           0 :       e_cNL = theta_alpha
     662             : !
     663             : !     Fourier transform the theta_alpha_i in order to get theta_alpha_k.
     664             : !
     665           0 :       call DATE_AND_TIME(values=time1)
     666             : !
     667           0 :       do i=1,n_alpha
     668             : 
     669           0 :          call inplfft( theta_alpha(:,i), -1 )
     670             : 
     671             :       enddo
     672             : !
     673           0 :       call DATE_AND_TIME(values=time2)
     674             : !
     675           0 :       write(jnlout,'(A)') 'time spent to Fourier transform Thetas:'
     676           0 :       call timing( time2 - time1 )
     677             : !
     678             : !     now we can evaluate the integral in eqn. (8) of Soler.
     679             : !
     680           0 :       call calc_ecnl(Ecnl, theta_alpha, u_a, e_cNL)
     681             : !
     682             : #ifdef POTENTIAL
     683             : 
     684           0 :       allocate(dq0_dn(n_grid))
     685           0 :       allocate(dq0_dgrad_n(n_grid))
     686             : 
     687           0 :       call calc_dq0(n, grad_n, dq0_dn, dq0_dgrad_n)
     688             : 
     689           0 :       call calc_potential(v_nl, u_a, q_0, dq0_dn, dq0_dgrad_n, n, grad_n)
     690             : 
     691           0 :       deallocate( dq0_dn, dq0_dgrad_n)
     692             : 
     693             : #endif
     694             : 
     695             : !     In vdW-DF GGA-correlation is been replaced by LDA-correlation + non local correlation. Thus
     696             : !     we have to calculate Ec_LDA and Ec_PBE now to output th energy term which has to be added
     697             : !     to the total energy. We already have the LDA correlation energy density from calc_q0 so we
     698             : !     just have to integrate that array and afterwards call a subroutine calculating the PBE.
     699             : !     omega/n_grid is the according volume element.
     700             : !
     701           0 :       call calc_PBE_correlation(n,grad_n,Ec_PBE,Ec_LDA,e_cLDA,e_cSL)
     702             : !
     703           0 :       write(jnlout,'(/,A)')     'PBE_c (Ha/Ry/eV):'
     704           0 :       write(jnlout,'(3F24.15)') Ec_PBE,Ec_PBE*2.0,Ec_PBE*Ha2eV
     705             : !
     706           0 :       write(jnlout,'(/,A)')     'LDA_c (Ha/Ry/eV):'
     707           0 :       write(jnlout,'(3F24.15)') Ec_LDA,Ec_LDA*2.0,Ec_LDA*Ha2eV
     708             : !
     709             : !     for testing write the sum over LDA energy density.
     710             : !
     711           0 :       Ec_LDA_rsp = sum(e_cLDA)*omega/real(n_grid)
     712             :       write(jnlout,'(/,A)')     &
     713           0 :                    '(check) LDA_c evaluated in real space (Ha/Ry):'
     714           0 :       write(jnlout,'(2F24.15)') Ec_LDA_rsp,Ec_LDA_rsp*2.0
     715             : !
     716           0 :       write(jnlout,'(/,A)')     'E_c,nl (Ha/Ry/eV):'
     717           0 :       write(jnlout,'(3F24.15)') Ecnl,Ecnl*2.0,Ecnl*Ha2eV
     718             : !
     719             : !     for testing write the sum over non local energy density.
     720             : !     it should be equal to Ec_nl
     721             : !
     722           0 :       Ecnl_rsp=sum(e_cNL)*omega/real(n_grid)
     723             :       write(jnlout,'(/,A)')     &
     724           0 :                    '(check) E_c,nl evaluated in real space (Ha/Ry):'
     725           0 :       write(jnlout,'(2F24.15)') 0.5*Ecnl_rsp,Ecnl_rsp
     726             : !
     727           0 :       call write_xcrysden_LDA('Ec_LDA.xsf',e_cLDA)
     728           0 :       call write_xcrysden_LDA('Ec_SL.xsf' ,e_cSL)
     729           0 :       call write_xcrysden_NL ('Ec_NL.xsf' ,e_cNL)
     730             : !
     731           0 :       call calc_GGA_exchange(n,grad_n,Ex_PBE,Ex_revPBE,Ex_PW86,Ex_LDA)
     732             : !
     733           0 :       write(jnlout,'(/,A)')     'Ex_LDA (Ha/Ry/eV):'
     734           0 :       write(jnlout,'(3F24.15)') Ex_LDA,Ex_LDA*2.0,Ex_LDA*Ha2eV
     735             : !
     736           0 :       write(jnlout,'(/,A)') '-------------- PBE exchange --------------'
     737           0 :       write(jnlout,'(A)')      'Ex_PBE (Ha/Ry/eV):'
     738           0 :       write(jnlout,'(3F24.15)') Ex_PBE,Ex_PBE*2.0,Ex_PBE*Ha2eV
     739             : !
     740           0 :       write(jnlout,'(A)') 'Ex_PBE+Ec_LDA+Ecnl (Ha/Ry/eV):'
     741           0 :       sum_test=Ex_PBE+Ec_LDA+Ecnl
     742           0 :       write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
     743             : !
     744           0 :       write(jnlout,'(/,A)') '-------------- revPBE exchange --------------'
     745           0 :       write(jnlout,'(A)')      'Ex_revPBE (Ha/Ry/eV):'
     746           0 :       write(jnlout,'(3F24.15)') Ex_revPBE,Ex_revPBE*2.0,Ex_revPBE*Ha2eV
     747             : !
     748           0 :       write(jnlout,'(A)') 'Ex_revPBE+Ec_LDA+Ecnl (Ha/Ry/eV):'
     749           0 :       sum_test=Ex_revPBE+Ec_LDA+Ecnl
     750           0 :       write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
     751             : !
     752           0 :       if (trim(option) == 'vdw2' .or. trim(option) == 'vdW2') then
     753             : !
     754             : ! in the case of vdW-DF2 use refit PW86 exchange
     755             : !
     756           0 :         write(jnlout,'(/,A)') '-------------- for vdW-DF2 --------------'
     757           0 :         write(jnlout,'(A)') 'refit PW86 exchange:'
     758           0 :         write(jnlout,'(A)')      'Ex_PW86 (Ha/Ry/eV):'
     759           0 :         write(jnlout,'(3F24.15)') Ex_PW86,Ex_PW86*2.0,Ex_PW86*Ha2eV
     760             : !
     761           0 :         write(jnlout,'(A)') 'Ex_PW86+Ec_LDA+Ecnl (Ha/Ry/eV):'
     762           0 :         sum_test=Ex_PW86+Ec_LDA+Ecnl
     763           0 :         write(jnlout,'(3F24.15)') sum_test,sum_test*2.0,sum_test*Ha2eV
     764             : !
     765             :       endif
     766             : !
     767           0 :       call calc_ehartree(n)
     768             : !
     769           0 :       deallocate(grad_n,G,G_ind)
     770           0 :       deallocate(e_cNL,e_cSL,e_cLDA)
     771           0 :       deallocate(q_0,theta_alpha,u_a)
     772             : !
     773           0 :       END SUBROUTINE soler
     774             : !
     775             : !==========================================================================================
     776             : !
     777             : !     We take the kernel from the Thonhauser implementation (Espresso). This subroutine to
     778             : !     read the kernel thus is also taken from there.
     779             : 
     780           0 :       SUBROUTINE read_kernel()
     781             :       USE param,        ONLY: pi,ftable,jnlout
     782             :       USE nonlocal_data,ONLY: n_alpha,n_k,          &
     783             :                               r_max,q_cut,dk,       &
     784             :                               q_alpha,phi,d2phi_dk2
     785             :       implicit none
     786             : 
     787             :       character(len=30)    :: double_format = '(1p4e23.14)'
     788             :       integer              :: q1, q2
     789           0 :       if (allocated(q_alpha)) return
     790           0 :       write(jnlout,*) 'reading kernel table'
     791             : 
     792           0 :       open(ftable, file='vdW_kernel_table')
     793             : 
     794           0 :       read(ftable, '(2i5)' ) n_alpha, n_k
     795           0 :       read(ftable, double_format) r_max
     796             : 
     797           0 :       allocate(q_alpha(n_alpha), phi(0:n_k,n_alpha,n_alpha), d2phi_dk2(0:n_k,n_alpha,n_alpha))
     798             : 
     799           0 :       read(ftable, double_format) q_alpha
     800             : 
     801           0 :       q_cut = q_alpha(n_alpha)
     802             : 
     803           0 :       write(jnlout,*)'n_a:', n_alpha
     804           0 :       write(jnlout,*)'q_c:', q_cut
     805           0 :       write(jnlout,*)
     806             : 
     807           0 :       do q1 = 1, n_alpha
     808           0 :          do q2 = 1, q1
     809             : 
     810           0 :             read(ftable, double_format ) phi(0:n_k, q1, q2)
     811           0 :             phi(0:n_k,q2, q1) = phi(0:n_k, q1, q2)
     812             : 
     813             :          end do
     814             :       end do
     815             : 
     816           0 :       do q1 = 1, n_alpha
     817           0 :          do q2 = 1, q1
     818             : 
     819           0 :             read(ftable, double_format ) d2phi_dk2(0:n_k,q1, q2)
     820           0 :             d2phi_dk2(0:n_k,q2, q1) = d2phi_dk2(0:n_k,q1, q2)
     821             : 
     822             :          end do
     823             :       end do
     824             : 
     825           0 :       dk = 2.0*pi/r_max
     826             : !
     827           0 :       close(ftable)
     828             : !
     829             :       END SUBROUTINE read_kernel
     830             : !
     831             : !==========================================================================================
     832             : !
     833           0 :       SUBROUTINE setup_g_vectors()
     834             :       USE param,        ONLY: jnlout
     835             :       USE nonlocal_data,ONLY: nx,ny,nz,n_gvectors,n_grid, &
     836             :                               G_cut,b1,b2,b3,             &
     837             :                               G,G_ind
     838             :       implicit none
     839             : !
     840             :       integer             :: igx,igy,igz
     841             :       integer             :: tmpx,tmpy,tmpz
     842             :       integer             :: ind
     843             : !
     844             :       real            :: G_squ
     845             :       real            :: gx, gy, gz
     846             : 
     847             : !     We need the G vectors in a number of subroutines of this code. So the idea is to set them up
     848             : !     once and for all to minimize possible sources of errors. The array G(n_gvectors,3) will hold
     849             : !     the components of the g_vectors and G_ind their index on the fft mesh. We have to do the loop
     850             : !     over all G vectors twice. First time to count number of G_vectors inside cut off sphere and
     851             : !     second time to set them up.
     852             : 
     853           0 :       n_gvectors = 0
     854             : 
     855           0 :       do igz = -(nz-1)/2,(nz-1)/2
     856           0 :          do igy = -(ny-1)/2,(ny-1)/2
     857           0 :             do igx = -(nx-1)/2,(nx-1)/2
     858             : 
     859           0 :                gx = igx * b1(1) + igy * b2(1) + igz * b3(1)
     860           0 :                gy = igx * b1(2) + igy * b2(2) + igz * b3(2)
     861           0 :                gz = igx * b1(3) + igy * b2(3) + igz * b3(3)
     862             : 
     863           0 :                G_squ =  gx**2 + gy**2 + gz**2
     864             : 
     865           0 :                if ( G_squ .le. G_cut ) n_gvectors = n_gvectors + 1
     866             : 
     867             :             enddo
     868             :          enddo
     869             :       enddo
     870           0 :       if (allocated(g)) return
     871           0 :       allocate(G(n_gvectors, 3), G_ind(n_gvectors))
     872             : 
     873           0 :       write(jnlout,*) "#g-vectors:",n_gvectors," outof ",nz*ny*nx
     874             : 
     875           0 :       G     = 0.0
     876           0 :       G_ind = 0
     877             : 
     878           0 :       n_gvectors = 0
     879             : 
     880           0 :       do igz = -(nz-1)/2,(nz-1)/2
     881           0 :          do igy = -(ny-1)/2,(ny-1)/2
     882           0 :             do igx = -(nx-1)/2,(nx-1)/2
     883             : 
     884           0 :                gx = igx * b1(1) + igy * b2(1) + igz * b3(1)
     885           0 :                gy = igx * b1(2) + igy * b2(2) + igz * b3(2)
     886           0 :                gz = igx * b1(3) + igy * b2(3) + igz * b3(3)
     887             : 
     888           0 :                G_squ =  gx**2 + gy**2 + gz**2
     889             : 
     890           0 :                tmpx = 0
     891           0 :                tmpy = 0
     892           0 :                tmpz = 0
     893             : 
     894           0 :                IF (igx .LT. 0) tmpx = nx
     895           0 :                IF (igy .LT. 0) tmpy = ny
     896           0 :                IF (igz .LT. 0) tmpz = nz
     897             : 
     898             :                ind = 1 + ( igx + tmpx ) + &
     899             :                       nx*( igy + tmpy ) + &
     900           0 :                    nx*ny*( igz + tmpz )
     901             : 
     902           0 :                if ( G_squ .le. G_cut ) then
     903             : 
     904           0 :                    n_gvectors = n_gvectors + 1
     905             : 
     906           0 :                    G(n_gvectors, 1) = gx
     907           0 :                    G(n_gvectors, 2) = gy
     908           0 :                    G(n_gvectors, 3) = gz
     909             : 
     910           0 :                    G_ind(n_gvectors) = ind
     911             : 
     912             :                endif
     913             : 
     914             :             enddo
     915             :          enddo
     916             :       enddo
     917             : 
     918             :       end subroutine setup_g_vectors
     919             : !
     920             : !==========================================================================================
     921             : !
     922           0 :       SUBROUTINE calc_gradient(n,grad_n)
     923             :       USE param,        ONLY: jnlout
     924             :       USE nonlocal_data,ONLY: n_grid,n_gvectors,nx,ny,nz, &
     925             :                               a1,a2,a3,G,G_ind
     926             :       USE driver_fft
     927             : !
     928             :       implicit none
     929             : !
     930             :       real, intent(in)     :: n(n_grid)
     931             :       real, intent(inout)  :: grad_n(n_grid,3)
     932             : 
     933           0 :       complex, allocatable :: n_cmplx(:)
     934           0 :       complex, allocatable :: grad_n_cmplx(:,:)
     935             : 
     936             :       integer                  :: i
     937             : 
     938             :       integer                  :: ind
     939             : 
     940             :       integer*8                :: plan
     941             : 
     942           0 :       grad_n = 0.0
     943             : 
     944             : !     A derivative d/dx f(x) in real space corresponds to i*G*f(G) in reciprocal space. So we have
     945             : !     to fourier transform the density, find the according G vectors, calculate the product and back
     946             : !     fourier transform the gradient
     947             : 
     948           0 :       allocate(n_cmplx(n_grid), grad_n_cmplx(n_grid,3))
     949             : 
     950           0 :       n_cmplx      = (0.0, 0.0)
     951           0 :       grad_n_cmplx = (0.0, 0.0)
     952             : 
     953           0 :       do i = 1, n_grid
     954             : 
     955           0 :          n_cmplx(i) = cmplx( n(i), 0.0)
     956             : 
     957             :       enddo
     958             : 
     959           0 :       call inplfft( n_cmplx, -1 )
     960             : 
     961           0 :       do i = 1, n_gvectors
     962             : 
     963           0 :          ind = G_ind(i)
     964             : 
     965           0 :          grad_n_cmplx(ind,:) = (0.0,1.0) * G(i,:) * n_cmplx(ind)
     966             : 
     967             :       enddo
     968             : 
     969           0 :       do i=1,3
     970             : 
     971           0 :         call inplfft( grad_n_cmplx(:,i), 1)
     972             : 
     973             :       enddo
     974             : !
     975           0 :       grad_n(:,:) = real(grad_n_cmplx(:,:))
     976             : 
     977           0 :       deallocate(n_cmplx, grad_n_cmplx)
     978             : !
     979           0 :       END SUBROUTINE calc_gradient
     980             : !
     981             : !==========================================================================================
     982             : !
     983             : !     This subroutine calculates q_0 following eqns. 11 and 12 of Dion and saturates it
     984             : !     following eqn. 7 of Soler
     985             : !
     986           0 :       SUBROUTINE calc_q0(n, grad_n, q_0)
     987             :       USE param,        ONLY: pi, &
     988             :                               A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
     989             :       USE nonlocal_data,ONLY: n_grid,q_cut,m_c,Zab
     990             : !
     991             :       implicit none
     992             : !
     993             :       real, intent(in)    ::  n(n_grid)        ! charge density
     994             :       real, intent(in)    ::  grad_n(n_grid,3) ! gradient of charge density
     995             :       real, intent(out)   ::  q_0(n_grid)      ! array holding g_0 on the grid
     996             : 
     997             :       real                ::  q                ! q before saturation
     998             : 
     999             :       real                ::  k_F                ! fermi wave vector
    1000             :       real                ::  r_s                !
    1001             :       real                ::  sqrt_r_s           !
    1002             :       real                ::  eLDA_c             ! LDA correlation
    1003             :       real                ::  eLDA_x             ! LDA exchange
    1004             :       real                ::  grad_n_squ         ! squared gradient of n
    1005             : 
    1006             :       real                ::  LDA_dummy_1
    1007             :       real                ::  LDA_dummy_2
    1008             : 
    1009             :       real                ::  sum_m
    1010             : 
    1011             :       integer                 ::  i, m         ! counters
    1012             : !
    1013             : !     loop over all grid points
    1014             : !
    1015           0 :       q_0(:) = q_cut
    1016             : !
    1017           0 :       do i=1,n_grid
    1018             : 
    1019             : !     check if charge density is negative. If so treat it like zero charge density. Zero
    1020             : !     density corresponds to high q_0. Thats why we set it to q_cut the highest possible
    1021             : !     q_0 and continue with the next grid point.
    1022             : 
    1023           0 :          if ( n(i) < 1e-12 ) then
    1024           0 :             q_0(i) = q_cut
    1025           0 :             cycle
    1026             :          end if
    1027             : !
    1028             : !     calculate q according to eqns (11) and (12) from Dion.
    1029             : !
    1030           0 :          k_F = (3.0*(pi**2.0)*n(i))**(1.0/3.0)
    1031           0 :          r_s = (3.0/(4.0*pi*n(i)))**(1.0/3.0)
    1032           0 :          sqrt_r_s = r_s**(1.0/2.0)
    1033             : 
    1034           0 :          grad_n_squ = grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0
    1035             : 
    1036           0 :          LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
    1037             :          LDA_dummy_2 =  2.0*A*(beta_1*sqrt_r_s     + beta_2*r_s + &
    1038           0 :                                   beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
    1039             : 
    1040           0 :          eLDA_x = (-3.0/(4.0*pi))*k_F
    1041           0 :          eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
    1042             : !
    1043             : !     LDA correlation energy density is epsilon_c[n(r)]*n(r)
    1044             : !
    1045             : !         e_cLDA(i) = eLDA_c*n(i)
    1046             : !
    1047             :          q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ &
    1048           0 :                      /(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F
    1049             : 
    1050             : !     now saturate q according to eq. (7) from Soler.
    1051             : 
    1052           0 :          sum_m = 0.0
    1053             : 
    1054           0 :          do m=1,m_c
    1055             : 
    1056           0 :             sum_m = sum_m + (q/q_cut)**m / real(m)
    1057             : 
    1058             :          enddo
    1059             : 
    1060           0 :          q_0(i) = q_cut*(1.0 - exp(-sum_m))
    1061             : 
    1062             :       enddo
    1063             : 
    1064           0 :       END SUBROUTINE calc_q0
    1065             : !
    1066             : !==========================================================================================
    1067             : !
    1068             : !
    1069           0 :       SUBROUTINE calc_theta_i(n, q_alpha, q_0, theta_alpha)
    1070             : 
    1071             :       USE nonlocal_data,ONLY: n_grid, n_alpha
    1072             : !
    1073             :       implicit none
    1074             : !
    1075             :       real, intent(in)      ::  n(n_grid)        ! charge density
    1076             :       real, intent(in)      ::  q_alpha(n_alpha) ! q-mesh for interpolation
    1077             :       real, intent(in)      ::  q_0(n_grid)      ! q values on the grid
    1078             : 
    1079             :       complex, intent(inout) ::  theta_alpha(n_grid,n_alpha) ! thetas on real space grid
    1080             : 
    1081             :       integer                    ::  i               ! counter
    1082             : 
    1083             : !     this call will interpolate dirac delta which gives p_alpha(q_i) according to a method
    1084             : !     from numerical recipes. Whithin this subroutine also the initial setup of the second
    1085             : !     derivatives needed for spline interpolation will be done once.
    1086             : 
    1087           0 :       call splint( q_alpha, q_0, theta_alpha )
    1088             : 
    1089           0 :       do i=1,n_grid
    1090             : 
    1091             : !     theta_alpha will hold at this stage the p_alpha(q_i)
    1092             : 
    1093           0 :             theta_alpha(i,:) = n(i)*theta_alpha(i,:)
    1094             : 
    1095             :       enddo
    1096             : 
    1097           0 :       END SUBROUTINE calc_theta_i
    1098             : !
    1099             : !==========================================================================================
    1100             : !
    1101           0 :       SUBROUTINE calc_ecnl(Ecnl, theta_alpha, u_a, e_cNL)
    1102             :       USE param        ,ONLY: pi,jnlout
    1103             :       USE nonlocal_data,ONLY: nx,ny,nz,n_grid,n_alpha,n_gvectors, &
    1104             :                               omega,time1,time2,                  &
    1105             :                               G, G_ind
    1106             :       USE driver_fft
    1107             : !
    1108             :       implicit none
    1109             : !
    1110             :       real,   intent(out)   :: Ecnl
    1111             :       complex,intent(in)    :: theta_alpha(n_grid,n_alpha)
    1112             :       complex,intent(out)   :: u_a(n_grid,n_alpha)   ! function u_a(r) needed for the construction
    1113             :                                                          ! of the potential and also for non local
    1114             :       real,   intent(inout) :: e_cNL(n_grid,n_alpha) ! non local correlation energy density                                                      ! correlation energy density.
    1115             : !
    1116             :       integer                 :: i, ind
    1117             :       integer                 :: alpha,beta
    1118             :       real                :: k
    1119           0 :       real                :: phi_k(n_alpha,n_alpha)
    1120             :       complex             :: integral
    1121             : !
    1122           0 :       u_a = (0.0, 0.0)
    1123           0 :       Ecnl = 0.0
    1124             : !
    1125             : !     integration of theta_a*theta_b*phi_ab
    1126             : !
    1127           0 :       write (jnlout,*) 'calculating E_c,nl:'
    1128             : !
    1129           0 :       call DATE_AND_TIME(values=time1)
    1130             : !
    1131             : !     The difference to older versions of this code is the way how the G-vectors are set up. The
    1132             : !     index is not any more calculated on the fly but stored in G_ind.
    1133             : !
    1134           0 :       do i = 1, n_gvectors
    1135             : 
    1136           0 :           ind = G_ind(i)
    1137             : 
    1138           0 :           k = sqrt(dble( G(i,1)**2 + G(i,2)**2 + G(i,3)**2))
    1139             : 
    1140           0 :           call interpolate_kernel(k, phi_k)
    1141             : 
    1142           0 :           integral = (0.0, 0.0)
    1143             : 
    1144           0 :           do alpha = 1,n_alpha
    1145           0 :              do beta = 1,n_alpha
    1146             : 
    1147             :                 integral = integral + conjg(theta_alpha(ind,alpha)) * &
    1148           0 :                     theta_alpha(ind,beta)*cmplx(phi_k(alpha,beta),0.0) !&
    1149             : !
    1150             : !    the array u_a(k,a) = sum_b theta_b(k) * phi_ab(k) equals the fourier transform
    1151             : !    of the function u_a(r) = sum_b int d^3 r' theta_b(r') phi_ab(r - r') which we need for
    1152             : !    the calculation of the potential and the non local correlation energy density.
    1153             : !
    1154             :                 u_a(ind,alpha) = u_a(ind,alpha) + &
    1155           0 :                     theta_alpha(ind,beta)*cmplx(phi_k(alpha,beta),0.0)
    1156             : 
    1157             :              enddo
    1158             :           enddo
    1159             : 
    1160           0 :           Ecnl = Ecnl + real(integral)
    1161             : 
    1162             :       enddo
    1163             : !
    1164             : !     backward fourier transform of u_a(k,a) in order to get u_a(r) the non local correlation
    1165             : !     energy density.
    1166             : !
    1167           0 :       do i=1,n_alpha
    1168             : 
    1169           0 :          call inplfft( u_a(:,i), 1 )
    1170             : 
    1171             :       enddo
    1172             : !
    1173             : !     at this moment e_cNL(i,alpha) holds the theta_alpha_i which we have to multiply by the convolution
    1174             : !     u_a(r) in order to get the non local energy density as a function of alpha
    1175             : !
    1176           0 :       e_cNL = e_cNL*real(u_a)
    1177             : !
    1178           0 :       call DATE_AND_TIME(values=time2)
    1179             : !
    1180           0 :       write(jnlout,'(A)') 'time for integration:'
    1181           0 :       call timing( time2 - time1 )
    1182             : !
    1183             : !     taking care of the units. (2*pi)^3/omega is the volume element of the reciprocal unit
    1184             : !     cell. omega^2 we get when replacing the integrals by sums over a real space grid. The
    1185             : !     (2*pi)^3 cancels with an according factor from the radial fourier transform on the kernel
    1186             : !
    1187           0 :       Ecnl = Ecnl*0.5*omega
    1188             : !
    1189           0 :       END SUBROUTINE calc_ecnl
    1190             : !
    1191             : !==========================================================================================
    1192             : !
    1193             :       SUBROUTINE make_q_alpha(q_alpha)
    1194             : 
    1195             : 
    1196             :       USE nonlocal_data,ONLY: n_alpha,q_cut,lambda ! n_alpha: number of q points, q_cut: maximum
    1197             :                                                    ! q value, lambda: parameter for logarithmic
    1198             :                                                    ! mesh.
    1199             :       implicit none
    1200             : !
    1201             :       real, intent(out) ::  q_alpha(n_alpha) ! array holding the qs
    1202             : 
    1203             :       real              ::  q1               ! auxiliary vairable holding the first
    1204             :                                                  ! q which can be calculated explicitely
    1205             : 
    1206             :       integer               ::  i                ! counter
    1207             : 
    1208             : !     first we have to setup the qs starting from the maximum value q_cut on a logarithmic mesh
    1209             : !     which satisfies ( q_a+1 - q_a ) = lambda*( q_a - q_a-1 ). Lambda > 1 is a parameter.
    1210             : !     In GPAW lambda = 1.2 is used.
    1211             : 
    1212             :       q1 = q_cut * ( lambda - 1.0) / (lambda**( n_alpha - 1.0 ) - 1.0 )
    1213             : 
    1214             :       do i = 1,n_alpha
    1215             : 
    1216             :          q_alpha(i) = q1 * (lambda**( i - 1.0 ) - 1.0 ) / ( lambda - 1.0 )
    1217             : 
    1218             :       enddo
    1219             : 
    1220             :       END SUBROUTINE make_q_alpha
    1221             : !
    1222             : !==========================================================================================
    1223             : !
    1224             : !     From Numerical Recipes
    1225             : 
    1226           0 :       SUBROUTINE splint( x_i, x, p_iofx  )
    1227             : 
    1228             : 
    1229             :       USE nonlocal_data,ONLY: n_grid, n_alpha
    1230             : !
    1231             :       implicit none
    1232             : !
    1233             :       real, intent(in)  ::  x_i(n_alpha) ! the x_i values where the function is known
    1234             :       real, intent(in)  ::  x(n_grid)    ! the x values for which the function
    1235             :                                              ! shall be interpolated
    1236             : 
    1237             :       complex, intent(inout) ::  p_iofx(n_grid,n_alpha) ! the function values for each x
    1238             : 
    1239           0 :       real, allocatable :: d2y_dx2(:,:) ! second derivatives needed for the splines
    1240             : 
    1241           0 :       real, allocatable :: y(:)
    1242             : 
    1243             :       real              :: a, b, c, d, h ! some intermediate variables for the interpolation
    1244             : 
    1245             :       integer               :: i, j               ! counters
    1246             :       integer               :: u_bound, l_bound   ! variables for finding the section of x_i
    1247             :                                                   ! in which x is located.
    1248             :       integer               :: alpha              ! index of the found section
    1249             : 
    1250           0 :       allocate(y(n_alpha))
    1251             : 
    1252           0 :       allocate(d2y_dx2(n_alpha,n_alpha))
    1253             : 
    1254           0 :       call setup_spline(x_i,d2y_dx2)
    1255             : 
    1256           0 :       do i=1,n_grid
    1257             : 
    1258             : !     first find the correct section of x_i in which x is located by bisectioning. According to
    1259             : !     numerical recipes this is efficient if the x values are random. In our case there might
    1260             : !     be some correlation as the density and its gradient are smooth.
    1261             : 
    1262           0 :          l_bound = 1
    1263           0 :          u_bound = n_alpha
    1264             : 
    1265           0 :          do while ((u_bound - l_bound) > 1)
    1266             : 
    1267           0 :             alpha = (u_bound + l_bound) / 2
    1268             : 
    1269           0 :             if ( x(i) > x_i(alpha)) then
    1270             :                l_bound = alpha
    1271             :             else
    1272             :                u_bound = alpha
    1273             :             endif
    1274             :          enddo
    1275             : 
    1276           0 :          h = x_i(u_bound) - x_i(l_bound)
    1277             : 
    1278           0 :          a = (x_i(u_bound) - x(i))/h
    1279           0 :          b = (x(i) - x_i(l_bound))/h
    1280           0 :          c = ((a**3.0 - a)*h**2.0)/6.0
    1281           0 :          d = ((b**3.0 - b)*h**2.0)/6.0
    1282             : 
    1283           0 :          do j=1,n_alpha
    1284             : 
    1285           0 :             y    = 0.0
    1286           0 :             y(j) = 1.0
    1287             : 
    1288             :             p_iofx(i,j) = a*y(l_bound) + b*y(u_bound) + &
    1289           0 :                c*d2y_dx2(j,l_bound) + d*d2y_dx2(j,u_bound)
    1290             : 
    1291             :          enddo
    1292             :       enddo
    1293             : 
    1294           0 :       deallocate(y, d2y_dx2)
    1295             : 
    1296           0 :       END SUBROUTINE splint
    1297             : !
    1298             : !==========================================================================================
    1299             : !
    1300             : !     From Numerical Recipes
    1301             : 
    1302           0 :       SUBROUTINE setup_spline(x_i,d2y_dx2)
    1303             : 
    1304             : 
    1305             :       USE nonlocal_data,ONLY: n_alpha
    1306             : !
    1307             :       implicit none
    1308             : !
    1309             :       real, intent(in)     ::  x_i(n_alpha)
    1310             :       real, intent(inout)  ::  d2y_dx2(n_alpha,n_alpha)
    1311             : 
    1312           0 :       real, allocatable    ::  y(:)       ! this array holds the function values at x_i which
    1313             :                                               ! are going to be interpolated.
    1314           0 :       real, allocatable    ::  dy_dx(:)   ! temporary array holding the first derivative
    1315             : 
    1316             :       real                 ::  sig, p     ! temporary values for the interpolation
    1317             : 
    1318             :       integer                  ::  i, j       ! counter
    1319             : 
    1320           0 :       allocate(y(n_alpha), dy_dx(n_alpha))
    1321             : 
    1322           0 :       do i=1,n_alpha
    1323             : 
    1324             : !     as we are interpolating dirac delta set the according function values:
    1325             : 
    1326           0 :          y=0.0
    1327           0 :          y(i)=1.0
    1328             : 
    1329             : !     now according to numerical recipes we set the boundary conditions which will give
    1330             : !     so called "natural" splines.
    1331             : 
    1332           0 :          d2y_dx2(i,1) = 0.0
    1333           0 :          dy_dx(1) = 0.0
    1334             : 
    1335           0 :          do j=2,n_alpha-1
    1336             : 
    1337           0 :             sig = (x_i(j) - x_i(j-1)) / (x_i(j+1) - x_i(j-1))
    1338           0 :             p = sig*d2y_dx2(i,j-1) + 2.0
    1339           0 :             d2y_dx2(i,j) = (sig - 1.0)/p
    1340           0 :             dy_dx(j) = (y(j+1) - y(j))/(x_i(j+1) - x_i(j)) - (y(j) - y(j-1))/(x_i(j) - x_i(j-1))
    1341           0 :             dy_dx(j) = (6.0*dy_dx(j)/(x_i(j+1) - x_i(j-1)) - sig*dy_dx(j-1))/p
    1342             : 
    1343             :          enddo
    1344             : 
    1345           0 :          d2y_dx2(i,n_alpha) = 0.0
    1346             : 
    1347           0 :          do j=n_alpha - 1, 1, -1
    1348             : 
    1349           0 :             d2y_dx2(i,j) = d2y_dx2(i,j)*d2y_dx2(i,j+1) + dy_dx(j)
    1350             : 
    1351             :          enddo
    1352             :       enddo
    1353             : 
    1354           0 :       deallocate(y, dy_dx)
    1355             : 
    1356           0 :       END SUBROUTINE setup_spline
    1357             : !
    1358             : !===========================================================================================
    1359             : !
    1360             : !     similar to the Thonhauser implementation
    1361             : 
    1362           0 :       SUBROUTINE interpolate_kernel(k, phi_k)
    1363             : 
    1364             : 
    1365             :       USE nonlocal_data,ONLY: n_alpha, dk, n_k, phi, d2phi_dk2
    1366             : !
    1367             :       implicit none
    1368             : !
    1369             :       real, intent(in)    :: k
    1370             :       real, intent(inout) :: phi_k(n_alpha,n_alpha)
    1371             : 
    1372             :       real                :: a, b, c, d
    1373             : 
    1374             :       integer                 :: q1, q2, k_i
    1375             : 
    1376             : !     the algorithm for interpolation will be more or less the same as in splint().
    1377           0 :       phi_k = 0.0
    1378             : 
    1379             : !     find the interval in which our k lies
    1380             : 
    1381           0 :       k_i = int(k/dk)
    1382             : 
    1383             : !     simple case when k equals one of the values we have tabulated the kernel for
    1384             : 
    1385           0 :       if (mod(k,dk) == 0.0) then
    1386             : 
    1387           0 :          do q1=1,n_alpha
    1388           0 :             do q2=1,q1
    1389             : 
    1390           0 :                phi_k(q1, q2) = phi(k_i, q1, q2)
    1391           0 :                phi_k(q2, q1) = phi(k_i, q2, q1)
    1392             : 
    1393             :             enddo
    1394             :          enddo
    1395             : 
    1396             :          return
    1397             : 
    1398             :       endif
    1399             : 
    1400           0 :       a = (dk*(k_i+1.0) - k)/dk
    1401           0 :       b = (k - dk*k_i)/dk
    1402           0 :       c = (a**3.0-a)*dk**2.0/6.0
    1403           0 :       d = (b**3.0-b)*dk**2.0/6.0
    1404             : 
    1405           0 :       do q1 = 1, n_alpha
    1406           0 :          do q2 = 1, q1
    1407             : 
    1408             :             phi_k(q1, q2) = a*phi(k_i, q1, q2) + b*phi(k_i+1, q1, q2) &
    1409           0 :             +(c*d2phi_dk2(k_i, q1, q2) + d*d2phi_dk2(k_i+1,q1, q2))
    1410             : 
    1411           0 :             phi_k(q2, q1) = phi_k(q1, q2)
    1412             : 
    1413             :          end do
    1414             :       end do
    1415             : 
    1416             :       END SUBROUTINE interpolate_kernel
    1417             : !
    1418             : !=================================================================================================
    1419             : !
    1420             : #ifdef POTENTIAL
    1421             : 
    1422           0 :       subroutine calc_dq0(n, grad_n, dq0_dn, dq0_dgrad_n)
    1423             : 
    1424             :       USE param,        ONLY: pi, &
    1425             :                               A,alpha_1,beta_1,beta_2,beta_2,beta_3,beta_4 ! Parameters for LDA_c
    1426             :       USE nonlocal_data,ONLY: n_grid,q_cut,m_c,Zab
    1427             : 
    1428             : 
    1429             :       implicit none
    1430             : 
    1431             :       real, intent(in)    ::  n(n_grid)        ! charge density
    1432             :       real, intent(in)    ::  grad_n(n_grid,3) ! gradient of charge density
    1433             : 
    1434             :       real, intent(out)   ::  dq0_dn(n_grid)   ! derivative of q0 w.r.t. the density
    1435             :       real, intent(out)   ::  dq0_dgrad_n(n_grid) ! derivative of q0 w.r.t. the gradient
    1436             : 
    1437             :       real                ::  q                ! q before saturation
    1438             : 
    1439             :       real                ::  dq0_dq           ! derivative needed for dq0_dn and dq0_dgrad_n
    1440             : 
    1441             :       real                ::  k_F                ! fermi wave vector
    1442             :       real                ::  r_s                !
    1443             :       real                ::  sqrt_r_s           !
    1444             :       real                ::  eLDA_c             ! LDA correlation
    1445             :       real                ::  eLDA_x             ! LDA exchange
    1446             :       real                ::  grad_n_squ         ! squared gradient of n
    1447             : 
    1448             :       real                ::  LDA_dummy_1
    1449             :       real                ::  LDA_dummy_2
    1450             : 
    1451             :       real                ::  sum_m
    1452             : 
    1453             :       integer                 ::  i, m         ! counters
    1454             : 
    1455             : !     loop over all grid points
    1456             : 
    1457           0 :       dq0_dn(:) = 0.0
    1458           0 :       dq0_dgrad_n(:) = 0.0
    1459             : 
    1460           0 :       do i=1,n_grid
    1461             : 
    1462           0 :          if ( n(i) < 1e-12 ) then
    1463             :             cycle     ! NOTE: The derivatives will be zero at these points
    1464             :          end if
    1465             : !
    1466             : !     calculate q according to eqns (11) and (12) from Dion.
    1467             : !
    1468           0 :          k_F = (3.0*(pi**2.0)*n(i))**(1.0/3.0)
    1469           0 :          r_s = (3.0/(4.0*pi*n(i)))**(1.0/3.0)
    1470           0 :          sqrt_r_s = r_s**(1.0/2.0)
    1471             : 
    1472           0 :          grad_n_squ = grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0
    1473             : 
    1474           0 :          LDA_dummy_1 = (8.0*pi/3.0)*A*(1.0+ alpha_1*r_s)
    1475             :          LDA_dummy_2 =  2.0*A*(beta_1*sqrt_r_s     + beta_2*r_s + &
    1476           0 :                                   beta_3*sqrt_r_s*r_s + beta_4*r_s*r_s)
    1477             : 
    1478           0 :          eLDA_x = (-3.0/(4.0*pi))*k_F
    1479           0 :          eLDA_c = (-3.0/(4.0*pi))*LDA_dummy_1*log(1.0+1.0/LDA_dummy_2)
    1480             : 
    1481             :          q = (1.0 + ( eLDA_c / eLDA_x) - (Zab/9.0)*grad_n_squ &
    1482           0 :                      /(4.0*(n(i)**2.0)*(k_F**2.0)))*k_F
    1483             : 
    1484             : !     now saturate q according to eq. (7) from Soler.
    1485             : 
    1486           0 :          sum_m = 0.0
    1487           0 :          dq0_dq = 0.0
    1488             : 
    1489           0 :          do m=1,m_c
    1490             : 
    1491           0 :             sum_m = sum_m + (q/q_cut)**m / real(m)
    1492             : 
    1493           0 :             dq0_dq = dq0_dq + ((q/q_cut)**(m-1))
    1494             : 
    1495             :          enddo
    1496             : 
    1497           0 :          dq0_dq = dq0_dq * exp(-sum_m)
    1498             : 
    1499             : !     here we calculate the derivatives needed for the potential later. i.e. we calculate
    1500             : !     n*dq0/dn and n*dq0/dgrad_n so that we do not have to divide by n which might be very
    1501             : !     small at some points and thus produce numerical errors.
    1502             : 
    1503             :          dq0_dn(i) = dq0_dq * ( k_F/3.0 &
    1504             :                      + k_F*7.0/3.0*(Zab/9.0)*grad_n_squ &
    1505             :                      /(4.0*(n(i)**2.0)*(k_F**2.0)) &
    1506             :                      - (8.0*pi/9.0)*A*alpha_1*r_s*log(1.0+1.0/LDA_dummy_2) &
    1507             :                      + LDA_dummy_1/(LDA_dummy_2*(1.0 + LDA_dummy_2)) &
    1508             :                      * (2.0*A*(beta_1/6.0*sqrt_r_s + beta_2/3.0*r_s &
    1509           0 :                      + beta_3/2.0*r_s*sqrt_r_s + 2.0*beta_4/3.0*r_s**2)))
    1510             : 
    1511             :          dq0_dgrad_n(i) = dq0_dq * 2.0 * (-1.0*Zab/9.0) &
    1512           0 :                           * sqrt(grad_n_squ) / (4.0*n(i)*k_F)
    1513             : 
    1514             :       enddo
    1515             : 
    1516           0 :       end subroutine calc_dq0
    1517             : !
    1518             : !=======================================================================================
    1519             : !
    1520             : !     similar to the Thonhauser implementation
    1521             : 
    1522           0 :       subroutine calc_potential(v_nl, u_a, q_0, dq0_dn, dq0_dgrad_n, n, grad_n)
    1523             : 
    1524             : 
    1525             :       use nonlocal_data, only : n_grid, n_alpha, q_alpha, q_cut, nx, ny, nz, G, G_ind, n_gvectors
    1526             : 
    1527             :       USE driver_fft
    1528             : 
    1529             :       implicit none
    1530             : 
    1531             :       real, intent(out)    :: v_nl(n_grid) ! non local part of the potential
    1532             : 
    1533             :       real, intent(in)     :: q_0(n_grid)
    1534             : 
    1535             :       real, intent(in)     :: dq0_dn(n_grid)      ! n*dq0/dn
    1536             :       real, intent(in)     :: dq0_dgrad_n(n_grid) ! n*dq0/dgrad_n
    1537             : 
    1538             :       real, intent(in)     :: n(n_grid)
    1539             :       real, intent(in)     :: grad_n(n_grid,3)
    1540             : 
    1541             :       complex, intent(in)  :: u_a(n_grid, n_alpha)
    1542             : 
    1543           0 :       complex, allocatable :: h_j(:,:)
    1544             : 
    1545           0 :       real, allocatable    :: d2y_dx2(:,:) ! second derivatives needed for the splines
    1546             : 
    1547           0 :       real, allocatable    :: y(:)
    1548             : 
    1549             :       real                 :: a, b, c, d, d1, d2, h ! some intermediate variables
    1550             :                                                         ! for the interpolation of p_a
    1551             : 
    1552             :       real                 :: p_a
    1553             :       real                 :: dpa_dq0
    1554             : 
    1555             :       real                 :: grad_n_abs ! | grad_n |
    1556             : 
    1557             :       integer                  :: i, alpha
    1558             : 
    1559             :       integer                  :: l_bound, u_bound
    1560             : 
    1561             :       integer                  :: ind
    1562             : 
    1563           0 :       allocate(h_j(n_grid,3))
    1564           0 :       allocate(y(n_alpha))
    1565             : 
    1566           0 :       h_j = (0.0, 0.0)
    1567             : 
    1568           0 :       allocate( d2y_dx2(n_alpha, n_alpha) )
    1569             : 
    1570           0 :       call setup_spline( q_alpha, d2y_dx2 )
    1571             : 
    1572             : !     The potential will be calculated using FFT following White and Bird. PRB 50, 4954
    1573             : !     (1994) eqn. (10). v^xc(r) = d f_xc / d n(r) - (1/N) * sum_G,r' i * G * ( grad_n(r')
    1574             : !     / |grad_n(r')| ) * ( d f_xc / d |grad_n(r')| ) * e^(i*G*(r - r')) with E_c,NL =
    1575             : !     int d^3 r f_xc .
    1576             : 
    1577           0 :       do i=1,n_grid
    1578             : 
    1579             : !     first we have to interpolate the polynomials p_a and their derivatives with respect
    1580             : !     to q0 again as we need them for the calculation of the potential. This will be the
    1581             : !     same procedure as in splint.
    1582             : 
    1583           0 :          l_bound = 1
    1584           0 :          u_bound = n_alpha
    1585             : 
    1586           0 :          do while ((u_bound - l_bound) > 1)
    1587             : 
    1588           0 :             ind = (u_bound + l_bound) / 2
    1589             : 
    1590           0 :             if ( q_0(i) > q_alpha(ind)) then
    1591             :                l_bound = ind
    1592             :             else
    1593             :                u_bound = ind
    1594             :             endif
    1595             :          enddo
    1596             : 
    1597           0 :          h = q_alpha(u_bound) - q_alpha(l_bound)
    1598             : 
    1599           0 :          a = (q_alpha(u_bound) - q_0(i))/h
    1600           0 :          b = (q_0(i) - q_alpha(l_bound))/h
    1601           0 :          c = ((a**3.0 - a)*h**2.0)/6.0
    1602           0 :          d = ((b**3.0 - b)*h**2.0)/6.0
    1603           0 :          d1 = (3.0*a**2.0 - 1.0)*h/6.0
    1604           0 :          d2 = (3.0*b**2.0 - 1.0)*h/6.0
    1605             : 
    1606           0 :          do alpha = 1, n_alpha
    1607             : 
    1608           0 :             y = 0.0
    1609           0 :             y(alpha) = 1.0
    1610             : 
    1611             :             p_a = a*y(l_bound) + b*y(u_bound) + &
    1612           0 :                   c*d2y_dx2(alpha,l_bound) + d*d2y_dx2(alpha,u_bound)
    1613             : 
    1614             :             dpa_dq0 = (y(u_bound) - y(l_bound))/h &
    1615           0 :                       - d1*d2y_dx2(alpha,l_bound) + d2*d2y_dx2(alpha,u_bound)
    1616             : 
    1617             : !     first term sum_a u_ai * ( p_ai + n_i * dpai/dqi * dqi/dni ). the factor n(i) is already
    1618             : !     contained in dq0_dn
    1619             : 
    1620           0 :             v_nl(i) = v_nl(i) + u_a(i,alpha) * (p_a + dpa_dq0 * dq0_dn(i)  )
    1621             : 
    1622             : !     the following sum we need for h_j which will be fourier transformed later. The IF
    1623             : !     condition excludes the contributions of high q values.
    1624             : 
    1625           0 :             grad_n_abs = sqrt(grad_n(i,1)**2.0 + grad_n(i,2)**2.0 + grad_n(i,3)**2.0)
    1626             : 
    1627           0 :             if ( q_0(i) .ne. q_cut ) then
    1628             : 
    1629             :                h_j(i,:) = h_j(i,:) + u_a(i,alpha)*cmplx((grad_n(i,:) / grad_n_abs) &
    1630           0 :                                    * dpa_dq0 * dq0_dgrad_n(i), 0.0)
    1631             : 
    1632             :             endif
    1633             :          enddo
    1634             :       enddo
    1635             : 
    1636             : !     now fourier transform h_j and carry out the sum over G.
    1637             : 
    1638           0 :       do i = 1,3
    1639             : 
    1640           0 :           call inplfft( h_j(:,i), -1 )
    1641             : 
    1642             :       enddo
    1643             : 
    1644           0 :       do i = 1, n_gvectors
    1645             : 
    1646           0 :          ind = G_ind(i)
    1647             : 
    1648           0 :          h_j(ind,:) = (0.0,1.0) * G(i,:) * h_j(ind,:)
    1649             : 
    1650             :       enddo
    1651             : 
    1652             : !     back fourier transform h_j and add it to the potential
    1653             : 
    1654           0 :       do i = 1,3
    1655             : 
    1656           0 :          call inplfft( h_j(:,i), 1 )
    1657             : 
    1658             :       enddo
    1659             : 
    1660           0 :       v_nl(:) = v_nl(:) - real( h_j(:,1) + h_j(:,2) + h_j(:,3) )
    1661             : 
    1662           0 :       deallocate(h_j, y,d2y_dx2)
    1663             : 
    1664           0 :       end subroutine
    1665             : !
    1666             : !==========================================================================================
    1667             : !
    1668             : #endif
    1669             : 
    1670           0 :       SUBROUTINE timing ( time )
    1671             :       USE param,ONLY: jnlout
    1672             : !
    1673             :       implicit none
    1674             : !
    1675             :       integer, intent(in)     :: time(8)
    1676             :       integer                 :: tmp(8), time_out(8)
    1677             : 
    1678             :       integer                 :: i
    1679             : 
    1680             : !     In the older versions of this program the time written to output file could be negative. This
    1681             : !     subroutine will fix this problem.
    1682             : 
    1683           0 :       tmp = 0
    1684             : 
    1685           0 :       tmp(8) = 1000
    1686           0 :       tmp(7) = 60
    1687           0 :       tmp(6) = 60
    1688           0 :       tmp(5) = 24
    1689             : 
    1690           0 :       do i=1,8
    1691             : 
    1692           0 :          if (time(i) .lt. 0) then
    1693           0 :             time_out(i) = time(i) + tmp(i)
    1694           0 :             time_out(i-1) = time(i-1) - 1
    1695             :          else
    1696           0 :             time_out(i) = time(i)
    1697             :          endif
    1698             : 
    1699             :       enddo
    1700             : !
    1701           0 :       write(jnlout, '(I4,A,I4,A,I4,A,I4,A)') time_out(5)," h  ", &
    1702           0 :                                              time_out(6)," min", &
    1703           0 :                                              time_out(7)," s  ", &
    1704           0 :                                              time_out(8)," ms "
    1705             : !
    1706           0 :       END SUBROUTINE timing
    1707             : !
    1708             : !     As the hole thing is intended to be interfaced with different codes I put here just an
    1709             : !     example program for testing, which is capable of reading an input file, reading a
    1710             : !     charge density, which otherwise has to be given by the calling program, and then call
    1711             : !     the subroutine soler(n, ... ).
    1712             : !
    1713             :       END MODULE nonlocal_funct

Generated by: LCOV version 1.13