1!!@LICENSE 2! 3!****************************************************************************** 4! MODULE m_vdwxc 5! Implements the nonlocal correlation energy part of the van der Waals density 6! functional of M.Dion et al, PRL 92, 246401 (2004): 7! Enlc = (1/2) Int Int dr1 dr2 rho(r1) phi(q1*r12,q2*r12) rho(r2) 8! with r12=|r2-r1|, q1=q0(rho(r1),grad_rho(r1)), q2=q0(rho(r2),grad_rho(r2)), 9! where rho(r) is the electron density at point r, grad_rho(r) its gradient, 10! and q0(rho,grad_rho) and phi(d1,d2) are universal functions defined in 11! Eqs.(11-12) and (14-16) of Dion et al. 12! Using separate module m_vv_vdwxc, it implements also the similar functional 13! of Vydrov and Van Voorhis. See that module for more information. 14! Refs: M.Dion et al, PRL 92, 246401 (2004) 15! J.Klimes et al, JPCM 22, 022201 (2009) 16! V.R.Cooper, PRB 81, 161104(R) (2010) 17! K.Lee et al, PRB 82, 081101 (2010) 18! O.A.Vydrov and T.VanVoorhis, JCP 133, 244103 (2010) 19! G.Roman-Perez and J.M.Soler, PRL 103, 096102 (2009) 20! Written by J.M.Soler. July 2007 - April 2010 and July 2012. 21!------------------------------------------------------------------------------ 22! Used module procedures: 23! use sys, only: die ! termination routine 24! use mesh1D, only: derivative ! Derivative of a function in a mesh 25! use m_ggaxc only: ggaxc ! General GGA XC routine 26! use m_ldaxc, only: ldaxc ! General LDA XC routine 27! use mesh1D, only: integral ! Integral of a function in a mesh 28! use mesh1D, only: get_mesh ! Returns the mesh points 29! use mesh1D, only: get_n ! Returns the number of mesh points 30! use m_radfft, only: radfft ! Radial fast Fourier transform 31! use mesh1D, only: set_interpolation ! Sets interpolation method 32! use mesh1D, only: set_mesh ! Sets a 1D mesh 33! use interpolation, only: spline ! Sets spline interpolation 34! use interpolation, only: splint ! Calculates spline interpolation 35! use m_vv_vdwxc, only: vv_vdw_beta ! Parameter beta of VV2010 functionl 36! use m_vv_vdwxc, only: vv_vdw_theta ! Func. theta of VV2010 functional 37! use m_vv_vdwxc, only: vv_vdw_get_kmesh ! Size and values of (kf,kg) mesh 38! use m_vv_vdwxc, only: vv_vdw_phi ! Interpolate rho1*rho2*phi(k1,k2,k) 39! use m_vv_vdwxc, only: vv_vdw_set_kcut ! Sets the planewave cutoff kc 40!------------------------------------------------------------------------------ 41! Used module parameters: 42! use precision, only: dp ! Real double precision type 43!------------------------------------------------------------------------------ 44! Public procedures available from this module: 45! vdw_decusp : Energy due to the softening of the VdW kernel cusp 46! vdw_get_qmesh : Returns size and values of q-mesh 47! vdw_localxc : LDA/GGA xc part apropriate for the used vdW flavour 48! vdw_phi : Finds and interpolates phi(q1,q2,k) 49! vdw_set_author: Sets the vdW functional flavour (author initials) 50! vdw_set_kcut : Sets the planewave cutoff kc of the integration grid 51! vdw_theta : Finds function theta_q(rho,grad_rho) 52!------------------------------------------------------------------------------ 53! Public types, parameters, variables, and arrays: 54! None 55!------------------------------------------------------------------------------ 56! Units: 57! All lengths and energies in atomic units (Bohr, Hartree) 58!****************************************************************************** 59! subroutine vdw_decusp( nspin, rhos, grhos, eps, dedrho, dedgrho ) 60! Finds the local energy correction due to the softening of the VdW kernel 61! cusp, defined as eps(rho,grad_rho) = 62! (1/2) * rho * Int 4*pi*r**2*dr * (phi(q1*r,q2*r) - phi_soft(q1*r,q2*r)) 63! where q1=q2=q0(rho,grad_rho). Notice that grad_rho is included in the value 64! of q0 at the origin but not in the change of rho(r) in the integrand. 65! phi_soft(d1,d2) is a softened version of the nonlocal VdW kernel phi(d1,d2) 66! (Eq.(14) of Dion et al) in which the logarithmic divergence at d1=d2=0 is 67! substituted by a smooth analytic function of the form defined in vdw_phi 68! Arguments: 69! integer, intent(in) :: nspin ! Number of spin components 70! real(dp),intent(in) :: rhos(nspin) ! Electron spin density 71! real(dp),intent(in) :: grhos(3,nspin) ! Spin density gradient 72! real(dp),intent(out):: eps ! Energy correction per electron 73! real(dp),intent(out):: dedrho(nspin) ! d(rho*eps)/d(rho) 74! real(dp),intent(out):: dedgrho(3,nspin) ! d(rho*eps)/d(grad_rho) 75! Notes: 76! - Requires a previous call to vdw_set_kcut. Otherwise stops with an error msg. 77!------------------------------------------------------------------------------ 78! subroutine vdw_exchng( iRel, nSpin, D, GD, epsX, dEXdD, dEXdGD ) 79! Finds the exchange energy density and its derivatives, using the GGA 80! functional apropriate for the previously-set vdW functional flavour 81! Arguments: 82! integer, intent(in) :: iRel ! Relativistic exchange? 0=no, 1=yes 83! integer, intent(in) :: nSpin ! Number of spin components 84! real(dp),intent(in) :: D(nSpin) ! Local electron (spin) density 85! real(dp),intent(in) :: GD(3,nSpin) ! Gradient of electron density 86! real(dp),intent(out):: epsX ! Exchange energy per electron 87! real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 88! real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 89! Sample usage: 90! integer,parameter:: iRel=0, nSpin=2 91! real(dp):: D(nSpin), dEXdD(nSpin), dEXdGD(3,nSpin), epsX, GD(3,nSpin) 92! call vdw_set_author('DRSLL') 93! do r points 94! Find D and GD at r 95! call vdw_exchng( iRel, nSpin, D, GD, epsX, dEXdD, dEXdGD ) 96! Ex = Ex + dVolume*sum(D)*epsX 97! end do 98!------------------------------------------------------------------------------ 99! subroutine vdw_get_qmesh( n, q ) 100! Returns size and values of q-mesh 101! Arguments: 102! integer, intent(out) :: n ! Number of q mesh points 103! real(dp),optional,intent(out) :: q(:) ! Values of q mesh points 104! Sample usage: 105! integer :: nq 106! real(dp):: kcut 107! real(dp),allocatable:: qmesh(:) 108! kcut = 10._dp ! 10 Bohr^-1 => 100 Ry (this should be adapted to your mesh) 109! call vdw_set_kcut( kcut ) 110! call vdw_get_qmesh( nq ) 111! allocate( qmesh(nq) ) 112! call vdw_get_qmesh( nq, qmesh ) 113! Notes: 114! - Requires a previous call to vdw_set_kcut. Otherwise stops with an error msg. 115! - If the size of array q is smaller than that of the stored qmesh, it is 116! filled with the first size(q) values of qmesh 117! - The size and values of the logarithmic q mesh are set by internal 118! parameters that can be changed only by editing them in this module: 119! nq : Number of q mesh points 120! qcut=qmesh(nq) : Max. value of q mesh 121! dqmaxdqmin : (qmesh(nq) - qmesh(nq-1)) / (qmesh(2) - qmesh(1)) 122! Although the presently-set values have been found to yield good accuracy 123! in preliminary calculations, more tests may be required to guarantee 124! convergence in other systems. The value of nq is particularly important: 125! larger nq increases accuracy but CPU time increases between nq and nq**2. 126!------------------------------------------------------------------------------ 127! subroutine vdw_phi( k, phi, dphidk ) 128! Finds phi_soft(q1,q2,k) (Fourier transform of phi_soft(q1,q2,r)) for all 129! values of q1 and q2 in qmesh. The q's are local wavevectors defined in 130! eqs.(11-12), and phi_soft is a smoothed version of Eq.(14) of Dion et al, 131! defined as phi_soft(d1,d2) = phi_soft(d,a) = phi0 + phi2*d**2 + phi4*d**4, 132! where phi0 is a parameter, d=sqrt(d1**2+d2**2), a=atan(d2/d1), and phi2, 133! phi4 are chosen so that phi_soft(d,a) matches phi(d,a) in value and slope 134! at d=dsoft (another parameter). 135! Arguments: 136! real(dp),intent(in) :: k ! Modulus of actual k vector 137! real(dp),intent(out):: phi(:,:) ! phi(q1,q2,k) at given k 138! ! for all q1,q2 in qmesh 139! real(dp),intent(out):: dphidk(:,:) ! dphi(q1,q2,k)/dk at given k 140! Sample usage: 141! integer :: nq 142! real(dp):: k, kcut 143! real(dp),allocatable:: phi(:,:), dphidk(:,:) 144! kcut = 10._dp ! 10 Bohr^-1 => 100 Ry (this should be adapted to your mesh) 145! call vdw_set_kcut( kcut ) 146! call vdw_get_qmesh( nq ) 147! allocate( phi(nq,nq), dphidk(nq,nq) ) 148! do k points 149! call vdw_phi( k, phi, dphidk ) 150! Notes: 151! - Requires a previous call to vdw_set_kcut. Otherwise stops with an error msg. 152! - Stops with an error message if size of array phi is smaller than nq*nq. 153!----------------------------------------------------------------------------- 154! subroutine vdw_set_author( author ) 155! Sets the functional flavour (author initials) and subsequent parameters 156! Arguments: 157! character(len=*),intent(in):: author ! Functnl flavour 158! ('DRSLL'|'LMKLL'|'KBM'|'C09'|'BH'|'VV') 159! Notes: 160! - If vdw_set_author is not called, author='DRSLL' is set by default 161! - Stops with an error message if author has not an allowed value 162!----------------------------------------------------------------------------- 163! subroutine vdw_set_kcut( kc ) 164! Sets the reciprocal planewave cutoff kc of the integration grid, and finds 165! the interpolation table to be used by vdw_phi to obtain the vdW kernel phi 166! at the reciprocal mesh wavevectors. 167! Arguments: 168! real(dp),intent(in):: kc ! Planewave cutoff: k>kcut => phi=0 169! Notes: 170! - An interpolation table to calculate the VdW kernel phi is read from file 171! 'vdw_kernel.table'. If the file does not exist, the table is calculated 172! and the file written when vdw_set_kcut is called. 173!------------------------------------------------------------------------------ 174! subroutine vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho ) 175! Finds the value and derivatives of theta_i(rho,grad_rho) = rho*p_i(q0), 176! where q0(rho,grad_rho) is the local wavevector defined in eqs.(11-12) 177! of Dion et al, PRL 92, 246401 (2004). 178! p_i(q0) are the cubic polynomials such that 179! y(q0) = Sum_i p_i(q0) * y_i 180! is the cubic spline interpolation at q0 of (any) function y(q) with 181! values y_i at mesh points qmesh_i 182! Arguments: 183! integer, intent(in) :: nspin ! Number of spin components 184! real(dp),intent(in) :: rhos(nspin) ! Electron spin density 185! real(dp),intent(in) :: grhos(3,nspin) ! Spin density gradient 186! real(dp),intent(out):: theta(nq) ! Expansion of rho*q in qmesh 187! real(dp),intent(out):: dtdrho(nq,nspin) ! dtheta(iq)/drhos 188! real(dp),intent(out):: dtdgrho(3,nq,nspin) ! dtheta(iq)/dgrhos(ix) 189! Notes: 190! - Requires a previous call to vdw_set_kcut 191! - The value of q0(rho,grad_rho) is saturated at qmesh(nq)=qcut (a parameter) 192! to avoid incorrect 'interpolation' of phi(q1,q2,r12) from qmesh points. 193!****************************************************************************** 194! Algorithms: 195! Although Eqs.(14-16) of Dion et al provide a straightforward definition of 196! the nonlocal VdW kernel phi(d1,d2), its direct evaluation with the required 197! accuracy turns out to be too time consuming. In addition, the kernel has a 198! logarithmic divergence at d1=d2=0, which prevents its direct numerical 199! Fourier transform. Therefore, an elaborate layered procedure is followed: 200! A smoothed version of phi(d1,d2) is defined as 201! phi_soft(d1,d2) = phi_soft(d,a) = phi0 + phi2*d**2 + phi4*d**4, 202! where phi0 is a parameter, d=sqrt(d1**2+d2**2), a=atan(d2/d1), and phi2, 203! phi4 are chosen so that phi_soft(d,a) matches phi(d,a) in value and slope 204! at d=dsoft (another parameter). 205! The difference in energy between phi_soft and the right phi (called DEcusp) 206! is approximated in a local-density approximation as eps(rho,grad_rho) = 207! (1/2) * rho * Int 4*pi*r**2*dr * (phi(q1*r,q2*r) - phi_soft(q1*r,q2*r)) 208! where q1=q2=q0(rho,grad_rho). Notice that grad_rho is included in the value 209! of q0 at the origin but not in the change of rho(r) in the integrand. This 210! could be improved in the future, although preliminary tests suggest that 211! the corrections would be quite small. 212! A one-dimensional table of phi(d,d) is first calculated and stored as a 213! function of d. Then, phi(d1,d2) is calculated as phi(dmax,dmax)+dphi(d1,d2), 214! with dmax=max(d1,d2). The reason is that dphi converges better for large d, 215! and therefore it requires smaller integration limits in Eq.(14) of Dion et al, 216! while phi(dmax,dmax) is interpolated from the stored table. 217! Using the above methods, a two-dimensional table of phi_soft(d1,d2) is 218! calculated and stored for all values of d1 and d2 in a logarithmic dmesh 219! of size nd with a cutoff dcut (two internal parameters). This table is 220! written on disk file 'vdw_kernel.table' for reuse in subsequent runs. 221! If the file exists already, the table is simply read and stored in memory. 222! Once we set the cutoff kcut (expressed as a wavevector), of the integration 223! mesh to be used in evaluating the VdW nonlocal energy Enlc, we calculate and 224! store (in memory) another interpolation table of phi_soft(q1,q2,k) for all 225! values of q1 and q2 in qmesh, and of k in a fine radial grid of nk=nr points. 226! This is done by first calculating phi_soft(q1,q2,r)=phi_soft(q1*r,q2*r) in a 227! radial grid in real space and then Fourier-transforming it to k space. 228! This table is then used to interpolate phi_soft(q1,q2,k) for any desired k. 229! In order to ensure that values of q are within the interpolation range, 230! they are 'saturated' smoothly to a cutoff qcut=qmesh(nq). This implies an 231! approximation when either rho is very large (i.e. near the nucleus) or when 232! (grad_rho/rho)**2/kF -> infinity, what tipically occurs in the tails of 233! the electron density. In the first case, Enlc is neglegible compared with 234! Ex and the local part of Ec. In the second case, it is neglegible because of 235! the factor rho in the integrand. Thus, the preliminary tests suggest that the 236! presently-set value of qcut gives sufficiently accurate values. 237!****************************************************************************** 238 239MODULE m_vdwxc 240 241! Used module procedures 242 use sys, only: die ! termination routine 243 use mesh1D, only: derivative ! Derivative of a function in a mesh 244 use mesh1D, only: integral ! Integral of a function in a mesh 245 use mesh1D, only: get_mesh ! Returns the mesh points 246 use mesh1D, only: get_n ! Returns the number of mesh points 247 use m_ggaxc, only: ggaxc ! General GGA XC routine 248 use m_ldaxc, only: ldaxc ! General LDA XC routine 249 use m_radfft, only: radfft ! Radial fast Fourier transform 250 use alloc, only: re_alloc ! Re-allocation routine 251 use mesh1D, only: set_interpolation ! Sets interpolation method 252 use mesh1D, only: set_mesh ! Sets a 1D mesh 253 use interpolation,only: spline ! Sets spline interpolation 254 use interpolation,only: splint ! Calculates spline interpolation 255 use m_vv_vdwxc, only: vv_vdw_beta ! Parameter beta of VV2010 functional 256 use m_vv_vdwxc, only: vv_vdw_theta ! Func. theta of VV2010 functional 257 use m_vv_vdwxc, only: vv_vdw_get_kmesh ! Size and values of (kf,kg) mesh 258 use m_vv_vdwxc, only: vv_vdw_phi ! Interpolates rho1*rho2*phi(k1,k2,k) 259 use m_vv_vdwxc, only: vv_vdw_set_kcut ! Sets the planewave cutoff kc 260 261! Used module parameters 262 use precision, only: dp ! Real double precision type 263 264#ifdef DEBUG_XC 265 use m_vv_vdwxc, only: vv_vdw_phiofr ! Interpolates rho1*rho2*phi(k1,k2,r) 266 use debugXC, only: udebug ! File unit for debug output 267! use plot_module, only: plot 268#endif /* DEBUG_XC */ 269 270 implicit none 271 272! Called by xc routines 273PUBLIC :: & 274 vdw_decusp, &! Energy due to the softening of the VdW kernel cusp 275 vdw_get_qmesh, &! Returns size and values of q-mesh 276 vdw_localxc, &! LDA/GGA exchange-corr apropriate for the used vdW flavour 277 vdw_phi, &! Finds and interpolates phi(q1,q2,k) 278 vdw_set_author, &! Sets the vdW functional flavour (author initials) 279 vdw_set_kcut, &! Sets the planewave cutoff kc of the integration grid 280 vdw_theta ! Finds function theta_q(rho,grad_rho) 281 282#ifdef DEBUG_XC 283! Called by debugging test programs 284PUBLIC :: & 285 phiofr, &! Finds the kernel phi(q1,q2,r) at tabulated q-mesh values 286 phi_interp, &! Finds soft phi(d1,d2) kernel by interpolation of phi_table 287 phi_soft, &! Finds phi(d1,d2) softened near d1=d2=0 288 phi_val, &! Finds kernel phi(d1,d2) by direct integration 289 pofq, &! Finds polynomials p(q) from cubic spline interpolation 290 qofrho ! Finds the local wavevector parameter q(rho,grad_rho) 291#endif /* DEBUG_XC */ 292 293PRIVATE ! Nothing is declared public beyond this point 294 295! integer, parameter:: dp = kind(1.d0) 296 297 ! Precision parameters for the integral defining phi in routine phi_val 298 real(dp),parameter:: acutmin = 10.0_dp ! Upper integration limit 299 real(dp),parameter:: acutbyd = 30.0_dp ! Upper integration limit / d 300 real(dp),parameter:: damax = 0.5_dp ! Max. integration interval 301 real(dp),parameter:: damin = 1.e-2_dp ! Min. integration interval 302 real(dp),parameter:: dabyd = 0.1_dp ! Min. integration interval / d 303 304 ! Precision parameter for the integral in routine dphi 305 real(dp),parameter:: ashortbyd = 2.5_dp ! Shorter integration limit / d 306 307 ! Parameters for phi_soft. Some recommended pairs of values are 308 ! (dsoft,phi0_soft)=(0.5,0.8)|(0.7|0.6)|(1.0|0.4)|(1.5,0.22)|(2.0,0.12) 309 real(dp),parameter:: dsoft = 1.0_dp ! Softening matching radius 310 real(dp),parameter:: phi0_soft = 0.40_dp ! phi_soft(0,0) (depends on dsoft) 311 real(dp),parameter:: dd = 0.01_dp ! Delta(d) for derivatives 312 313 ! Mesh parameters for table phi(d1,d2) 314 integer, parameter:: nd = 20 ! Number of d mesh points 315 real(dp),parameter:: dcut = 30.0_dp ! Max. value of d mesh 316 real(dp),parameter:: ddmaxddmin = 20.0_dp ! Last d mesh interval / first one 317 318 ! Use routine dphi for better efficiency in setting phi table? 319 logical,parameter:: use_dphi =.true. !(.true.=>efficiency|.false.=>accuracy) 320 321 ! Set derivation methods to use for interpolation table 322 character(len=*),parameter:: deriv_method = 'numeric' !('numeric'|'interp') 323 character(len=*),parameter:: interp_method= 'Spline' !('Lagrange'|'Spline') 324 325 ! Parameters to find numerical derivatives of phi, to set interpolation 326 real(dp),parameter:: ddbyd = 0.01_dp ! Delta to find phi derivs / d 327 real(dp),parameter:: ddmin = 0.001_dp ! Min. delta to find phi derivs 328 329 ! Mesh parameters for table of phi(q1,q2,r) and its Fourier transform 330 integer, parameter:: nr = 2048 ! Radial mesh points (power of 2) 331 integer, parameter:: mq = 30 ! Total number of q mesh points 332! integer, parameter:: nq = mq-1 ! Effective number of q mesh points 333 real(dp),parameter:: qcut = 5.0_dp ! Max. value of q mesh 334 real(dp),parameter:: dqmaxdqmin = 20.0_dp ! Last q mesh interval / first one 335 real(dp),parameter:: rcut = 100._dp ! Radial cutoff: r>rcut => phi=0 336 real(dp),parameter:: rmin = 1.e-6_dp ! Min. radius as denominator 337 real(dp),parameter:: rsoft = 0.0_dp ! Soften kernel in r<rsoft 338 339 ! Parameters for cutoff function, used in radial Fourier transforms of phi 340 integer, parameter:: ncut1 = 8 ! cutoff(x)=(1-x**ncut1)**ncut2 341 integer, parameter:: ncut2 = 4 342 343 ! Parameters for saturate function, used to enforce that q<qcut 344 integer, parameter:: nsat = 12 ! xsat(x)=1-exp(-sum_n=1:nsat x**n/n) 345 346 ! Parameters for saturate_inverse function 347 real(dp),parameter:: xmaxbyxc = 1.5_dp ! qmax/qcut 348 real(dp),parameter:: ytol = 1.e-15_dp ! Tol. for saturated q 349 350 ! Private module variables and arrays 351 character(len=5),save:: vdw_author='DRSLL' ! Functional 'flavour' name 352 real(dp),save:: dmesh(nd) ! Mesh points for phi(d1,d2) table 353 real(dp),save:: qmesh(mq) ! Mesh points for phi(q1,q2,r) 354 real(dp),save:: phi_table(0:3,0:3,nd,nd) ! Coefs. for bicubic interpolation 355 logical, save:: phi_table_set=.false. ! Has phi_table been set? 356 logical, save:: qmesh_set=.false. ! Has qmesh been set? 357 logical, save:: kcut_set=.false. ! Has kcut been set? 358 real(dp),save:: dr ! r-mesh interval 359 real(dp),save:: dk ! k-mesh interval 360 real(dp),save:: kcut ! Planewave cutoff: k>kcut => phi=0 361 integer, save:: nk ! # k points within kcut 362 real(dp),save:: zab=-0.8491_dp ! Parameter of the vdW functional 363 real(dp),pointer,save:: & 364 phir(:,:,:)=>null(), &! Table of phi(r) 365 phik(:,:,:)=>null(), &! Table of phi(k) 366 d2phidr2(:,:,:)=>null(),&! Table of d^2(phi)/dr^2 367 d2phidk2(:,:,:)=>null() ! Table of d^2(phi)/dk^2 368 369! real(dp),save:: dqmaxdqmin, qcut 370 371CONTAINS 372 373! ----------------------------------------------------------------------------- 374 375SUBROUTINE bcucof( n1, n2, x1, x2, y, dydx1, dydx2, d2ydx1dx2, c ) 376! Finds coefficients for bicubic interpolation 377! Adapted from Numerical recipes 378 IMPLICIT NONE 379 INTEGER, INTENT(IN) :: n1, n2 380 REAL(dp),INTENT(IN) :: x1(n1) 381 REAL(dp),INTENT(IN) :: x2(n2) 382 REAL(dp),INTENT(IN) :: y(n1,n2) 383 REAL(dp),INTENT(IN) :: dydx1(n1,n2) 384 REAL(dp),INTENT(IN) :: dydx2(n1,n2) 385 REAL(dp),INTENT(IN) :: d2ydx1dx2(n1,n2) 386 REAL(dp),INTENT(OUT):: c(0:3,0:3,n1,n2) 387 388 INTEGER :: i1, i11, i12, i13, i14, i2, i21, i22, i23, i24 389 REAL(dp) :: dx1, dx2, wt(16,16), z(16) 390 DATA wt /1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,& 391 8*0,3,0,-9,6,-2,0,6,-4,10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,& 392 2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,& 393 2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,0,1,-2,1,5*0,& 394 -3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,10*0,-3,3,2*0,2,-2,2*0,& 395 -1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,& 396 -1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/ 397 398! Set coefs. for i1<n1 and i2<n2 399 do i2 = 1,n2-1 400 do i1 = 1,n1-1 401 dx1 = x1(i1+1) - x1(i1) 402 dx2 = x2(i2+1) - x2(i2) 403 i11 = i1 404 i12 = i1+1 405 i13 = i1+1 406 i14 = i1 407 i21 = i2 408 i22 = i2 409 i23 = i2+1 410 i24 = i2+1 411 z( 1) = y(i11,i21) 412 z( 2) = y(i12,i22) 413 z( 3) = y(i13,i23) 414 z( 4) = y(i14,i24) 415 z( 5) = dydx1(i11,i21) * dx1 416 z( 6) = dydx1(i12,i22) * dx1 417 z( 7) = dydx1(i13,i23) * dx1 418 z( 8) = dydx1(i14,i24) * dx1 419 z( 9) = dydx2(i11,i21) * dx2 420 z(10) = dydx2(i12,i22) * dx2 421 z(11) = dydx2(i13,i23) * dx2 422 z(12) = dydx2(i14,i24) * dx2 423 z(13) = d2ydx1dx2(i11,i21) * dx1 * dx2 424 z(14) = d2ydx1dx2(i12,i22) * dx1 * dx2 425 z(15) = d2ydx1dx2(i13,i23) * dx1 * dx2 426 z(16) = d2ydx1dx2(i14,i24) * dx1 * dx2 427 z = matmul(wt,z) 428 c(0:3,0:3,i1,i2) = reshape(z,(/4,4/),order=(/2,1/)) 429 end do ! i1 430 end do ! i2 431 432! Set c for i1=n1 and i2=n2 (valid only at the very border) using 433! sum_i,j c(i,j,n1,i2) * 0**i * x2**j = sum_i,j c(i,j,n1-1,i2) * 1**i * x2**j 434! sum_i,j c(i,j,i1,n2) * x1**i * 0**j = sum_i,j c(i,j,i1,n2-1) * x1**i * 1**j 435! sum_i,j c(i,j,n1,n2) * 0**i * 0**j = sum_i,j c(i,j,n1-1,n2-1) * 1**i * 1**j 436 c(:,:,n1,:) = 0 437 c(:,:,:,n2) = 0 438 c(0,:,n1,1:n2-1) = sum(c(:,:,n1-1,1:n2-1),dim=1) 439 c(:,0,1:n1-1,n2) = sum(c(:,:,1:n1-1,n2-1),dim=2) 440 c(0,0,n1,n2) = sum(c(:,:,n1-1,n2-1)) 441 442END SUBROUTINE bcucof 443 444! ----------------------------------------------------------------------------- 445 446function cutoff( x ) 447 448 implicit none 449 real(dp),intent(in):: x 450 real(dp) :: cutoff 451 452 if (x<=0._dp) then 453 cutoff = 1 454 else if (x>=1._dp) then 455 cutoff = 0 456 else 457 cutoff = (1-x**ncut1)**ncut2 458 end if 459 460end function cutoff 461 462!----------------------------------------------------------------------------- 463 464real(dp) function dphi( d1, d2 ) 465 466! Finds phi(d1,d2) - phi(dmax,dmax), with dmax=max(d1,d2) 467 468 real(dp),intent(in) :: d1, d2 469 470 integer :: ia, ib, n, nmesh, nshort 471 real(dp) :: a, acut, b, da1, dan, deff, dmax, dmin, gamma, & 472 pi, t, t0, w 473 real(dp),allocatable :: amesh(:), c(:), dphida(:), dphidb(:), & 474 s(:), nu0(:), nu1(:), nu2(:) 475 476! Find integration mesh 477 dmax = max( abs(d1), abs(d2) ) 478 dmin = min( abs(d1), abs(d2) ) 479 deff = dmax 480! deff = sqrt(d1**2+d2**2) 481 acut = max( acutmin, acutbyd * deff ) 482 da1 = min( damax, dabyd * deff ) 483 da1 = max( da1, damin ) 484 dan = damax 485 n = get_n( 0._dp, acut, da1, dan ) 486 allocate( amesh(n), c(n), dphida(n), dphidb(n), s(n), & 487 nu0(n), nu1(n), nu2(n) ) 488 call set_mesh( n, xmax=acut, dxndx1=dan/da1 ) 489 call get_mesh( n, nmesh, amesh ) 490 491! Find limit of shorter mesh 492 nshort = n 493 do ia = n-1,1,-1 494 if (amesh(ia) > ashortbyd*deff) nshort = ia 495 end do 496 497! Find cos(a), sin(a), nu1(a), and nu2(a) 498 pi = acos(-1._dp) 499 gamma = 4*pi/9 500 do ia = 1,n 501 a = amesh(ia) 502 c(ia) = cos(a) 503 s(ia) = sin(a) 504 if (a==0._dp) then 505 nu0(ia) = dmax**2 / 2 / gamma 506 nu1(ia) = d1**2 / 2 / gamma 507 nu2(ia) = d2**2 / 2 / gamma 508 else ! (a/=0) 509 if (d1==0._dp) then 510 nu1(ia) = a*a / 2 511 else 512 nu1(ia) = a*a / 2 / (1-exp(-gamma*(a/d1)**2)) 513 end if 514 if (d2==0._dp) then 515 nu2(ia) = a*a / 2 516 else 517 nu2(ia) = a*a / 2 / (1-exp(-gamma*(a/d2)**2)) 518 end if 519 if (dmax<=0._dp) then 520 nu0(ia) = a*a / 2 521 else 522 nu0(ia) = a*a / 2 / (1-exp(-gamma*(a/dmax)**2)) 523 end if 524 end if ! (a==0) 525 end do 526 527! Make integral on variable a 528 dphida(1) = 0 529 do ia = 2,nshort 530 a = amesh(ia) 531 532 ! Make integral on variable b 533 dphidb(1) = 0 534 do ib = 2,n 535 b = amesh(ib) 536 537 w = 2*( (3-a*a)*b*c(ib)*s(ia) + (3-b*b)*a*c(ia)*s(ib) & 538 + (a*a+b*b-3)*s(ia)*s(ib) - 3*a*b*c(ia)*c(ib) )/(a*b)**3 539 540 t = 0.5_dp * ( 1/(nu1(ia)+nu1(ib)) + 1/(nu2(ia)+nu2(ib)) ) & 541 * ( 1/(nu1(ia)+nu2(ia))/(nu1(ib)+nu2(ib)) & 542 + 1/(nu1(ia)+nu2(ib))/(nu2(ia)+nu1(ib)) ) 543 544 t0 = 0.5_dp * ( 1/(nu0(ia)+nu0(ib)) + 1/(nu0(ia)+nu0(ib)) ) & 545 * ( 1/(nu0(ia)+nu0(ia))/(nu0(ib)+nu0(ib)) & 546 + 1/(nu0(ia)+nu0(ib))/(nu0(ia)+nu0(ib)) ) 547 548 dphidb(ib) = a*a * b*b * w * (t-t0) 549 550 end do ! ib 551 dphida(ia) = integral( n, dphidb ) - integral( ia, dphidb ) 552 553 end do ! ia 554 dphi = 2/pi**2 * 2*integral( nshort, dphida ) 555 556 deallocate( amesh, c, dphida, dphidb, s, nu0, nu1, nu2 ) 557 558#ifdef DEBUG_XC 559! print'(a,2f12.6,i8,f12.6)', 'dphi: d1,d2, na, phi =', d1, d2, n, dphi 560#endif /* DEBUG_XC */ 561 562end function dphi 563 564! ----------------------------------------------------------------------------- 565 566function dphi_fast( d1, d2 ) 567 568! Finds phi(d1,d2)-phi(dmax,dmax), with dmax=max(d1,d2), by 569! - Direct integration if d < dsoft, where d=sqrt(d1**2+d2**2) 570! - Interpolation of phi_table if d > dsoft 571 572 implicit none 573 real(dp),intent(in) :: d1, d2 574 real(dp) :: dphi_fast 575 576 real(dp):: d, dmax 577 578 if (.not.phi_table_set) call set_phi_table() 579 580 d = sqrt( d1**2 + d2**2 ) 581 dmax = max( d1, d2 ) 582 583 if (d < dsoft) then 584 dphi_fast = dphi( d1, d2 ) 585 else 586 dphi_fast = phi_interp( d1, d2 ) - phi_interp( dmax, dmax ) 587 end if 588 589end function dphi_fast 590 591! ----------------------------------------------------------------------------- 592 593function dphi_soft( d1, d2 ) 594 595! Finds phi_soft(d1,d2)-phi_soft(dmax,dmax), with dmax=max(d1,d2) 596 597 implicit none 598 real(dp),intent(in) :: d1, d2 599 real(dp) :: dphi_soft 600 601 real(dp):: d, dmax 602 603 d = sqrt( d1**2 + d2**2 ) 604 dmax = max( d1, d2 ) 605 606 if (d < dsoft) then 607 dphi_soft = phi_soft( d1, d2 ) - phi_soft( dmax, dmax ) 608 else 609 dphi_soft = dphi( d1, d2 ) 610 end if 611 612end function dphi_soft 613 614!----------------------------------------------------------------------------- 615 616integer function iofd( d ) 617 618! Finds index i such that dmesh(i) <= d < dmesh(i+1) 619 620 implicit none 621 real(dp), intent(in) :: d 622 623 real(dp),parameter :: amin = 1.e-12_dp 624 real(dp),save:: a, b 625 logical, save:: first_call = .true. 626 627 if (first_call) then 628 a = log( (dmesh(nd)-dmesh(nd-1)) / (dmesh(2)-dmesh(1)) ) / (nd-2) 629 a = max( a, amin ) 630 b = (dmesh(2) - dmesh(1)) / (exp(a) - 1) 631 first_call = .false. 632 end if 633 634 iofd = 1 + log( 1 + (d-dmesh(1))/b ) / a 635 iofd = max( 1, iofd ) 636 iofd = min( nd-1, iofd ) 637 638end function iofd 639 640!----------------------------------------------------------------------------- 641 642integer function iofq( q ) 643 644! Finds index i such that qmesh(i) <= q < qmesh(i+1) 645 646 implicit none 647 real(dp), intent(in) :: q 648 649 real(dp),parameter :: amin = 1.e-12_dp 650 real(dp),save:: a, b 651 logical, save:: first_call = .true. 652 653 if (first_call) then 654 a = log( (qmesh(mq)-qmesh(mq-1)) / (qmesh(2)-qmesh(1)) ) / (mq-2) 655 a = max( a, amin ) 656 b = (qmesh(2) - qmesh(1)) / (exp(a) - 1) 657 first_call = .false. 658 end if 659 660 iofq = 1 + log( 1 + (q-qmesh(1))/b ) / a 661 iofq = max( 1, iofq ) 662 iofq = min( mq-1, iofq ) 663 664end function iofq 665 666! ----------------------------------------------------------------------------- 667 668function phi_fast( d1, d2 ) 669 670! Finds hard phi(d1,d2) kernel by 671! - Direct integration if d < dsoft, where d=sqrt(d1**2+d2**2) 672! - Interpolation of phi_table if d > dsoft 673 674 implicit none 675 real(dp),intent(in) :: d1, d2 676 real(dp) :: phi_fast 677 678 real(dp):: d 679 680 if (.not.phi_table_set) call set_phi_table() 681 682 d = sqrt( d1**2 + d2**2 ) 683 684 if (d < dsoft) then 685 phi_fast = phi_val( d1, d2 ) 686 else 687 phi_fast = phi_interp( d1, d2 ) 688 end if 689 690end function phi_fast 691 692! ----------------------------------------------------------------------------- 693 694function phi_interp( d1, d2 ) 695 696! Finds soft phi(d1,d2) kernel by interpolation of phi_table 697 698 implicit none 699 real(dp),intent(in) :: d1, d2 700 real(dp) :: phi_interp 701 702 integer :: i1, i2, id1, id2 703 real(dp):: dd1(0:3), dd2(0:3) 704 705 if (.not.phi_table_set) call set_phi_table() 706 707 if (d1>=dcut .or. d2>=dcut) then 708 phi_interp = 0 709 return 710 end if 711 712 id1 = iofd( d1 ) 713 id2 = iofd( d2 ) 714 715 dd1(0) = 1 716 dd2(0) = 1 717 dd1(1) = (d1 - dmesh(id1)) / (dmesh(id1+1) - dmesh(id1)) 718 dd2(1) = (d2 - dmesh(id2)) / (dmesh(id2+1) - dmesh(id2)) 719 dd1(2) = dd1(1)**2 720 dd2(2) = dd2(1)**2 721 dd1(3) = dd1(1)**3 722 dd2(3) = dd2(1)**3 723 724 phi_interp = 0 725 do i2 = 0,3 726 do i1 = 0,3 727 phi_interp = phi_interp + phi_table(i1,i2,id1,id2) * dd1(i1) * dd2(i2) 728 end do 729 end do 730 731end function phi_interp 732 733! ----------------------------------------------------------------------------- 734 735subroutine phiofr( r, phi ) 736 737! Finds phi(q1,q2,r) = phi(q1*r,q2*r) with q1=qmesh(iq1), q2=qmesh(iq2), 738! by interpolation from phi_table 739! Notice that q is a density parameter, related to the Fermi wavevector 740 741 implicit none 742 real(dp),intent(in) :: r 743 real(dp),intent(out):: phi(:,:) 744 745 integer :: iq1, iq2 746 real(dp):: dphidr 747 748 ! Trap VV-version exception 749#ifdef DEBUG_XC 750 if (vdw_author=='VV') then 751 call vv_vdw_phiofr( r, phi ) 752 return 753 end if 754#endif /* DEBUG_XC */ 755 756 if (size(phi,1)<mq .or. size(phi,2)<mq) & 757 stop 'phiofr: ERROR: size(phi) too small' 758 if (.not.qmesh_set) call set_qmesh() 759 760 if (r >= rcut) then 761 phi(:,:) = 0 762 else 763 do iq2 = 1,mq 764 do iq1 = 1,iq2 765 ! Use unfiltered kernel 766! phi(iq1,iq2) = phi_interp( qmesh(iq1)*r, qmesh(iq2)*r ) 767 ! Use filtered kernel 768 call splint( dr, phir(:,iq1,iq2), d2phidr2(:,iq1,iq2), nr+1, r, & 769 phi(iq1,iq2), dphidr ) 770 phi(iq2,iq1) = phi(iq1,iq2) 771 end do ! iq1 772 end do ! iq2 773 end if ! (r>=rcut) 774 775end subroutine phiofr 776 777!----------------------------------------------------------------------------- 778 779function phi_soft( d1, d2 ) 780 781! Finds phi(d1,d2) softened near d1=d2=0 782 783 implicit none 784 real(dp),intent(in) :: d1, d2 785 real(dp) :: phi_soft 786 787 real(dp):: d, d1m, d1p, d2m, d2p, dphidd, & 788 phi, phi0, phi2, phi4, phim, phip 789 790 d = sqrt( d1**2 + d2**2 ) 791 792 if (d<=0._dp) then 793 phi_soft = phi0_soft 794 else if (d > dsoft) then 795 phi_soft = phi_val( d1, d2 ) 796 else ! (0<d<dsoft) 797 798 d1p = d1/d * (dsoft+dd) 799 d2p = d2/d * (dsoft+dd) 800 d1m = d1/d * (dsoft-dd) 801 d2m = d2/d * (dsoft-dd) 802 phip = phi_val( d1p, d2p ) 803 phim = phi_val( d1m, d2m ) 804 phi = (phip + phim) / 2 805 dphidd = (phip - phim) / (2*dd) 806! phi0 = phi - dphidd*dsoft/2 807 phi0 = phi0_soft 808 phi2 = (4*(phi-phi0) - dphidd*dsoft) / (2*dsoft**2) 809 phi4 = (2*(phi0-phi) + dphidd*dsoft) / (2*dsoft**4) 810 phi_soft = phi0 + phi2*d**2 + phi4*d**4 811#ifdef DEBUG_XC 812! print'(a,5f8.3)', 'phi_soft: d,delta,phi0,phi2,phi4=', & 813! d, (d1-d2)/(d1+d2), phi0, phi2, phi4 814#endif /* DEBUG_XC */ 815 816 end if ! (d<=0) 817 818end function phi_soft 819 820! ----------------------------------------------------------------------------- 821 822real(dp) function phi_val( d1, d2 ) 823 824! Finds kernel phi by direct integration of Eq.(14) of 825! Dion et al, PRL 92, 246401 (2004) 826 827 real(dp),intent(in) :: d1, d2 828 829 integer :: ia, ib, n, nmesh 830 real(dp) :: a, acut, b, da1, dan, deff, dmax, dmin, gamma, & 831 pi, t, w 832 real(dp),allocatable :: amesh(:), c(:), dphida(:), dphidb(:), & 833 s(:), nu1(:), nu2(:) 834 835! Find integration mesh 836 dmax = max( abs(d1), abs(d2) ) 837 dmin = min( abs(d1), abs(d2) ) 838! deff = dmax 839 deff = sqrt(d1**2+d2**2) 840 acut = max( acutmin, acutbyd * deff ) 841 da1 = min( damax, dabyd * deff ) 842 da1 = max( da1, damin ) 843 dan = damax 844 n = get_n( 0._dp, acut, da1, dan ) 845 allocate( amesh(n), c(n), dphida(n), dphidb(n), s(n), nu1(n), nu2(n) ) 846 call set_mesh( n, xmax=acut, dxndx1=dan/da1 ) 847 call get_mesh( n, nmesh, amesh ) 848 849#ifdef DEBUG_XC 850! print'(a,i6,/,(10f8.3))', 'phi_val: size, amesh =', n, amesh 851#endif /* DEBUG_XC */ 852 853! Find cos(a), sin(a), nu1(a), and nu2(a) 854 pi = acos(-1._dp) 855 gamma = 4*pi/9 856 do ia = 1,n 857 a = amesh(ia) 858 c(ia) = cos(a) 859 s(ia) = sin(a) 860 if (a==0._dp) then 861 nu1(ia) = d1**2 / 2 / gamma 862 nu2(ia) = d2**2 / 2 / gamma 863 else ! (a/=0) 864 if (d1==0._dp) then 865 nu1(ia) = a*a / 2 866 else 867 nu1(ia) = a*a / 2 / (1-exp(-gamma*(a/d1)**2)) 868 end if 869 if (d2==0._dp) then 870 nu2(ia) = a*a / 2 871 else 872 nu2(ia) = a*a / 2 / (1-exp(-gamma*(a/d2)**2)) 873 end if 874 end if ! (a==0) 875 end do 876 877! Make integral on variable a 878 dphida(1) = 0 879 do ia = 2,n 880 a = amesh(ia) 881 882 ! Make integral on variable b 883 dphidb(1) = 0 884 do ib = 2,n 885 b = amesh(ib) 886 887 w = 2*( (3-a*a)*b*c(ib)*s(ia) + (3-b*b)*a*c(ia)*s(ib) & 888 + (a*a+b*b-3)*s(ia)*s(ib) - 3*a*b*c(ia)*c(ib) )/(a*b)**3 889 890 t = 0.5_dp * ( 1/(nu1(ia)+nu1(ib)) + 1/(nu2(ia)+nu2(ib)) ) & 891 * ( 1/(nu1(ia)+nu2(ia))/(nu1(ib)+nu2(ib)) & 892 + 1/(nu1(ia)+nu2(ib))/(nu2(ia)+nu1(ib)) ) 893 894 dphidb(ib) = a*a * b*b * w * t 895 896 end do ! ib 897 dphida(ia) = integral( n, dphidb ) 898 899 end do ! ia 900 phi_val = 2/pi**2 * integral( n, dphida ) 901 902 deallocate( amesh, c, dphida, dphidb, s, nu1, nu2 ) 903 904#ifdef DEBUG_XC 905! print'(a,2f12.6,i8,f12.6)', 'phi_val: d1,d2, na, phi =', d1, d2, n, phi_val 906#endif /* DEBUG_XC */ 907 908end function phi_val 909 910! ----------------------------------------------------------------------------- 911 912subroutine pofq( q0, p0, dp0dq0 ) 913 914! Finds the values and derivatives, at q0, of the cubic polynomials 915! p_i(q0) such that 916! y(q0) = Sum_i p_i(q0) * y_i 917! is the cubic spline interpolation at q0 of (any) function y(q) with 918! values y_i at mesh points qmesh_i 919 920 implicit none 921 real(dp),intent(in) :: q0 922 real(dp),intent(out):: p0(mq) 923 real(dp),intent(out):: dp0dq0(mq) 924 925 integer :: iq, iq0 926 real(dp):: a, b, dq 927 logical, save :: first_call=.true. 928 real(dp),save :: p(mq,mq), d2pdq2(mq,mq) 929 930! Set up spline polynomial basis 931 if (first_call) then 932 p = 0 933 do iq = 1,mq 934 p(iq,iq) = 1 935 call spline( qmesh, p(:,iq), mq, 1.e30_dp, 1.e30_dp, d2pdq2(:,iq) ) 936! call spline( qmesh, p(:,iq), mq, 0._dp, 0._dp, d2pdq2(:,iq) ) 937 end do 938 first_call = .false. 939 end if 940 941! Find interval of qmesh in which q0 is included 942 if (q0>qmesh(mq)) then ! q0 out of range 943 p0 = 0 944 dp0dq0 = 0 945 return 946 end if 947 iq0 = iofq( q0 ) 948 949! Evaluate polynomials of spline basis 950 dq = qmesh(iq0+1) - qmesh(iq0) 951 a = (qmesh(iq0+1) - q0) / dq ! dadq0 = -1/dq 952 b = (q0 - qmesh(iq0)) / dq ! dbdq0 = +1/dq 953 do iq = 1,mq 954 p0(iq) = a*p(iq0,iq) + b*p(iq0+1,iq) & 955 + ((a**3-a)*d2pdq2(iq0,iq) + (b**3-b)*d2pdq2(iq0+1,iq)) * dq**2/6 956 dp0dq0(iq) = - (p(iq0,iq) - p(iq0+1,iq)) / dq & 957 - ((3*a**2-1)*d2pdq2(iq0,iq) - (3*b**2-1)*d2pdq2(iq0+1,iq)) * dq/6 958 end do 959 960end subroutine pofq 961 962!----------------------------------------------------------------------------- 963 964subroutine qofrho( rho, grho, q, dqdrho, dqdgrho ) 965 966! Finds the local wavevector parameter q0 defined in eqs.(11-12) of 967! M.Dion et al, PRL 92, 246401 (2004) 968 969 implicit none 970 real(dp), intent(in) :: rho ! Electron density 971 real(dp), intent(in) :: grho(3) ! Density gradient 972 real(dp), intent(out):: q ! Local wave vector parameter q0 973 real(dp), intent(out):: dqdrho ! d_q/d_rho 974 real(dp), intent(out):: dqdgrho(3) ! d_q/d_grho 975 976 character(len=*),parameter :: author = 'PW92' ! Perdew-Wang'92 for LDA 977 integer, parameter :: irel = 0 ! Non-relativistic exchange 978 integer, parameter :: nspin = 1 ! Unpolarized electron gas 979 real(dp):: decdrho, dexdrho, dkfdrho, dq0dgrho(3), dq0dgrho2, dq0drho, & 980 dqdq0, dvxdrho(nspin), dvcdrho(nspin), ex, ec, grho2, & 981 kf, pi, q0, rhos(nspin), vx(nspin), vc(nspin) 982 983! Trap exception for zero density 984 if (rho <= 1.e-15_dp) then 985 q = qcut 986 dqdrho = 0 987 dqdgrho = 0 988 return 989 end if 990 991 pi = acos(-1._dp) 992 kf = (3*pi**2 * rho)**(1._dp/3) 993 994! Find exchange and correlation energy densities 995 rhos(1) = rho 996 call ldaxc( author, irel, nspin, rhos, ex, ec, vx, vc, dvxdrho, dvcdrho ) 997 998! Find q 999 grho2 = sum(grho**2) 1000 q0 = ( 1 + ec/ex - zab/9 * grho2 / (2*kf*rho)**2 ) * kf 1001 1002! Find derivatives 1003 dkfdrho = kf / (3*rho) 1004 dexdrho = (vx(1) - ex) / rho ! Since vx = d(rho*ex)/d_rho = ex + rho*dexdrho 1005 decdrho = (vc(1) - ec) / rho 1006 dq0drho = ( decdrho/ex - ec/ex**2*dexdrho + 2 * zab/9 * grho2 / & 1007 (2*kf*rho)**3 * (2*dkfdrho*rho + 2*kf) ) * kf & 1008 + q0/kf * dkfdrho 1009 dq0dgrho2 = -(zab/9) / (2*kf*rho)**2 * kf 1010 dq0dgrho(:) = dq0dgrho2 * 2*grho(:) ! Since d(vector**2)/d_vector = 2*vector 1011 1012! Saturate q to qcut smoothly 1013 call saturate( q0, qcut, q, dqdq0 ) 1014 dqdrho = dqdq0 * dq0drho 1015 dqdgrho = dqdq0 * dq0dgrho 1016 1017end subroutine qofrho 1018 1019!----------------------------------------------------------------------------- 1020 1021subroutine saturate( x, xc, y, dydx ) 1022 1023 ! Defines a function y(x) = xc * (1 - exp(-Sum_n=1:nsat (x/xc)**n/n)) 1024 ! It is approx. equal to x for x<xc and it saturates to xc when x->infinity 1025 1026 implicit none 1027 real(dp),intent(in) :: x ! Independent variable 1028 real(dp),intent(in) :: xc ! Saturation value 1029 real(dp),intent(out):: y ! Function value 1030 real(dp),intent(out):: dydx ! Derivative dy/dx 1031 1032 integer :: n 1033 real(dp):: dpdx, p 1034 1035 if (nsat >= 100) then 1036 if (x < xc) then 1037 y = x 1038 dydx = 1 1039 else ! (x >= xc) 1040 y = xc 1041 dydx = 0 1042 end if ! (x < xc) 1043 else ! (nsat < 100) 1044! This is the straightforward polynomial evaluation 1045! p = x/xc 1046! dpdx = 1/xc 1047! do n = 2,nsat 1048! p = p + (x/xc)**n / n 1049! dpdx = dpdx + (x/xc)**(n-1) / xc 1050! end do 1051! And this is a more accurate way to evaluate it 1052 p = (x/xc)/nsat 1053 dpdx = 1._dp/xc 1054 do n = nsat-1,1,-1 1055 p = (p + 1._dp/n) * x/xc 1056 dpdx = (dpdx*x + 1) / xc 1057 end do 1058 y = xc * (1 - exp(-p)) 1059 dydx = xc * dpdx * exp(-p) 1060 end if ! (nsat >= 100) 1061 1062end subroutine saturate 1063 1064!----------------------------------------------------------------------------- 1065 1066subroutine saturate_inverse( y, xc, x, dydx ) 1067 1068! Finds the inverse of the function defined in saturate subroutine 1069 1070 implicit none 1071 real(dp),intent(in) :: y ! Independent variable 1072 real(dp),intent(in) :: xc ! Saturation value 1073 real(dp),intent(out):: x ! Inverse function value 1074 real(dp),intent(out):: dydx ! Derivative dy/dx 1075 1076 real(dp):: x1, x2, yx 1077 1078 if (y<0._dp .or. y>xc) stop 'vdw:saturate_inverse: y out of range' 1079 x1 = 0 1080 x2 = xmaxbyxc * xc 1081 do 1082 x = (x1+x2)/2 1083 call saturate( x, xc, yx, dydx ) 1084 if (abs(y-yx)<ytol) then 1085 return 1086 else if (yx < y) then 1087 x1 = x 1088 else 1089 x2 = x 1090 end if 1091 end do 1092 1093end subroutine saturate_inverse 1094 1095!----------------------------------------------------------------------------- 1096 1097subroutine set_phi_table() 1098 1099! Finds and writes in disk the interpolation table (mesh points and function 1100! values) for the kernel phi(d1,d2). If the table file already exists, it is 1101! simply read and stored in memory. 1102 1103 implicit none 1104 1105 logical :: file_found 1106 integer :: id, id1, id2, nmesh 1107 real(dp):: d, d1, d1m, d1p, d2, d2m, d2p, dd, & 1108 dphidd1(nd,nd), dphidd2(nd,nd), d2phidd1dd2(nd,nd), & 1109 phi(nd,nd), phi1, phim, phip, phimm, phimp, phipm, phipp 1110 1111! Read file with table, if present 1112 inquire( file='vdw_kernel.table', exist=file_found ) 1113 if (file_found) then 1114 open( unit=1, file='vdw_kernel.table', form='unformatted' ) 1115 read(1,end=1) nmesh 1116 if (nmesh==nd) then 1117 read(1,end=1) dmesh 1118 read(1,end=1) phi_table 1119 end if 1120 close( unit=1 ) 1121 phi_table_set = .true. 1122 return 1123 end if 11241 continue ! Come here if end-of-file found 1125 1126! Set d-mesh points 1127 call set_mesh( nd, xmax=dcut, dxndx1=ddmaxddmin ) 1128 call get_mesh( nd, nmesh, dmesh ) 1129#ifdef DEBUG_XC 1130 write(udebug,'(a,/,(10f8.4))') 'set_phi_table: dmesh =', dmesh 1131#endif /* DEBUG_XC */ 1132 1133! Find function at mesh points 1134 do id1 = 1,nd 1135 d1 = dmesh(id1) 1136 phi1 = phi_soft( d1, d1 ) 1137 phi(id1,id1) = phi1 1138 do id2 = 1,id1-1 1139 d2 = dmesh(id2) 1140 d = sqrt( d1**2 + d2**2 ) 1141 if (d < dsoft) then 1142 phi(id1,id2) = phi_soft( d1, d2 ) 1143 else 1144 if (use_dphi) then ! Use dphi for better efficiency 1145 phi(id1,id2) = phi1 + dphi( d1, d2 ) 1146 else ! Use only phi_val, to eliminate uncertainties 1147 phi(id1,id2) = phi_val( d1, d2 ) 1148 end if 1149 end if 1150 phi(id2,id1) = phi(id1,id2) 1151 end do ! id2 1152 end do ! id1 1153 1154#ifdef DEBUG_XC 1155! open( unit=44, file='phi.out' ) 1156! do id2 = 1,nd 1157! write(44,'(/,(f12.6))') phi(:,id2) 1158! end do 1159! close( unit=44 ) 1160#endif /* DEBUG_XC */ 1161 1162 if (deriv_method == 'numeric') then 1163! print*, 'set_phi_table: Using numerical derivatives' 1164 1165 ! Find derivatives at mesh points 1166 do id1 = 1,nd 1167 d1 = dmesh(id1) 1168 dd = ddbyd * d1 1169 dd = max( dd, ddmin ) 1170! d1 = max( d1, dd ) 1171 d1m = d1 - dd 1172 d1p = d1 + dd 1173 phim = phi_soft( d1m, d1m ) 1174 phip = phi_soft( d1p, d1p ) 1175 do id2 = 1,id1 1176 d2 = dmesh(id2) 1177! d2 = max( d2, dd ) 1178 d = sqrt( d1**2 + d2**2 ) 1179 d1m = d1 - dd 1180 d1p = d1 + dd 1181 d2m = d2 - dd 1182 d2p = d2 + dd 1183 1184 if (d < dsoft) then 1185 phimm = phi_soft( d1m, d2m ) 1186 phipm = phi_soft( d1p, d2m ) 1187 phipp = phi_soft( d1p, d2p ) 1188 phimp = phi_soft( d1m, d2p ) 1189 else ! (d>dsoft) 1190 if (use_dphi) then 1191 phimm = phim + dphi( d1m, d2m ) 1192 phipm = phip + dphi( d1p, d2m ) 1193 phipp = phip + dphi( d1p, d2p ) 1194 if (id1==id2) then 1195 phimp = phip + dphi( d1m, d2p ) 1196 else 1197 phimp = phim + dphi( d1m, d2p ) 1198 end if 1199 else ! (.not.use_dphi) 1200 phimm = phi_val( d1m, d2m ) 1201 phipm = phi_val( d1p, d2m ) 1202 phipp = phi_val( d1p, d2p ) 1203 phimp = phi_val( d1m, d2p ) 1204 end if ! (use_dphi) 1205 end if ! (d<dsoft) 1206 1207 dphidd1(id1,id2) = (phipp+phipm-phimp-phimm) / (4*dd) 1208 dphidd2(id1,id2) = (phipp-phipm+phimp-phimm) / (4*dd) 1209 d2phidd1dd2(id1,id2) = (phipp-phipm-phimp+phimm) / (2*dd)**2 1210 dphidd1(id2,id1) = dphidd2(id1,id2) 1211 dphidd2(id2,id1) = dphidd1(id1,id2) 1212 d2phidd1dd2(id2,id1) = d2phidd1dd2(id1,id2) 1213 1214 end do ! id2 1215 end do ! id1 1216 1217 else if (deriv_method == 'interp') then 1218! print*, 'set_phi_table: Using interpolation for derivatives' 1219 1220 ! Reset mesh, which has been changed by phi_val 1221 call set_mesh( nd, xmax=dcut, dxndx1=ddmaxddmin ) 1222 1223 ! Set interpolation method 1224 if (interp_method=='Lagrange') then 1225 call set_interpolation( 'Lagrange' ) 1226 else if (interp_method=='Spline') then 1227! call set_interpolation( 'Spline', huge(phi), huge(phi) ) 1228 call set_interpolation( 'Spline', 0._dp, 0._dp ) 1229 else 1230 stop 'set_phi_val: ERROR: Unknown interp_method' 1231 end if 1232 1233 ! Find first partial derivative d_phi/d_d1 1234 do id = 1,nd 1235 dphidd1(:,id) = derivative( nd, phi(:,id) ) 1236 dphidd2(id,:) = dphidd1(:,id) 1237 end do 1238 1239 ! Find second cross partial derivative d_phi/d_d1/d_d2 1240 do id = 1,nd 1241 d2phidd1dd2(id,:) = derivative( nd, dphidd1(id,:) ) 1242 end do 1243 1244 ! Symmetrize d_phi/d_d1/d_d2 1245 do id2 = 2,nd 1246 do id1 = 1,id2-1 1247 d2phidd1dd2(id1,id2) = (d2phidd1dd2(id1,id2) + & 1248 d2phidd1dd2(id2,id1)) / 2 1249 d2phidd1dd2(id2,id1) = d2phidd1dd2(id1,id2) 1250 end do 1251 end do 1252 1253 else 1254 stop 'set_phi_table: ERROR: Unknown deriv_method' 1255 end if ! (deriv_method) 1256 1257! Make values and derivatives strictly zero when d1=dmax or d2=dmax 1258 phi(:,nd) = 0 1259 phi(nd,:) = 0 1260 dphidd1(:,nd) = 0 1261 dphidd1(nd,:) = 0 1262 dphidd2(:,nd) = 0 1263 dphidd2(nd,:) = 0 1264 d2phidd1dd2(:,nd) = 0 1265 d2phidd1dd2(nd,:) = 0 1266 1267! Make dphi(d1,d2)/dd1=0 for d1=0 and dphi(d1,d2)/dd2=0 for d2=0 1268 dphidd1(1,:) = 0 1269 dphidd2(:,1) = 0 1270 d2phidd1dd2(:,1) = 0 1271 d2phidd1dd2(1,:) = 0 1272 1273#ifdef DEBUG_XC 1274! Print values and derivatives for debugging 1275! print'(a,/,2a10,4a15)', & 1276! 'set_phi_table:', 'd1', 'd2', 'phi', 'dphi/dd1', 'dphi/dd2', 'd2phi/dd1dd2' 1277! do id1 = 1,nd 1278! do id2 = id1,nd 1279! d1 = dmesh(id1) 1280! d2 = dmesh(id2) 1281! print'(2f10.6,4e15.6)', d1, d2, phi(id1,id2), & 1282! dphidd1(id1,id2), dphidd2(id1,id2), d2phidd1dd2(id1,id2) 1283! end do 1284! end do 1285#endif /* DEBUG_XC */ 1286 1287! Set up bicubic interpolation coefficients 1288 call bcucof( nd, nd, dmesh, dmesh, phi, dphidd1, dphidd2, d2phidd1dd2, & 1289 phi_table ) 1290 1291! Save phi_table in file 1292 open( unit=1, file='vdw_kernel.table', form='unformatted' ) 1293 write(1) nd 1294 write(1) dmesh 1295 write(1) phi_table 1296 close( unit=1 ) 1297 1298#ifdef DEBUG_XC 1299! Save phi_table also in formatted file 1300 open( unit=1, file='vdw_kernel.table.formatted', form='formatted' ) 1301 write(1,*) nd 1302 write(1,*) dmesh 1303 do id2 = 1,nd 1304 do id1 = 1,nd 1305 write(1,*) phi_table(:,:,id1,id2) 1306 end do 1307 end do 1308 close( unit=1 ) 1309#endif /* DEBUG_XC */ 1310 1311! Mark table as set 1312 phi_table_set = .true. 1313 1314end subroutine set_phi_table 1315 1316! ----------------------------------------------------------------------------- 1317 1318subroutine set_qmesh() 1319 1320! Sets mesh of q values 1321 1322 implicit none 1323 integer :: nmesh 1324 1325 call set_mesh( mq, xmax=qcut, dxndx1=dqmaxdqmin ) 1326 call get_mesh( mq, nmesh, qmesh ) 1327 qmesh_set = .true. 1328 1329#ifdef DEBUG_XC 1330 write(udebug,'(/,a,/,(10f8.4))') 'vdw:set_qmesh: qmesh =', qmesh 1331#endif /* DEBUG_XC */ 1332 1333end subroutine set_qmesh 1334 1335! ----------------------------------------------------------------------------- 1336 1337subroutine vdw_decusp( nspin, rhos, grhos, eps, dedrho, dedgrho ) 1338 1339! Finds the local energy correction due to the softening of the VdW kernel 1340! cusp, defined as DEcusp(rho,grad_rho) = 1341! (1/2) * rho**2 * Int 4*pi*r**2*dr * (phi(q1*r,q2*r) - phi_soft(q1*r,q2*r)) 1342! where q1=q2=q0(rho,grad_rho). Notice that grad_rho is included in the value 1343! of q0 at the origin but not in the change of rho(r) in the integrand. 1344 1345 implicit none 1346 integer, intent(in) :: nspin ! Number of spin components 1347 real(dp),intent(in) :: rhos(nspin) ! Electron spin density 1348 real(dp),intent(in) :: grhos(3,nspin) ! Spin density gradient 1349 real(dp),intent(out):: eps ! Energy correction per electron 1350 real(dp),intent(out):: dedrho(nspin) ! d(rho*eps)/d(rho) 1351 real(dp),intent(out):: dedgrho(3,nspin) ! d(rho*eps)/d(grad_rho) 1352 1353 logical, save:: initialized=.false. 1354 real(dp),save:: table(mq,mq) 1355 integer :: iq1, iq2, ir, is, ix, ns 1356 real(dp):: dptpdq, dpdq(mq), dqdrho, dqdgrho(3), grho(3), p(mq), & 1357 ptp, phi22, phi(mq,mq), pi, pt(mq), q, r, rho 1358 1359 ! Trap VV-version exception 1360 if (vdw_author=='VV') then 1361 eps = vv_vdw_beta() 1362 dedrho = eps 1363 dedgrho = 0 1364 return 1365 end if 1366 1367 if (.not.initialized) then 1368 1369 if (.not.phi_table_set) call set_phi_table() 1370 if (.not.kcut_set) stop 'vdw_decusp: ERROR: kcut has not been set' 1371 1372 pi = acos(-1._dp) 1373 table = 0 1374 do ir = 1,nr 1375 r = ir * dr 1376 call phiofr( r, phi ) 1377 do iq2 = 1,mq 1378 do iq1 = 1,iq2 1379 table(iq1,iq2) = table(iq1,iq2) - 2*pi*dr * r**2 * phi(iq1,iq2) 1380 end do ! iq1 1381 end do ! iq2 1382 end do 1383 1384 do iq2 = 2,mq 1385 do iq1 = 1,iq2-1 1386 table(iq2,iq1) = table(iq1,iq2) 1387 end do 1388 end do 1389 1390 initialized = .true. 1391 1392 end if ! (.not.initialized) 1393 1394 ns = min(nspin,2) 1395 rho = sum(rhos(1:ns)) 1396 do ix = 1,3 1397 grho(ix) = sum(grhos(ix,1:ns)) 1398 end do 1399 1400 call qofrho( rho, grho, q, dqdrho, dqdgrho ) 1401 call pofq( q, p, dpdq ) 1402 1403 pt = matmul(p,table) 1404 ptp = sum( pt * p ) 1405 dptpdq = 2 * sum( pt * dpdq ) 1406 eps = rho * ptp 1407 dedrho(:) = 2 * rho * ptp + rho**2 * dptpdq * dqdrho 1408 do ix = 1,3 1409 dedgrho(ix,:) = rho**2 * dptpdq * dqdgrho(ix) 1410 end do 1411 1412end subroutine vdw_decusp 1413 1414!----------------------------------------------------------------------------- 1415 1416subroutine vdw_get_qmesh( n, q ) 1417 1418! Returns size and values of q-mesh 1419 1420 implicit none 1421 integer, intent(out) :: n ! Number of qmesh points 1422 real(dp),optional,intent(out) :: q(:) ! Values of qmesh points 1423 integer:: nmax, nkf, nkg 1424 1425 ! Trap VV-version exception 1426 if (vdw_author=='VV') then 1427 call vv_vdw_get_kmesh( nkf, nkg ) 1428 n = nkf*nkg 1429 if (present(q)) & 1430 call die('vdw_get_qmesh: ERROR q-mesh not available for author=VV') 1431 return 1432 end if 1433 1434 if (.not.qmesh_set) call set_qmesh() 1435 n = mq 1436 if (present(q)) then 1437 nmax = max( mq, size(q) ) 1438 q(1:nmax) = qmesh(1:nmax) 1439 end if 1440end subroutine vdw_get_qmesh 1441 1442! ----------------------------------------------------------------------------- 1443 1444subroutine vdw_localxc( iRel, nSpin, D, GD, epsX, epsC, & 1445 dEXdD, dECdD, dEXdGD, dECdGD ) 1446 1447! Finds the (semi)local part of the correlation energy density and its 1448! derivatives, using the GGA functional apropriate for the previously-set 1449! vdW functional flavour 1450 1451 implicit none 1452 integer, intent(in) :: iRel ! Relativistic exchange? 0=no, 1=yes 1453 integer, intent(in) :: nSpin ! Number of spin components 1454 real(dp),intent(in) :: D(nSpin) ! Local electron (spin) density 1455 real(dp),intent(in) :: GD(3,nSpin) ! Gradient of electron density 1456 real(dp),intent(out):: epsX ! Local exchange energy per electron 1457 real(dp),intent(out):: epsC ! Local correlation energy per electron 1458 real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 1459 real(dp),intent(out):: dECdD(nSpin) ! dEc/dDens, Ec=Int(dens*epsC) 1460 real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 1461 real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens) 1462 1463! Internal variables and arrays 1464 real(dp):: epsAux, dEdDaux(nSpin), dEdGDaux(3,nSpin) 1465 real(dp):: dVXdD(nSpin,nSpin), dVCdD(nSpin,nSpin) 1466 1467! Initialize output 1468 epsX = 0 1469 epsC = 0 1470 dEXdD = 0 1471 dECdD = 0 1472 dEXdGD = 0 1473 dECdGD = 0 1474 1475! Call the appropriate GGA functional. 1476! Use aux arrays to avoid overwritting the wrong ones 1477 if (vdw_author=='DRSLL' .or. vdw_author=='drsll' .or. & 1478 vdw_author=='DF1' .or. vdw_author=='df1') then 1479 ! Dion et al, PRL 92, 246401 (2004) 1480 ! GGA exchange and LDA correlation (we choose PW92 for the later) 1481 call GGAxc( 'revPBE', iRel, nSpin, D, GD, epsX, epsAux, & 1482 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1483 call LDAxc( 'PW92', iRel, nSpin, D, epsAux, epsC, & 1484 dEdDaux, dECdD, dVXdD, dVCdD ) 1485 else if (vdw_author=='LMKLL' .or. vdw_author=='lmkll' .or. & 1486 vdw_author=='DF2' .or. vdw_author=='df2') then 1487 ! Lee et al, PRB 82, 081101 (2010) 1488 call GGAxc( 'PW86R', iRel, nSpin, D, GD, epsX, epsAux, & 1489 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1490 call LDAxc( 'PW92', iRel, nSpin, D, epsAux, epsC, & 1491 dEdDaux, dECdD, dVXdD, dVCdD ) 1492 else if (vdw_author=='KBM' .or. vdw_author=='kbm') then 1493 ! optB88-vdW of Klimes et al, JPCM 22, 022201 (2009) 1494 call GGAxc( 'B88KBM', iRel, nSpin, D, GD, epsX, epsAux, & 1495 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1496 call LDAxc( 'PW92', iRel, nSpin, D, epsAux, epsC, & 1497 dEdDaux, dECdD, dVXdD, dVCdD ) 1498 else if (vdw_author=='C09' .or. vdw_author=='c09') then 1499 ! C09x-vdWc of Cooper, PRB 81, 161104(R) (2010) 1500 call GGAxc( 'C09', iRel, nSpin, D, GD, epsX, epsAux, & 1501 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1502 call LDAxc( 'PW92', iRel, nSpin, D, epsAux, epsC, & 1503 dEdDaux, dECdD, dVXdD, dVCdD ) 1504 else if (vdw_author=='BH' .or. vdw_author=='bh') then 1505 ! Berland and Hyldgaard, PRB 89, 035412 (2014) 1506 call GGAxc( 'BH', iRel, nSpin, D, GD, epsX, epsAux, & 1507 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1508 call LDAxc( 'PW92', iRel, nSpin, D, epsAux, epsC, & 1509 dEdDaux, dECdD, dVXdD, dVCdD ) 1510 else if (vdw_author=='VV' .or. vdw_author=='vv') then 1511 ! Vydrov and VanVoorhis, JCP 133, 244103 (2010) 1512 ! GGA for both exchange and correlation, but with different flavours 1513 call GGAxc( 'PW86R', iRel, nSpin, D, GD, epsX, epsAux, & 1514 dEXdD, dEdDaux, dEXdGD, dEdGDaux ) 1515 call GGAxc( 'PBE', iRel, nSpin, D, GD, epsAux, epsC, & 1516 dEdDaux, dECdD, dEdGDaux, dECdGD ) 1517 else 1518 stop 'vdw_exchng ERROR: unknown author' 1519 end if 1520 1521end subroutine vdw_localxc 1522 1523!----------------------------------------------------------------------------- 1524 1525subroutine vdw_phi( k, phi, dphidk ) 1526 1527! Finds by interpolation phi(q1,q2,k) (Fourier transform of phi(q1,q2,r)) 1528! for all values of q1 and q2 in qmesh. If the interpolation table does not 1529! exist, it is calculated in the first call to vdw_phi. It requires a 1530! previous call to vdw_set_kc to set qmesh. 1531 1532 implicit none 1533 real(dp),intent(in) :: k ! Modulus of actual k vector 1534 real(dp),intent(out):: phi(:,:) ! phi(q1,q2,k) at given k 1535 ! for all q1,q2 in qmesh 1536 real(dp),intent(out):: dphidk(:,:) ! dphi(q1,q2,k)/dk at given k 1537 1538 integer :: ik, iq1, iq2 1539 real(dp):: a, a2, a3, b, b2, b3 1540 1541 ! Trap VV-version exception 1542 if (vdw_author=='VV') then 1543 call vv_vdw_phi( k, phi, dphidk ) 1544 return 1545 end if 1546 1547 if (.not.kcut_set) stop 'vdw_phi: ERROR: kcut must be previously set' 1548 1549! Check argument sizes 1550 if (size(phi,1)<mq .or. size(phi,2)<mq) & 1551 stop 'vdw_phi: ERROR: size(phi) too small' 1552 1553! Find phi values at point k 1554 if (k >= kcut) then 1555 phi(:,:) = 0 1556 else 1557 ! Expand interpolation inline since this is the hottest point in VdW 1558 ik = k/dk 1559 a = ((ik+1)*dk-k)/dk 1560 b = 1 - a 1561 a2 = (3*a**2 -1) * dk / 6 1562 b2 = (3*b**2 -1) * dk / 6 1563 a3 = (a**3 - a) * dk**2 / 6 1564 b3 = (b**3 - b) * dk**2 / 6 1565 do iq2 = 1,mq 1566 do iq1 = 1,iq2 1567! call splint( dk, phik(:,iq1,iq2), d2phidk2(:,iq1,iq2), nk+1, k, & 1568! phi(iq1,iq2), dphidk(iq1,iq2), pr ) 1569 phi(iq1,iq2) = a*phik(ik,iq1,iq2) + b*phik(ik+1,iq1,iq2) & 1570 + a3*d2phidk2(ik,iq1,iq2) + b3*d2phidk2(ik+1,iq1,iq2) 1571 dphidk(iq1,iq2) = (-phik(ik,iq1,iq2) + phik(ik+1,iq1,iq2)) / dk & 1572 - a2*d2phidk2(ik,iq1,iq2) + b2*d2phidk2(ik+1,iq1,iq2) 1573 phi(iq2,iq1) = phi(iq1,iq2) 1574 dphidk(iq2,iq1) = dphidk(iq1,iq2) 1575 end do 1576 end do 1577 end if 1578 1579end subroutine vdw_phi 1580 1581!----------------------------------------------------------------------------- 1582 1583subroutine vdw_set_author( author ) 1584 1585! Sets the functional flavour (author initials) and subsequent parameters 1586 1587 implicit none 1588 character(len=*),intent(in):: author ! Functional flavour 1589 ! ('DRSLL'|'LMKLL'|'KBM'|'C09'|'BH'|'VV') 1590 1591 if (author=='DRSLL') then 1592 ! Dion et al, PRL 92, 246401 (2004) 1593 zab = -0.8491_dp 1594 else if (author=='LMKLL') then 1595 ! Lee et al, PRB 82, 081101 (2010) 1596 zab = -1.887_dp 1597 else if (author=='KBM') then 1598 ! optB88-vdW of Klimes et al, JPCM 22, 022201 (2009) 1599 zab = -0.8491_dp 1600 else if (author=='C09') then 1601 ! Cooper, PRB 81, 161104(R) (2010) 1602 zab = -0.8491_dp 1603 else if (author=='BH') then 1604 ! Berland and Hyldgaard, PRB 89, 035412 (2014) 1605 zab = -0.8491_dp 1606 else if (author=='VV') then 1607 ! Vydrov and Van Voorhis, JCP 133, 244103 (2010) 1608 else 1609 stop 'vdw_set_author: ERROR: author not known' 1610 end if 1611 vdw_author = author 1612 1613end subroutine vdw_set_author 1614 1615!----------------------------------------------------------------------------- 1616 1617subroutine vdw_set_kcut( kc ) 1618 1619! Sets the reciprocal planewave cutoff kc of the integration grid, and finds 1620! the interpolation table to be used by vdw_phi to obtain the vdW kernel phi 1621! at the reciprocal mesh wavevectors. 1622 1623 implicit none 1624 real(dp),intent(in):: kc ! Planewave cutoff: k>kcut => phi=0 1625 1626 character(len=*),parameter:: myName = 'vdw_set_kcut ' 1627 integer :: ik, iq1, iq2, ir, nrs 1628 real(dp):: dphids, dphidk0, dphidkmax, dphidr0, dphidrmax, dqdq0, & 1629 k, kmax, phi(0:nr), phi0, phi2, phis, pi, q1, q2, r(0:nr), rs 1630 1631 ! Trap VV-version exception 1632 if (vdw_author=='VV') then 1633 call vv_vdw_set_kcut( kc ) 1634 return 1635 end if 1636 1637 if (kc == kcut) return ! Work alredy done 1638 if (.not.qmesh_set) call set_qmesh() 1639 1640 ! Allocate arrays 1641 call re_alloc( phir, 0,nr, 1,mq, 1,mq, myName//'phir' ) 1642 call re_alloc( phik, 0,nr, 1,mq, 1,mq, myName//'phik' ) 1643 call re_alloc( d2phidr2, 0,nr, 1,mq, 1,mq, myName//'d2phidr2' ) 1644 call re_alloc( d2phidk2, 0,nr, 1,mq, 1,mq, myName//'d2phidk2' ) 1645 1646 pi = acos(-1._dp) 1647 dr = rcut / nr 1648 dk = pi / rcut 1649 kmax = pi / dr 1650 nk = int(kc/dk) + 1 1651 nrs = nint(rsoft/dr) 1652 rs = nrs * dr 1653 if (nk>nr) stop 'vdw_set_kcut: ERROR: nk>nr' 1654 1655#ifdef DEBUG_XC 1656 write(udebug,'(a,5f8.3)') 'vdw_set_kcut: dcut,qcut,rcut,kcut,kmax=', & 1657 dcut, qcut, rcut, kc, kmax 1658#endif /* DEBUG_XC */ 1659 1660 ! For each pair of values q1 and q2 1661 do iq2 = 1,mq 1662 do iq1 = 1,iq2 1663! print*, 'vdw_set_kcut: iq1,iq2=', iq1, iq2 1664 1665! Saturated q values 1666! q1 = qmesh(iq1) 1667! q2 = qmesh(iq2) 1668 1669! Find original (unsaturated) q values 1670 call saturate_inverse( qmesh(iq1), qcut, q1, dqdq0 ) 1671 call saturate_inverse( qmesh(iq2), qcut, q2, dqdq0 ) 1672 1673 ! Find kernel in real space 1674 do ir = 0,nr 1675 r(ir) = ir * dr 1676 phir(ir,iq1,iq2) = phi_interp( q1*r(ir), q2*r(ir) ) 1677 end do 1678 phi(:) = phir(:,iq1,iq2) 1679 1680 ! Change kernel near origin to a parabola matching at rs 1681 if (nrs>0) then 1682 phis = phi(nrs) 1683 dphids = (phi(nrs+1) - phi(nrs-1)) / (2*dr) 1684 phi0 = phis - dphids * rs/2 1685 phi2 = dphids / (2*rs) 1686 do ir = 0,nrs 1687 phir(ir,iq1,iq2) = phi0 + phi2 * r(ir)**2 1688 end do 1689 end if ! (nrs>0) 1690 1691 ! Kill kernel smoothly at r=rcut 1692 do ir = 0,nr 1693 phir(ir,iq1,iq2) = phir(ir,iq1,iq2) * cutoff( r(ir)/rcut ) 1694 end do 1695 1696 ! Optimized filter (warning: inaccurate for very large kcut*rcut) 1697! call filter( 0, nr+1, r(:), phir(:,iq1,iq2), kc, 0 ) 1698 1699 ! Find kernel in reciprocal space 1700 call radfft( 0, nr, rcut, phir(:,iq1,iq2), phik(:,iq1,iq2) ) 1701 phik(:,iq1,iq2) = phik(:,iq1,iq2) * (2*pi)**1.5_dp 1702 1703 ! Filter out above kcut 1704 phik(nk:nr,iq1,iq2) = 0 1705 1706 ! Soft filter below kcut 1707 do ik = 1,nk 1708 k = ik * dk 1709 phik(ik,iq1,iq2) = phik(ik,iq1,iq2) * cutoff(k/kc) 1710 end do 1711 1712 ! Find filtered kernel in real space 1713 call radfft( 0, nr, kmax, phik(:,iq1,iq2), phir(:,iq1,iq2) ) 1714 phir(:,iq1,iq2) = phir(:,iq1,iq2) / (2*pi)**1.5_dp 1715 1716 ! Set up spline interpolation tables 1717 dphidr0 = 0 1718 dphidrmax = 0 1719 dphidk0 = 0 1720 dphidkmax = 0 1721 call spline( dr, phir(:,iq1,iq2), nr+1, dphidr0, dphidrmax, & 1722 d2phidr2(:,iq1,iq2) ) 1723 call spline( dk, phik(:,iq1,iq2), nk+1, dphidk0, dphidkmax, & 1724 d2phidk2(:,iq1,iq2) ) 1725 1726 ! Fill symmetric elements 1727 phir(:,iq2,iq1) = phir(:,iq1,iq2) 1728 phik(:,iq2,iq1) = phik(:,iq1,iq2) 1729 d2phidr2(:,iq2,iq1) = d2phidr2(:,iq1,iq2) 1730 d2phidk2(:,iq2,iq1) = d2phidk2(:,iq1,iq2) 1731 1732! if (.false. .and. iq1==iq2) then 1733! print*, 'vdw_set_kcut: iq1,iq2=', iq1, iq2 1734! call window( 0._dp, 5._dp, -1._dp, 4._dp, 0 ) 1735! call axes( 0._dp, 1._dp, 0._dp, 1._dp ) 1736! call plot( nr+1, r, phi, phir(:,iq1,iq2) ) 1737! call window( 0._dp, 10._dp, -0.05_dp, 0.15_dp, 0 ) 1738! call axes( 0._dp, 1._dp, 0._dp, 0.05_dp ) 1739! call plot( nr+1, r, q1*q2*r**2*phi, q1*q2*r**2*phir(:,iq1,iq2) ) 1740! call show() 1741! end if 1742 1743 end do ! iq1 1744 end do ! iq2 1745 1746! print'(a,/,(2i6,f12.6))', 'vdw_set_kcut: iq1, iq2, phir(0,iq1,iq2) =', & 1747! ((iq1,iq2,phir(0,iq1,iq2),iq1=2,iq2),iq2=2,mq) 1748 1749 kcut = kc 1750 kcut_set = .true. 1751 1752end subroutine vdw_set_kcut 1753 1754!----------------------------------------------------------------------------- 1755 1756subroutine vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho ) 1757 1758! Finds the value and derivatives of theta_i(rho,grad_rho) = rho*p_i(q0), 1759! where q0(rho,grad_rho) is the local wavevector defined in eqs.(11-12) 1760! of Dion et al. p_i(q0) are the cubic polynomials such that 1761! y(q0) = Sum_i p_i(q0) * y_i 1762! is the cubic spline interpolation at q0 of (any) function y(q) with 1763! values y_i at mesh points qmesh_i 1764 1765 implicit none 1766 integer, intent(in) :: nspin ! Number of spin components 1767 real(dp),intent(in) :: rhos(nspin) ! Electron spin density 1768 real(dp),intent(in) :: grhos(3,nspin) ! Spin density gradient 1769 real(dp),intent(out):: theta(:) ! Expansion of rho*q in qmesh 1770 real(dp),intent(out):: dtdrho(:,:) ! dtheta(iq)/drhos 1771 real(dp),intent(out):: dtdgrho(:,:,:) ! dtheta(iq)/dgrhos(ix) 1772! real(dp),intent(out):: theta(mq) ! Expansion of rho*q in qmesh 1773! real(dp),intent(out):: dtdrho(mq,nspin) ! dtheta(iq)/drhos 1774! real(dp),intent(out):: dtdgrho(3,mq,nspin) ! dtheta(iq)/dgrhos(ix) 1775 1776 integer :: is, ix, ns 1777 real(dp):: rho, grho(3), dpdq(mq), dqdrho, dqdgrho(3), p(mq), q 1778 1779 ! Trap VV-version exception 1780 if (vdw_author=='VV') then 1781 call vv_vdw_theta( nspin, rhos, grhos, theta, dtdrho, dtdgrho ) 1782 return 1783 end if 1784 1785 ns = min(nspin,2) 1786 rho = sum(rhos(1:ns)) 1787 do ix = 1,3 1788 grho(ix) = sum(grhos(ix,1:ns)) 1789 end do 1790 1791 call qofrho( rho, grho, q, dqdrho, dqdgrho ) 1792 call pofq( q, p, dpdq ) 1793 1794 theta(1:mq) = rho * p(1:mq) 1795 dtdrho(:,:) = 0 1796 dtdgrho(:,:,:) = 0 1797 do is = 1,ns 1798 dtdrho(1:mq,is) = p(1:mq) + rho * dpdq(1:mq) * dqdrho 1799 do ix = 1,3 1800 dtdgrho(ix,1:mq,is) = rho * dpdq(1:mq) * dqdgrho(ix) 1801 end do 1802 end do 1803 1804end subroutine vdw_theta 1805 1806END MODULE m_vdwxc 1807