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