1#if defined HAVE_CONFIG_H 2#include "config.h" 3#endif 4 5MODULE m_atomXC 6 7 implicit none 8 9 PUBLIC:: atomXC ! Exchange and correlation of a spherical density 10 11 PRIVATE ! Nothing is declared public beyond this point 12 13CONTAINS 14 15!< Finds total exchange-correlation energy and potential for a 16! spherical electron density distribution. 17!``` 18! This version implements the Local (spin) Density Approximation and 19! the Generalized-Gradient-Aproximation with the 'explicit mesh 20! functional' approach of White & Bird, PRB 50, 4954 (1994). 21! Gradients are 'defined' by numerical derivatives, using 2*nn+1 mesh 22! points, where nn is a parameter defined below 23! Ref: L.C.Balbas et al, PRB 64, 165110 (2001) 24! Written by J.M.Soler using algorithms developed by 25! L.C.Balbas, J.L.Martins and J.M.Soler, Dec.1996 26! Van der Waals functional added by J.M.Soler, Jul.2008, as explained in 27! G.Roman-Perez and J.M.Soler, PRL 103, 096102 (2009) 28! ************************* INPUT *********************************** 29! INTEGER irel : Relativistic exchange? (0=>no, 1=>yes) 30! INTEGER nr : Number of radial mesh points 31! INTEGER maxr : Physical first dimension of Dens and Vxc 32! REAL*8 rmesh(nr) : Radial mesh points. Must be nr.le.maxr 33! INTEGER nSpin : nSpin=1 => unpolarized; nSpin=2 => polarized 34! REAL*8 Dens(maxr,nSpin) : Total (nSpin=1) or spin (nSpin=2) electron 35! density at mesh points 36! ************************* OUTPUT ********************************** 37! REAL*8 Ex : Total exchange energy 38! REAL*8 Ec : Total correlation energy 39! REAL*8 Dx : IntegralOf( rho * (eps_x - v_x) ) 40! REAL*8 Dc : IntegralOf( rho * (eps_c - v_c) ) 41! REAL*8 Vxc(maxr,nSpin) : (Spin) exch-corr potential 42! ************************ UNITS ************************************ 43! Distances in atomic units (Bohr). 44! Densities in atomic units (electrons/Bohr**3) 45! Energy unit depending of parameter Eunit below (currently Rydberg) 46!``` 47 48subroutine atomXC( irel, nr, maxr, rmesh, nSpin, Dens, Ex, Ec, Dx, Dc, Vxc ) 49 50! Module routines 51 use alloc, only: alloc_default ! Sets (re)allocation defaults 52 use alloc, only: de_alloc ! Deallocates arrays 53 use mesh1d, only: derivative ! Performs numerical derivative 54 use sys, only: die ! Termination routine 55 use xcmod, only: getXC ! Returns the XC functional to be used 56 use m_ggaxc, only: ggaxc ! General GGA XC routine 57 use mesh1d, only: interpolation=>interpolation_local ! Interpolation routine 58 use m_ldaxc, only: ldaxc ! General LDA XC routine 59! use m_filter,only: kcPhi ! Finds planewave cutoff of a radial func. 60 use m_radfft,only: radfft ! Radial fast Fourier transform 61 use alloc, only: re_alloc ! Reallocates arrays 62 use mesh1d, only: set_mesh ! Sets a one-dimensional mesh 63 use mesh1d, only: set_interpolation ! Sets the interpolation method 64 use m_vdwxc, only: vdw_decusp ! Cusp correction to VDW energy 65 use m_vdwxc, only: vdw_localxc ! Local LDA/GGA xc apropriate for vdW flavour 66 use m_vdwxc, only: vdw_get_qmesh ! Returns q-mesh for VDW integrals 67 use m_vdwxc, only: vdw_phi ! Returns VDW functional kernel 68 use m_vdwxc, only: vdw_set_kcut ! Fixes k-cutoff in VDW integrals 69 use m_vdwxc, only: vdw_theta ! Returns VDW theta function 70 71! Module types and variables 72 use precision, only: dp ! Double precision real type 73 use alloc, only: allocDefaults ! Derived type for allocation defaults 74 use gridxc_config, only: myNode => gridxc_myNode 75 76#ifdef DEBUG_XC 77! use m_vdwxc, only: qofrho ! Returns q(rho,grad_rho) 78 79 use debugXC, only: udebug ! Output file unit for debug info 80 use debugXC, only: setDebugOutputUnit ! Sets udebug 81#endif /* DEBUG_XC */ 82 83#ifdef HAVE_LIBXC 84 use xc_f90_types_m 85 use xc_f90_lib_m 86#endif /* HAVE_LIBXC */ 87 88 implicit none 89 90! Argument types and dimensions 91 integer, intent(in) :: irel ! Relativistic exchange? (0=>no, 1=>yes) 92 integer, intent(in) :: maxr ! First dimension of arrays dens and Vxc 93 integer, intent(in) :: nr ! Number of radial mesh points 94 integer, intent(in) :: nSpin ! Number of spin components 95 real(dp),intent(in) :: Dens(maxr,nSpin) ! (spin) Density at mesh points 96 real(dp),intent(in) :: rmesh(nr) ! Radial mesh points 97 real(dp),intent(out):: Ex ! Total exchange energy 98 real(dp),intent(out):: Ec ! Total correlation energy 99 real(dp),intent(out):: Dx ! IntegralOf(rho*(eps_x-v_x)) 100 real(dp),intent(out):: Dc ! IntegralOf(rho*(eps_c-v_c)) 101 real(dp),intent(out):: Vxc(maxr,nSpin) ! (spin) xc potential 102 103! Subroutine name 104 character(len=*),parameter :: myName = 'atomXC ' 105 character(len=*),parameter :: errHead = myName//'ERROR: ' 106 107! Fix energy unit: Eunit=1.0 => Hartrees, 108! Eunit=0.5 => Rydbergs, 109! Eunit=0.03674903 => eV 110 real(dp), parameter :: Eunit = 0.5_dp 111 112! Fix order of the numerical derivatives: used radial points = 2*nn+1 113 integer, parameter :: nn = 5 114 115! Fix number of points for radial FFTs (VDW only) 116 integer, parameter :: nk = 512 117 118! Fix the interpolation method to change to the uniform FFT mesh (VDW only) 119 character(len=*),parameter:: interp_method = 'lagrange' !('lagrange'|'spline') 120 121! Fix kin. energy leak to find the planewave cutoff of the radial density (VDW) 122! real(dp), parameter :: EtolKc = 0.003_dp ! Ry 123 124! Fix limits to the planewave cutoff of the radial density (VDW only) 125 real(dp), parameter :: kcMin = 20._dp ! Bohr^-1 126 real(dp), parameter :: kcMax = 50._dp ! Bohr^-1 127 real(dp), parameter :: Dmin = 1.e-9_dp ! Min density when estimating kc 128 real(dp), parameter :: Dcut = 1.e-9_dp ! Min density for nonzero Vxc 129 130! Fix the maximum number of functionals to be combined 131 integer, parameter :: maxFunc = 10 132 133! Max number of spin components 134 integer, parameter :: maxSpin = 4 135 136! Local variables and arrays 137 logical :: & 138 GGA, GGAfunc, VDW, VDWfunc 139 integer :: & 140 ik, in, in1, in2, iq, ir, is, ix, jn, ndSpin, nf, nq, nXCfunc 141 real(dp) :: & 142 dEcdD(maxSpin), dEcdGD(3,maxSpin), dEcuspdD(maxSpin), & 143 dExdD(maxSpin), dExdGD(3,maxSpin), dEdDaux(maxSpin), & 144 dEcuspdGD(3,maxSpin), dVcdD(maxSpin,maxSpin), dVxdD(maxSpin,maxSpin) 145 real(dp) :: & 146 dGdm(-nn:nn), d2ydx2, dk, dr, Dtot, & 147 Eaux, Ecut, epsC, epsCusp, epsNL, epsX, f1, f2, & 148 k(0:nk), kc, kmax, pi, r(0:nk), rmax, x0, xm, xp, y0, ym, yp, & 149 XCweightC(maxFunc), XCweightX(maxFunc) 150 character(len=20):: & 151 XCauth(maxFunc), XCfunc(maxFunc) 152 153 real(dp), pointer :: & 154 D(:,:)=>null(), dGDdD(:,:)=>null(), drdm(:)=>null(), dVol(:)=>null(), & 155 GD(:,:,:)=>null() 156 real(dp), pointer:: & 157 dphidk(:,:)=>null(), dtdgd(:,:,:)=>null(), dtdd(:,:)=>null(), & 158 phi(:,:)=>null(), tk(:,:)=>null(), tr(:,:)=>null(), & 159 uk(:,:)=>null(), ur(:,:)=>null() 160 type(allocDefaults):: & 161 prevAllocDefaults 162#ifdef DEBUG_XC 163! real(dp):: q, dqdrho, dqdgrho(3) 164! real(dp):: epsCtmp, dEcdDtmp(nSpin), dEcdGDtmp(3,nSpin) 165! integer:: fileUnit 166! logical:: fileOpened 167! character(len=32):: fileName 168#endif /* DEBUG_XC */ 169 170#ifdef HAVE_LIBXC 171 type(xc_f90_pointer_t), allocatable :: xc_func(:), xc_info(:) 172 logical, allocatable :: is_libxc(:) 173 integer :: xc_ispin, xc_id, iostat 174#endif 175 176#ifdef DEBUG_XC 177 call setDebugOutputUnit(myNode) ! Initialize udebug variable 178#endif /* DEBUG_XC */ 179 180! Check dimension of arrays Dens and Vxc 181 if (maxr<nr) call die(errHead//'maxr<nr') 182 183! Find number of diagonal spin components 184 ndSpin = min( nSpin, 2 ) 185 186! Get the functional(s) to be used 187 call getXC( nXCfunc, XCfunc, XCauth, XCweightX, XCweightC ) 188 189 ! Set GGA and VDW switches 190 GGA = .false. 191 VDW = .false. 192 do nf = 1,nXCfunc 193 if ( XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then 194 VDW = .true. 195 GGA = .true. 196 else if ( XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then 197 GGA = .true. 198 else if ( XCfunc(nf).ne.'LDA' .and. XCfunc(nf).ne.'lda' .and. & 199 XCfunc(nf).ne.'LSD' .and. XCfunc(nf).ne.'lsd' ) then 200 call die(errHead//'Unknown functional '//XCfunc(nf)) 201 endif 202 enddo ! nf 203 204 ! Set the possible libxc objects to avoid call overhead in grid 205#ifdef HAVE_LIBXC 206 allocate(is_libxc(nXCfunc)) 207 is_libxc(:) = .false. 208 allocate(xc_func(nXCfunc), xc_info(nXCfunc)) 209#endif 210 211 do nf = 1,nXCfunc 212 if ( XCauth(nf)(1:5) == 'LIBXC') then 213#ifdef HAVE_LIBXC 214 read(XCauth(nf)(7:10),iostat=iostat,fmt=*) xc_id 215 if (iostat /= 0) call die("Bad libxc code in " // XCauth(nf)) 216 IF (nspin == 1) THEN 217 xc_ispin = XC_UNPOLARIZED 218 ELSE 219 xc_ispin = XC_POLARIZED 220 ENDIF 221 222#include "xc_version.h" 223#if XC_MAJOR_VERSION >= 4 224 if ((irel == 1) .and. (xc_id == 1)) then 225 ! Change LDA_X to LDA_X_REL (532) 226 ! This was introduced in libXC v4 227 xc_id = 532 228 endif 229 call xc_f90_func_init(xc_func(nf), xc_info(nf), xc_id, xc_ispin) 230#else 231 call xc_f90_func_init(xc_func(nf), xc_info(nf), xc_id, xc_ispin) 232 if ((irel == 1) .and. (xc_id == 1)) then 233 ! To set the relativistic flag we need to call this 234 ! general function, and include the default 4/3 factor 235 ! that gives standard exchange (implicit in the above 236 ! initialization call) 237 call xc_f90_lda_x_set_par(xc_func(nf), 4.0_dp/3.0_dp, & 238 XC_RELATIVISTIC, 0.0_dp) 239 endif 240#endif 241 is_libxc(nf) = .true. 242#else 243 call die("Libxc not compiled in. Cannot handle "// trim(XCauth(nf))) 244#endif /* HAVE_LIBXC */ 245 endif 246 enddo ! nf 247 248! Set routine name for allocations 249! call alloc_default( old=prevAllocDefaults, routine=myName ) 250 251! Allocate temporary arrays 252 call re_alloc( D, 1,nSpin, 1,nr, myName//'D' ) 253 call re_alloc( dGDdD, -nn,nn, 1,nr, myName//'dGDdD' ) 254 call re_alloc( drdm, 1,nr, myName//'drdm' ) 255 call re_alloc( dVol, 1,nr, myName//'dVol' ) 256 call re_alloc( GD, 1,3, 1,nSpin, 1,nr, myName//'GD' ) 257 258! Find some parameters for the FFT mesh 259 pi = 4.0_dp * atan(1.0_dp) 260 rmax = rmesh(nr) 261 dr = rmax / nk 262 kmax = pi / dr 263 dk = kmax / nk 264 265! Find density and gradient of density at mesh points 266 do ir = 1,nr 267 268 ! Find interval of neighbour points to calculate derivatives 269 in1 = max( 1, ir-nn ) - ir 270 in2 = min( nr, ir+nn ) - ir 271 272 ! Find weights of numerical derivation, dGdm = d(dF(n)/dn)/dF_m, 273 ! from Lagrange interpolation formula: 274 ! F(n) = Sum_m F_m * Prod_{j/=m} (n-j)/(m-j) 275 ! dF(n)/dn = Sum_m F_m * Prod_{j/=m,j/=n} (n-j) / Prod_{j/=m} (m-j) 276 do in = in1,in2 277 if (in==0) then 278 dGdm(in) = 0 279 do jn = in1,in2 280 if (jn/=0) dGdm(in) = dGdm(in) + 1._dp / (0 - jn) 281 enddo 282 else 283 f1 = 1 284 f2 = 1 285 do jn = in1,in2 286 if (jn/=in .and. jn/=0) f1 = f1 * (0 - jn) 287 if (jn/=in) f2 = f2 * (in - jn) 288 enddo 289 dGdm(in) = f1 / f2 290 endif 291 enddo 292 293 ! Find dr(m)/dm 294 drdm(ir) = 0.0_dp 295 do in = in1,in2 296 drdm(ir) = drdm(ir) + rmesh(ir+in) * dGdm(in) 297 enddo 298 299 ! Find differential of volume. Use trapezoidal integration rule 300 dVol(ir) = 4.0_dp * pi * rmesh(ir)**2 * drdm(ir) 301 if (ir==1 .or. ir==nr) dVol(ir) = 0.5_dp*dVol(ir) 302 303 ! Find the weights for the derivative d(gradF(n))/d(F(m)), of 304 ! the gradient at point n with respect to the value at point m 305 ! dGDdD = d((dD(r)/dr)(n)/dD_m = d( (dD(n)/dn) / (dr(n)/dn) ) / dD_m 306 ! = d(dD(n)/dn)/dD_m / (dr(n)/dn) 307 if (GGA) then 308 do in = in1,in2 309 dGDdD(in,ir) = dGdm(in) / drdm(ir) 310 enddo 311 endif 312 313 ! Find density and gradient of density at this point 314 do is = 1,nSpin 315 D(is,ir) = Dens(ir,is) 316 enddo 317 if (GGA) then 318 do is = 1,nSpin 319 GD(1:3,is,ir) = 0.0_dp 320 do in = in1,in2 321 GD(3,is,ir) = GD(3,is,ir) + dGDdD(in,ir) * Dens(ir+in,is) 322 enddo 323 enddo 324 endif ! GGA 325 326 end do ! ir 327 328! Van der Waals initializations 329 if (VDW) then 330 331 ! Allocate VdW arrays 332 call vdw_get_qmesh( nq ) 333 call re_alloc( dtdd, 1,nq, 1,nSpin, myName//'dtdd' ) 334 call re_alloc( dtdgd, 1,3, 1,nq, 1,nSpin, myName//'dtdgd' ) 335 call re_alloc( dphidk, 1,nq, 1,nq, myName//'dphidk') 336 call re_alloc( phi, 1,nq, 1,nq, myName//'phi' ) 337 call re_alloc( tk, 0,nk, 1,nq, myName//'tk' ) 338 call re_alloc( tr, 1,nr, 1,nq, myName//'tr' ) 339 call re_alloc( uk, 0,nk, 1,nq, myName//'uk' ) 340 call re_alloc( ur, 1,nr, 1,nq, myName//'ur' ) 341 342 ! Find planewave cutoff of density 343! kc = kcPhi( 0, nr, rmesh(:), sum(dens(:,1:ndSpin),2), EtolKc ) 344! kc = min( kc, kcmax ) 345#ifdef DEBUG_XC 346! write(udebug,'(a,f12.6)') myName//'kc =', kc 347#endif /* DEBUG_XC */ 348 349 ! Estimate the planewave cutoff of the density 350 Ecut = 0 351 do ir = 2,nr-1 352 Dtot = sum(dens(ir,1:ndSpin)) 353 if (Dtot < Dmin) cycle 354 xm = rmesh(ir-1) 355 x0 = rmesh(ir) 356 xp = rmesh(ir+1) 357 ym = rmesh(ir-1)*sum(dens(ir-1,1:ndSpin)) 358 y0 = rmesh(ir) *sum(dens(ir ,1:ndSpin)) 359 yp = rmesh(ir+1)*sum(dens(ir+1,1:ndSpin)) 360 d2ydx2 = ( (yp-y0)/(xp-x0) - (y0-ym)/(x0-xm) ) * 2 / (xp-xm) 361 Ecut = max( Ecut, abs(d2ydx2/y0) ) 362 end do ! ir 363 kc = sqrt(Ecut) 364 kc = max( kc, kcMin ) 365 kc = min( kc, kcMax ) 366#ifdef DEBUG_XC 367! write(udebug,'(a,f12.6)') myName//'kc =', kc 368#endif /* DEBUG_XC */ 369 370 ! Set mesh cutoff to filter VdW kernel 371 call vdw_set_kcut( kc ) 372 373 ! Find expansion of theta(q(ir)) 374 do ir = 1,nr 375 call vdw_theta( nSpin, D(:,ir), GD(:,:,ir), tr(ir,1:nq), dtdd, dtdgd ) 376 end do ! ir 377 378 ! Find uniform meshes for FFTs 379 forall(ir=0:nk) r(ir) = ir*dr 380 forall(ik=0:nk) k(ik) = ik*dk 381 382 ! Find theta_iq(r) in a uniform radial mesh 383 call set_interpolation( interp_method ) 384 call set_mesh( nr, rmesh ) ! 'From' mesh of tr array 385 do iq = 1,nq 386 tk(0:nk,iq) = interpolation( nk+1, r(0:nk), nr, tr(1:nr,iq) ) 387 end do 388 389#ifdef DEBUG_XC 390! ! Write density 391! open(1,file='d.out') 392! do ir = 1,nr 393! write(1,'(3f15.9)') rmesh(ir), D(:,ir) 394! end do 395! close(1) 396! ! Write q(rho,grad_rho) 397! open(1,file='q.out') 398! do ir = 1,nr 399! call qofrho( sum(d(:,ir)), sum(gd(:,:,ir),2), q, dqdrho, dqdgrho ) 400! write(1,'(3f15.9)') rmesh(ir), q 401! end do 402! close(1) 403! ! Write theta 404! open(1,file='trmesh.out') 405! do ir = 1,nr 406! write(1,'(30f15.9)') rmesh(ir), tr(ir,1:nq) 407! end do 408! close(1) 409! open(1,file='tr.out') 410! do ir = 0,nk 411! write(1,'(30f15.9)') r(ir), tk(ir,1:nq) 412! end do 413! close(1) 414#endif /* DEBUG_XC */ 415 416 ! Fourier-transform theta_iq(r) for each iq 417 do iq = 1,nq 418 call radfft( 0, nk, rmax, tk(0:nk,iq), tk(0:nk,iq) ) 419 end do ! iq 420 421 ! Find u_iq(r) = Sum_iq' Int_dr' phi_iq_iq'(r,r')*theta_iq'(r') 422 ! by convolution in reciprocal space 423 do ik = 0,nk 424 ! Find Fourier transform of VdW kernel phi_iq_iq'(r,r') 425 call vdw_phi( k(ik), phi, dphidk ) 426 ! Find Fourier transform of u_iq(r) 427 uk(ik,1:nq) = matmul( tk(ik,1:nq), phi(1:nq,1:nq) ) 428 end do ! ik 429 430 ! Inverse Fourier transform of u_iq(r) for each iq 431 do iq = 1,nq 432 call radfft( 0, nk, kmax, uk(0:nk,iq), uk(0:nk,iq) ) 433 end do ! iq 434 435 ! Find u_iq(r) in the original radial mesh 436 call set_mesh( nk+1, xmin=0._dp, xmax=rmax ) ! 'From' mesh of uk array 437 do iq = 1,nq 438 ur(1:nr,iq) = interpolation( nr, rmesh(1:nr), nk+1, uk(0:nk,iq) ) 439 end do ! iq 440 441#ifdef DEBUG_XC 442! ! Write u 443! open(1,file='ur.out') 444! do ir = 0,nk 445! write(1,'(30f15.9)') r(ir), uk(ir,1:nq) 446! end do 447! close(1) 448! open(1,file='urmesh.out') 449! do ir = 1,nr 450! write(1,'(30f15.9)') rmesh(ir), ur(ir,1:nq) 451! end do 452! close(1) 453#endif /* DEBUG_XC */ 454 455 end if ! (VDW) 456 457! Initialize output 458 Ex = 0.0_dp 459 Ec = 0.0_dp 460 Dx = 0.0_dp 461 Dc = 0.0_dp 462 Vxc(1:nr,1:nSpin) = 0.0_dp 463 464! Loop over mesh points 465 do ir = 1,nr 466 467 ! Find interval of neighbour points to calculate derivatives 468 in1 = max( 1, ir-nn ) - ir 469 in2 = min( nr, ir+nn ) - ir 470 471 ! Loop over exchange-correlation functions 472 do nf = 1,nXCfunc 473 ! Is this a GGA or VDW? 474 if (XCfunc(nf).eq.'VDW' .or. XCfunc(nf).eq.'vdw') then 475 VDWfunc = .true. 476 GGAfunc = .true. 477 else if (XCfunc(nf).eq.'GGA' .or. XCfunc(nf).eq.'gga') then 478 VDWfunc = .false. 479 GGAfunc = .true. 480 else 481 VDWfunc = .false. 482 GGAfunc = .false. 483 endif 484 485 ! Find exchange and correlation energy densities and their 486 ! derivatives with respect to density and density gradient 487 if (VDWfunc) then 488 489 ! Local exchange-corr. part from the apropriate LDA/GGA functional 490 call vdw_localxc( irel, nSpin, D(:,ir), GD(:,:,ir), epsX, epsC, & 491 dExdD, dEcdD, dExdGD, dEcdGD ) 492 493#ifdef DEBUG_XC 494! ! Select only non local correlation energy and potential 495! epsX = 0 496! epsC = 0 497! dExdD = 0 498! dEcdD = 0 499! dExdGD = 0 500! dEcdGD = 0 501#endif /* DEBUG_XC */ 502 503 ! Local cusp correction to nonlocal VdW energy integral 504 call vdw_decusp( nSpin, D(:,ir), GD(:,:,ir), & 505 epsCusp, dEcuspdD, dEcuspdGD ) 506 507#ifdef DEBUG_XC 508! ! Select only non local correlation energy and potential 509! epsCusp = 0 510! dEcuspdD = 0 511! dEcuspdGD = 0 512#endif /* DEBUG_XC */ 513 514 ! Find expansion of theta(q(r)) for VdW 515 call vdw_theta( nSpin, D(:,ir), GD(:,:,ir), tr(ir,1:nq), dtdd, dtdgd ) 516 517 ! Add nonlocal VdW energy contribution and its derivatives 518 Dtot = sum(D(1:ndSpin,ir)) 519 epsNL = epsCusp + 0.5_dp*sum(ur(ir,1:nq)*tr(ir,1:nq))/(Dtot+tiny(Dtot)) 520 epsC = epsC + epsNL 521 do is = 1,nSpin 522 dEcdD(is) = dEcdD(is) + dEcuspdD(is) & 523 + sum(ur(ir,1:nq)*dtdd(1:nq,is)) 524 do ix = 1,3 525 dEcdGD(ix,is) = dEcdGD(ix,is) + dEcuspdGD(ix,is) & 526 + sum(ur(ir,1:nq)*dtdgd(ix,1:nq,is)) 527 end do ! ix 528 end do ! is 529 530 else if (GGAfunc) then 531 call ggaxc( XCauth(nf), irel, nSpin, D(:,ir), GD(:,:,ir), & 532 epsX, epsC, dExdD, dEcdD, dExdGD, dEcdGD & 533#ifdef HAVE_LIBXC 534 , is_libxc(nf), xc_func(nf), xc_info(nf) ) 535#else 536 ) 537#endif 538 539#ifdef DEBUG_XC 540! call ldaxc( 'PW92', irel, nSpin, D(:,ir), Eaux, epsC, & 541! dEdDaux, dEcdD, dVxdD, dVcdD ) 542! call ggaxc( XCauth(nf), irel, nSpin, D(:,ir), GD(:,:,ir), & 543! epsX, epsCtmp, dExdD, dEcdDtmp, dExdGD, dEcdGDtmp ) 544#endif /* DEBUG_XC */ 545 else 546 call ldaxc( XCauth(nf), irel, nSpin, D(:,ir), epsX, epsC, & 547 dExdD, dEcdD, dVxdD, dVcdD & 548#ifdef HAVE_LIBXC 549 , is_libxc(nf), xc_func(nf), xc_info(nf) ) 550#else 551 ) 552#endif 553 endif 554 555 ! Scale terms by weights 556 epsX = epsX*XCweightX(nf) 557 epsC = epsC*XCweightC(nf) 558 dExdD(1:nSpin) = dExdD(1:nSpin)*XCweightX(nf) 559 dEcdD(1:nSpin) = dEcdD(1:nSpin)*XCweightC(nf) 560 if (GGAfunc) then 561 dExdGD(1:3,1:nSpin) = dExdGD(1:3,1:nSpin)*XCweightX(nf) 562 dEcdGD(1:3,1:nSpin) = dEcdGD(1:3,1:nSpin)*XCweightC(nf) 563 endif 564 565 ! Add contributions to exchange-correlation energy and its 566 ! derivatives with respect to density at all points 567 do is = 1,nSpin 568 Ex = Ex + dVol(ir)*D(is,ir)*epsX 569 Ec = Ec + dVol(ir)*D(is,ir)*epsC 570 Dx = Dx + dVol(ir)*D(is,ir)*(epsX - dExdD(is)) 571 Dc = Dc + dVol(ir)*D(is,ir)*(epsC - dEcdD(is)) 572 Vxc(ir,is) = Vxc(ir,is) + dVol(ir)*(dExdD(is) + dEcdD(is)) 573 if (GGAfunc) then 574 do in = in1,in2 575 Dx= Dx - dVol(ir)*D(is,ir+in)*dExdGD(3,is)*dGDdD(in,ir) 576 Dc= Dc - dVol(ir)*D(is,ir+in)*dEcdGD(3,is)*dGDdD(in,ir) 577 Vxc(ir+in,is) = Vxc(ir+in,is) + & 578 dVol(ir)*(dExdGD(3,is) + dEcdGD(3,is))*dGDdD(in,ir) 579 enddo ! in 580 endif ! (GGAfunc) 581 enddo ! is 582 583 enddo ! nf 584 585 enddo ! ir 586 587#ifdef HAVE_LIBXC 588 do nf = 1,nXCfunc 589 if (is_libxc(nf)) then 590 call xc_f90_func_end(xc_func(nf)) 591 endif 592 enddo 593 deallocate(xc_func,xc_info,is_libxc) 594#endif 595 596! Divide by volume element to obtain the potential (per electron) 597 do is = 1,nSpin 598 ! Avoid ir=1 => r=0 => dVol=0 599 Vxc(2:nr,is) = Vxc(2:nr,is) / dVol(2:nr) 600 ! Extrapolate to the origin from first two points, requiring dVxc/di=0 601 Vxc(1,is) = (4*Vxc(2,is) - Vxc(3,is)) / 3 602 enddo ! is 603 604! Make Vxc=0 if VDWfunctl and Dens<Dcut, to avoid singularities 605 if (VDW) then 606 do ir = 1,nr 607 Dtot = sum(Dens(ir,1:ndSpin)) 608 if (Dtot<Dcut) Vxc(ir,:) = 0 609 end do 610 end if ! (VDWfunctl) 611 612! Divide by energy unit 613 Ex = Ex / Eunit 614 Ec = Ec / Eunit 615 Dx = Dx / Eunit 616 Dc = Dc / Eunit 617 Vxc(1:nr,1:nSpin) = Vxc(1:nr,1:nSpin) / Eunit 618 619#ifdef DEBUG_XC 620! ! Write potential 621! open(1,file='v.out') 622! do ir = 1,nr 623! write(1,'(3f15.9)') rmesh(ir), Vxc(ir,1:nSpin) 624! end do 625! close(1) 626#endif /* DEBUG_XC */ 627 628! Deallocate VDW-related arrays 629 if (VDW) then 630 call de_alloc( ur, myName//'ur' ) 631 call de_alloc( uk, myName//'uk' ) 632 call de_alloc( tr, myName//'tr' ) 633 call de_alloc( tk, myName//'tk' ) 634 call de_alloc( phi, myName//'phi' ) 635 call de_alloc( dphidk, myName//'dphidk' ) 636 call de_alloc( dtdgd, myName//'dtdgd' ) 637 call de_alloc( dtdd, myName//'dtdd' ) 638 endif ! (VDW) 639 640! Deallocate temporary arrays 641 call de_alloc( GD, myName//'GD' ) 642 call de_alloc( dVol, myName//'dVol' ) 643 call de_alloc( drdm, myName//'drdm' ) 644 call de_alloc( dGDdD, myName//'dGDdD' ) 645 call de_alloc( D, myName//'D' ) 646 647! Restore previous allocation defaults 648! call alloc_default( restore=prevAllocDefaults ) 649 650#ifdef DEBUG_XC 651! fileUnit = 57 652! fileName = 'atomxc.vxc' 653! inquire(file=fileName,opened=fileOpened) 654! if (.not.fileOpened) open(unit=fileUnit,file=fileName) 655! write(fileUnit,*) nr, nSpin 656! do ir = 1,nr 657! write(fileUnit,'(5e15.6)') & 658! rmesh(ir), Dens(ir,1:nSpin), Vxc(ir,1:nSpin) 659! end do 660#endif /* DEBUG_XC */ 661 662end subroutine atomXC 663 664END MODULE m_atomXC 665 666 667