1#if defined HAVE_CONFIG_H 2#include "config.h" 3#endif 4 5!!@LICENSE 6! 7!****************************************************************************** 8! MODULE m_ggaxc 9! Provides routines for GGA XC functional evaluation 10!****************************************************************************** 11! 12! PUBLIC procedures available from this module: 13! ggaxc, ! General subroutine for all coded GGA XC functionals 14! blypxc, ! Becke-Lee-Yang-Parr (see subroutine blypxc) 15! pbexc, ! Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996) 16! pbesolxc, ! Perdew et al, PRL, 100, 136406 (2008) 17! pw91xc, ! Perdew & Wang, JCP, 100, 1290 (1994) 18! revpbexc, ! GGA Zhang & Yang, PRL 80,890(1998) 19! rpbexc, ! Hammer, Hansen & Norskov, PRB 59, 7413 (1999) 20! am05xc, ! Mattsson & Armiento, PRB, 79, 155101 (2009) 21! wcxc, ! Wu-Cohen (see subroutine wcxc) 22! pbeJsJrLOxc ! Reparametrizations of the PBE functional by 23! pbeJsJrHEGxc ! L.S.Pedroza et al, PRB 79, 201106 (2009) and 24! pbeGcGxLOxc ! M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009) 25! pbeGcGxHEGxc ! using 4 different combinations of criteria 26! pw86x, ! Perdew & Wang, PRB 33, 8800 (1986) (exchange only) 27! pw86rx, ! pw86x refitted by Murray, Lee & Langreth, JCTC 5, 2754 (2009) 28! b88x, ! Becke, PRA 38, 3098 (1988) (exchange only) 29! b88kbmx, ! Klimes et al, JPCM 22, 022201 (2009) (exchange only) 30! c09x, ! Cooper, PRB 81, 161104 (2010) (exchange only) 31! bhx ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exchange only) 32! 33! PUBLIC parameters, types, and variables available from this module: 34! none 35! 36!****************************************************************************** 37! 38! USED module procedures: 39! use sys, only: die ! Termination routine 40! use m_ldaxc, only: exchng ! Local exchange 41! use m_ldaxc, only: pw92c ! Perdew & Wang, PRB, 45, 13244 (1992) (Corr. only) 42! 43! USED module parameters: 44! use precision, only: dp ! Real double precision type 45! 46! EXTERNAL procedures used: 47! none 48! 49!****************************************************************************** 50 51MODULE m_ggaxc 52 53 ! Used module procedures 54 use sys, only: die ! Termination routine 55 USE m_ldaxc, only: exchng ! Local exchange 56 USE m_ldaxc, only: pw92c ! Perdew & Wang, PRB, 45, 13244 (1992) correl 57 58 ! Used module parameters 59 use precision, only : dp ! Double precision real kind 60 61#ifdef HAVE_LIBXC 62 use xc_f90_types_m 63 use xc_f90_lib_m 64#endif 65 66 implicit none 67 68 private ! everything is private, unless otherwise declared 69 70 public :: ggaxc ! General subroutine for all coded GGA XC functionals 71 public :: blypxc ! Becke-Lee-Yang-Parr (see subroutine blypxc) 72 public :: pbexc ! Perdew, Burke & Ernzerhof, PRL 77, 3865 (1996) 73 public :: pbesolxc ! Perdew et al, PRL, 100, 136406 (2008) 74 public :: pw91xc ! Perdew & Wang, JCP, 100, 1290 (1994) 75 public :: revpbexc ! GGA Zhang & Yang, PRL 80,890(1998) 76 public :: rpbexc ! Hammer, Hansen & Norskov, PRB 59, 7413 (1999) 77 public :: am05xc ! Mattsson & Armiento, PRB, 79, 155101 (2009) 78 public :: wcxc ! Wu-Cohen (see subroutine wcxc) 79 public :: pbeJsJrLOxc ! Reparametrizations of the PBE functional by 80 public :: pbeJsJrHEGxc ! L.S.Pedroza et al, PRB 79, 201106 (2009) and 81 public :: pbeGcGxLOxc ! M.M.Odashima et al, JCTC 5, 798 (2009) 82 public :: pbeGcGxHEGxc ! using 4 different combinations of criteria 83 public :: pw86x ! Perdew & Wang, PRB 33, 8800 (1986) (exchange only) 84 public :: pw86rx ! pw86x refit: Murray, Lee & Langreth, JCTC 5, 2754 (2009) 85 public :: b88x ! Becke, PRA 38, 3098 (1988) (exchange only) 86 public :: b88kbmx ! Klimes et al, JPCM 22, 022201 (2009) (exchange only) 87 public :: c09x ! Cooper, PRB 81, 161104 (2010) (exchange only) 88 public :: bhx ! Berland & Hyldgaard, PRB 89, 035412 (2014) (exch. only) 89 90contains 91 92 subroutine GGAXC( AUTHOR, IREL, nSpin, D, GD, & 93 EPSX, EPSC, dEXdD, dECdD, dEXdGD, dECdGD & 94#ifdef HAVE_LIBXC 95 , is_libxc_in, xc_func, xc_info ) 96#else 97 ) 98#endif 99 100 ! Finds the exchange and correlation energies at a point, and their 101 ! derivatives with respect to density and density gradient, in the 102 ! Generalized Gradient Correction approximation. 103 ! Lengths in Bohr, energies in Hartrees 104 ! Written by L.C.Balbas and J.M.Soler, Dec'96. 105 ! Modified by V.M.Garcia-Suarez to include non-collinear spin. June 2002 106 ! Non collinear part rewritten by J.M.Soler. Sept. 2009 107 ! Interface with LibXC added by Micael Oliveira, Jul.2014 108 109 ! Input 110 character(len=*),intent(in):: AUTHOR ! GGA flavour (author initials) 111 integer, intent(in) :: IREL ! Relativistic exchange? 0=no, 1=yes 112 integer, intent(in) :: nSpin ! Number of spin components 113 real(dp),intent(in) :: D(nSpin) ! Local electron (spin) density 114 real(dp),intent(in) :: GD(3,nSpin) ! Gradient of electron density 115 116 ! Output 117 real(dp),intent(out):: EPSX ! Exchange energy per electron 118 real(dp),intent(out):: EPSC ! Correlation energy per electron 119 real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 120 real(dp),intent(out):: dECdD(nSpin) ! dEc/dDens 121 real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 122 real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens) 123#ifdef HAVE_LIBXC 124 logical, intent(in), optional :: is_libxc_in 125 type(xc_f90_pointer_t), optional :: xc_func, xc_info 126#endif 127 128 ! Internal variables and arrays 129 integer :: NS, is, ix 130 real(dp) :: DD(2), dECdDD(2), dEXdDD(2), & 131 dDDdD(2,4), dDTOTdD(4), dDPOLdD(4), & 132 dECdGDD(3,2), dEXdGDD(3,2), & 133 dGDDdD(3,2,4), dGDDdGD(2,4), & 134 dGDTOTdD(3,4), dGDPOLdD(3,4), & 135 dGDTOTdGD(4), dGDPOLdGD(4), & 136 DPOL, DTOT, GDD(3,2), GDTOT(3), GDPOL(3), & 137 VPOL 138 real(dp), parameter :: DENMIN = 1.e-12_dp 139 140 real(dp):: theta, phi, c2, s2, st, ct, cp, sp, dpolz, dpolxy 141 142 logical , parameter :: old_scheme = .true. 143 144#ifdef HAVE_LIBXC 145 logical :: is_libxc 146 real(dp) :: eps, dedn(nspin), sgm(3), dedsgm(3), dedgd(3, nspin) 147#endif 148 149 ! Handle non-collinear spin case 150 if (nSpin==4) then 151 NS = 2 ! Diagonal spin components 152 153 if ( old_scheme ) then 154 DTOT = D(1) + D(2) 155 dpolz= D(1) - D(2) 156 dpolxy= 2.0d0*sqrt(d(3)**2+d(4)**2) 157 dpol = sqrt( dpolz**2 + dpolxy**2 ) 158 if ( dpol.gt.1.0d-12 ) then 159 theta = atan2(dpolxy,dpolz) 160 else 161 theta = 0.0_dp 162 endif 163 C2 = COS(THETA/2.0_dp) 164 S2 = SIN(THETA/2.0_dp) 165 ST = SIN(THETA) 166 CT = COS(THETA) 167 PHI = ATAN2(-D(4),D(3)) 168 CP = COS(PHI) 169 SP = SIN(PHI) 170 171 DD(1) = 0.5D0 * ( DTOT + DPOL ) 172 DD(2) = 0.5D0 * ( DTOT - DPOL ) 173 DO IX = 1,3 174 GDD(IX,1) = GD(IX,1)*C2**2 + GD(IX,2)*S2**2 + & 175 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP) 176 GDD(IX,2) = GD(IX,1)*S2**2 + GD(IX,2)*C2**2 - & 177 2.d0*C2*S2*(GD(IX,3)*CP - GD(IX,4)*SP) 178 ENDDO 179 180 else 181 182 ! Find eigenvalues of density matrix Dij (diagonal densities DD, i.e. 183 ! up and down densities along the spin direction). Note convention: 184 ! D(1)=D11, D(2)=D22, D(3)=Re(D12)=Re(D21), D(4)=Im(D12)=-Im(D21) 185 DTOT = D(1) + D(2) ! DensTot (DensUp+DensDn) 186 DPOL = SQRT( (D(1)-D(2))**2 & ! DensPol (DensUp-DensDn) 187 + 4*(D(3)**2+D(4)**2) ) 188 DD(1) = ( DTOT + DPOL ) / 2 ! DensUp 189 DD(2) = ( DTOT - DPOL ) / 2 ! DensDn 190 DPOL = max( DPOL , DENMIN ) ! Avoid division by zero 191 192 ! Find gradients of up and down densities 193 GDTOT(:) = GD(:,1) + GD(:,2) ! Grad(DensTot) 194 GDPOL(:) = ( (D(1)-D(2))*(GD(:,1)-GD(:,2)) & ! Grad(DensPol) 195 + 4*(D(3)*GD(:,3)+D(4)*GD(:,4)) ) / DPOL 196 GDD(:,1) = ( GDTOT(:) + GDPOL(:) ) / 2 ! Grad(DensUp) 197 GDD(:,2) = ( GDTOT(:) - GDPOL(:) ) / 2 ! Grad(DensDn) 198 199 ! Derivatives of Dup and Ddn with respect to input density matrix 200 dDTOTdD(1:2) = 1 ! dDensTot/dD(i) 201 dDTOTdD(3:4) = 0 202 dDPOLdD(1) = +( D(1) - D(2) ) / DPOL ! dDensPol/dD(i) 203 dDPOLdD(2) = -( D(1) - D(2) ) / DPOL 204 dDPOLdD(3) = 4 * D(3) / DPOL 205 dDPOLdD(4) = 4 * D(4) / DPOL 206 dDDdD(1,:) = ( dDTOTdD(:) + dDPOLdD(:) ) / 2 ! dDensUp/dD(i) 207 dDDdD(2,:) = ( dDTOTdD(:) - dDPOLdD(:) ) / 2 ! dDensDn/dD(i) 208 209 ! Derivatives of grad(Dup) and grad(Ddn) with respect to D(i) 210 dGDTOTdD(1:3,1:4) = 0 ! dGradDensTot/dD(i) 211 dGDPOLdD(:,1) = + (GD(:,1)-GD(:,2))/DPOL & ! dGradDensPol/dD(i) 212 - GDPOL(:) * dDPOLdD(1)/DPOL 213 dGDPOLdD(:,2) = - (GD(:,1)-GD(:,2))/DPOL & 214 - GDPOL(:) * dDPOLdD(2)/DPOL 215 dGDPOLdD(:,3) = 4*GD(:,3)/DPOL & 216 - GDPOL(:) * dDPOLdD(3)/DPOL 217 dGDPOLdD(:,4) = 4*GD(:,4)/DPOL & 218 - GDPOL(:) * dDPOLdD(4)/DPOL 219 dGDDdD(:,1,:) = ( dGDTOTdD(:,:) & ! dGradDensUp/dD(i) 220 + dGDPOLdD(:,:) ) / 2 221 dGDDdD(:,2,:) = ( dGDTOTdD(:,:) & ! dGradDensDn/dD(i) 222 - dGDPOLdD(:,:) ) / 2 223 224 ! Derivatives of grad(Dup) and grad(Ddn) with respect to grad(D(i)) 225 dGDTOTdGD(1:2) = 1 ! dGradDensTot/dGradD(i) 226 dGDTOTdGD(3:4) = 0 227 dGDPOLdGD(1) = +( D(1) - D(2) ) / DPOL ! dGradDensPol/dGradD(i) 228 dGDPOLdGD(2) = -( D(1) - D(2) ) / DPOL 229 dGDPOLdGD(3) = 4 * D(3) / DPOL 230 dGDPOLdGD(4) = 4 * D(4) / DPOL 231 dGDDdGD(1,:) = ( dGDTOTdGD(:) & ! dGradDensUp/dGradD(i) 232 + dGDPOLdGD(:) ) / 2 233 dGDDdGD(2,:) = ( dGDTOTdGD(:) & ! dGradDensDn/dGradD(i) 234 - dGDPOLdGD(:) ) / 2 235 endif 236 237 else if (nSpin==1 .or. nSpin==2) then ! Normal (collinear) spin 238 NS = nSpin 239 DD(1:NS) = max( D(1:NS), 0.0_dp ) ! ag: Avoid negative densities 240 GDD(1:3,1:NS) = GD(1:3,1:NS) 241 else 242 call die('ggaxc: ERROR: invalid value of nSpin') 243 end if ! (nSpin==4) 244 245 ! Select functional to find energy density and its derivatives 246#ifdef HAVE_LIBXC 247 is_libxc = .false. 248 if (present(is_libxc_in)) then 249 if (is_libxc_in) then 250 if ((.not. present(xc_func)) .or. & 251 (.not. present(xc_info))) then 252 call die("xc_func and xc_info not present") 253 endif 254 is_libxc = .true. 255 endif 256 endif 257 258 if (is_libxc) then 259 sgm(1) = sum(GDD(1:3,1)*GDD(1:3,1)) 260 IF (nspin == 1) THEN 261 ! do nothing 262 ELSE 263 sgm(2) = sum(GDD(1:3,1)*GDD(1:3,2)) 264 sgm(3) = sum(GDD(1:3,2)*GDD(1:3,2)) 265 ENDIF 266 IF (xc_f90_info_family(xc_info) /= XC_FAMILY_GGA) THEN 267 call die('GGAXC: Functional is not a GGA') 268 ENDIF 269 270 call xc_f90_gga_exc_vxc(xc_func, 1, DD(1), sgm(1), eps, dedn(1), & 271 dedsgm(1)) 272 IF (nspin == 1) THEN 273 dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1) 274 ELSE 275 dedgd(1:3, 1) = 2.0D0*dedsgm(1)*GDD(1:3, 1) + dedsgm(2)*GDD(1:3, 2) 276 dedgd(1:3, 2) = 2.0D0*dedsgm(3)*GDD(1:3, 2) + dedsgm(2)*GDD(1:3, 1) 277 ENDIF 278 IF (xc_f90_info_kind(xc_info) == XC_CORRELATION) THEN 279 EPSC = eps 280 dECdDD(1:nspin) = dedn(1:nspin) 281 dECdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin) 282 EPSX = 0.0_dp 283 dEXdDD(1:nspin) = 0.0_dp 284 dEXdGDD(1:3, 1:nspin) = 0.0_dp 285 ELSE if (xc_f90_info_kind(xc_info) == XC_EXCHANGE) THEN 286 EPSX = eps 287 dEXdDD(1:nspin) = dedn(1:nspin) 288 dEXdGDD(1:3, 1:nspin) = dedgd(1:3, 1:nspin) 289 EPSC = 0.0_dp 290 dECdDD(1:nspin) = 0.0_dp 291 dECdGDD(1:3, 1:nspin) = 0.0_dp 292 ELSE ! combined functional, use an arbitrary 50/50 split 293 EPSX = 0.5_dp * eps 294 dEXdDD(1:nspin) = 0.5_dp * dedn(1:nspin) 295 dEXdGDD(1:3, 1:nspin) = 0.5_dp * dedgd(1:3, 1:nspin) 296 ! 297 EPSC = EPSX 298 dECdDD(1:nspin) = dEXdDD(1:nspin) 299 dECdGDD(1:3, 1:nspin) = dEXdGDD(1:3,1:nspin) 300 ENDIF 301 302 else IF (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN 303#else 304 if (AUTHOR.EQ.'PBE' .OR. AUTHOR.EQ.'pbe') THEN 305#endif 306 CALL PBEXC( IREL, NS, DD, GDD, & ! JMS 307 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 308 309 ELSE IF (AUTHOR.EQ.'RPBE'.OR.AUTHOR.EQ.'rpbe') THEN 310 CALL RPBEXC( IREL, NS, DD, GDD, & ! MVFS 311 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 312 313 ELSE IF (AUTHOR.EQ.'WC'.OR.AUTHOR.EQ.'wc') THEN 314 CALL WCXC( IREL, NS, DD, GDD, & ! MVFS 315 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 316 317 ELSE IF (AUTHOR.EQ.'REVPBE'.OR.AUTHOR.EQ.'revpbe' & 318 .OR.AUTHOR.EQ.'revPBE') THEN 319 CALL REVPBEXC( IREL, NS, DD, GDD, & ! EA 320 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 321 322 ELSE IF (AUTHOR.EQ.'BLYP'.OR.AUTHOR.EQ.'blyp'.OR. & 323 AUTHOR.EQ.'LYP'.OR.AUTHOR.EQ.'lyp') THEN 324 CALL BLYPXC( NS, DD, GDD, & ! AG 325 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 326 327 ELSEIF (AUTHOR.EQ.'PW91' .OR. AUTHOR.EQ.'pw91') THEN 328 CALL PW91XC( IREL, NS, DD, GDD, & ! AG 329 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 330 331 ELSEIF (AUTHOR.EQ.'PBESOL' .OR. AUTHOR.EQ.'pbesol' .OR. & 332 AUTHOR.EQ.'PBEsol') THEN 333 CALL PBESOLXC( IREL, NS, DD, GDD, & 334 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 335 336 ELSEIF (AUTHOR.EQ.'AM05' .OR. AUTHOR.EQ.'am05') THEN 337 CALL AM05XC( IREL, NS, DD, GDD, & 338 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 339 340 ELSEIF (AUTHOR.EQ.'PBEJsJrLO' .OR. & 341 AUTHOR.EQ.'pbejsjrlo' .OR. & 342 AUTHOR.EQ.'PBEJSJRLO') THEN 343 CALL PBEJsJrLOxc( IREL, NS, DD, GDD, & 344 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 345 346 ELSEIF (AUTHOR.EQ.'PBEJsJrHEG' .OR. & 347 AUTHOR.EQ.'pbejsjrheg' .OR. & 348 AUTHOR.EQ.'PBEJSJRHEG') THEN 349 CALL PBEJsJrHEGxc( IREL, NS, DD, GDD, & 350 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 351 352 ELSEIF (AUTHOR.EQ.'PBEGcGxLO' .OR. & 353 AUTHOR.EQ.'pbegcgxlo' .OR. & 354 AUTHOR.EQ.'PBEGCGXLO') THEN 355 CALL PBEGcGxLOxc( IREL, NS, DD, GDD, & 356 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 357 358 ELSEIF (AUTHOR.EQ.'PBEGcGxHEG' .OR. & 359 AUTHOR.EQ.'pbegcgxheg' .OR. & 360 AUTHOR.EQ.'PBEGCGXHEG') THEN 361 CALL PBEGcGxHEGxc( IREL, NS, DD, GDD, & 362 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 363 364 ELSEIF (AUTHOR.EQ.'PW86' .OR. AUTHOR.EQ.'pw86') THEN 365 CALL PW86X( IREL, NS, DD, GDD, & 366 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 367 368 ELSEIF (AUTHOR.EQ.'PW86R' .OR. AUTHOR.EQ.'pw86r') THEN 369 CALL PW86RX( IREL, NS, DD, GDD, & 370 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 371 372 ELSEIF (AUTHOR.EQ.'B88' .OR. AUTHOR.EQ.'b88') THEN 373 CALL B88X( IREL, NS, DD, GDD, & 374 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 375 376 ELSEIF (AUTHOR.EQ.'B88KBM' .OR. AUTHOR.EQ.'b88kbm') THEN 377 CALL B88KBMX( IREL, NS, DD, GDD, & 378 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 379 380 ELSEIF (AUTHOR.EQ.'C09' .OR. AUTHOR.EQ.'c09') THEN 381 CALL C09X( IREL, NS, DD, GDD, & 382 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 383 384 ELSEIF (AUTHOR.EQ.'BH' .OR. AUTHOR.EQ.'bh') THEN 385 CALL BHX( IREL, NS, DD, GDD, & 386 EPSX, EPSC, dEXdDD, dECdDD, dEXdGDD, dECdGDD ) 387 388 ELSE 389 call die('GGAXC: Unknown author ' // trim(AUTHOR)) 390 ENDIF 391 392 ! Find dE/dD(i) and dE/dGradD(i). Note convention: 393 ! DEDD(1)=dE/dD11, DEDD(2)=dE/dD22, DEDD(3)=Re(dE/dD12)=Re(dE/dD21), 394 ! DEDD(4)=Im(dE/dD12)=-Im(dE/D21) 395 if (nSpin==4) then ! Non colinear spin 396 ! dE/dD(i) = dE/dDup * dDup/dD(i) + dE/dDdn * dDdn/dD(i) 397 ! + dE/dGDup * dGDup/dD(i) + dE/dGDdn * dGDdn/dD(i) 398 ! dE/dGradD(i) = dE/dGDup * dGDup/dGD(i) + dE/dGDdn * dGDdn/dGD(i) 399 400 if ( old_scheme ) then 401 VPOL = (dExdDD(1)-dExdDD(2)) * ct 402 dExdD(1) = 0.5D0 * ( dExdDD(1) + dExdDD(2) + VPOL ) 403 dExdD(2) = 0.5D0 * ( dExdDD(1) + dExdDD(2) - VPOL ) 404 dExdD(3) = 0.5d0 * (dExdDD(1)-dExdDD(2)) * st * cp 405 dExdD(4) =-0.5d0 * (dExdDD(1)-dExdDD(2)) * st * sp 406 VPOL = (dEcdDD(1)-dEcdDD(2)) * ct 407 dEcdD(1) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) + VPOL ) 408 dEcdD(2) = 0.5D0 * ( dEcdDD(1) + dEcdDD(2) - VPOL ) 409 dEcdD(3) = 0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * cp 410 dEcdD(4) =-0.5d0 * (dEcdDD(1)-dEcdDD(2)) * st * sp 411 ! Gradient terms 412 DO IX = 1,3 413 dExdGD(IX,1) = dExdGDD(IX,1)*C2**2 + dExdGDD(IX,2)*S2**2 414 dExdGD(IX,2) = dExdGDD(IX,1)*S2**2 + dExdGDD(IX,2)*C2**2 415 dExdGD(IX,3) = 0.5D0*(dExdGDD(IX,1) - dExdGDD(IX,2))*ST*CP 416 dExdGD(IX,4) =-0.5D0*(dExdGDD(IX,1) - dExdGDD(IX,2))*ST*SP 417 dEcdGD(IX,1) = dEcdGDD(IX,1)*C2**2 + dEcdGDD(IX,2)*S2**2 418 dEcdGD(IX,2) = dEcdGDD(IX,1)*S2**2 + dEcdGDD(IX,2)*C2**2 419 dEcdGD(IX,3) = 0.5D0*(dEcdGDD(IX,1) - dEcdGDD(IX,2))*ST*CP 420 dEcdGD(IX,4) =-0.5D0*(dEcdGDD(IX,1) - dEcdGDD(IX,2))*ST*SP 421 enddo 422 423 else 424 425 do is = 1,4 426 dEXdD(is) = sum( dEXdDD(:) * dDDdD(:,is) ) & 427 + sum( dEXdGDD(:,:) * dGDDdD(:,:,is) ) 428 dECdD(is) = sum( dECdDD(:) * dDDdD(:,is) ) & 429 + sum( dECdGDD(:,:) * dGDDdD(:,:,is) ) 430 do ix = 1,3 431 dEXdGD(ix,is) = sum( dEXdGDD(ix,:) * dGDDdGD(:,is) ) 432 dECdGD(ix,is) = sum( dECdGDD(ix,:) * dGDDdGD(:,is) ) 433 end do 434 end do 435 ! Divide by two the non-diagonal derivatives. This is necessary 436 ! because DEDD(3:4) intend to be Re(dE/dD12)=Re(dE/dD21) and 437 ! Im(dE/dD12)=-Im(dE/dD21), respectively. However, both D12 and D21 438 ! depend on D(3) and D(4), and we have derived directly dE/dD(3) and 439 ! dE/dD(4). Although less trivially, the same applies to dE/dGD(3:4). 440 dEXdD(3:4) = dEXdD(3:4) / 2 441 dECdD(3:4) = dECdD(3:4) / 2 442 dEXdGD(:,3:4) = dEXdGD(:,3:4) / 2 443 dECdGD(:,3:4) = dECdGD(:,3:4) / 2 444 445 endif 446 else ! Collinear spin => just copy derivatives to output arrays 447 dEXdD(1:nSpin) = dEXdDD(1:nSpin) 448 dECdD(1:nSpin) = dECdDD(1:nSpin) 449 dEXdGD(:,1:nSpin) = dEXdGDD(:,1:nSpin) 450 dECdGD(:,1:nSpin) = dECdGDD(:,1:nSpin) 451 end if ! (nSpin==4) 452 453 END SUBROUTINE GGAXC 454 455 456 457 SUBROUTINE PBEformXC( beta, mu, kappa, iRel, nSpin, Dens, GDens, & 458 EX, EC, dEXdD, dECdD, dEXdGD, dECdGD ) 459 460 ! ********************************************************************* 461 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation 462 ! functional form, but with variable values for parameters beta, mu, and 463 ! kappa. Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996) 464 ! Written by L.C.Balbas and J.M.Soler. December 1996. 465 ! ******** INPUT ****************************************************** 466 ! REAL*8 BETA : Parameter beta of the PBE functional 467 ! REAL*8 MU : Parameter mu of the PBE functional 468 ! REAL*8 KAPPA : Parameter kappa of the PBE functional 469 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 470 ! INTEGER nspin : Number of spin polarizations (1 or 2) 471 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 472 ! spin electron density (if nspin=2) 473 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 474 ! ******** OUTPUT ***************************************************** 475 ! REAL*8 EX : Exchange energy density 476 ! REAL*8 EC : Correlation energy density 477 ! REAL*8 DEXDD(nspin) : Partial derivative 478 ! d(DensTot*Ex)/dDens(ispin), 479 ! where DensTot = Sum_ispin( Dens(ispin) ) 480 ! For a constant density, this is the 481 ! exchange potential 482 ! REAL*8 DECDD(nspin) : Partial derivative 483 ! d(DensTot*Ec)/dDens(ispin), 484 ! where DensTot = Sum_ispin( Dens(ispin) ) 485 ! For a constant density, this is the 486 ! correlation potential 487 ! REAL*8 DEXDGD(3,nspin): Partial derivative 488 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 489 ! REAL*8 DECDGD(3,nspin): Partial derivative 490 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 491 ! ********* UNITS **************************************************** 492 ! Lengths in Bohr 493 ! Densities in electrons per Bohr**3 494 ! Energies in Hartrees 495 ! Gradient vectors in cartesian coordinates 496 ! ********* ROUTINES CALLED ****************************************** 497 ! EXCHNG, PW92C 498 ! ******************************************************************** 499 500 ! Input 501 real(dp),intent(in) :: beta ! Parameter of PBE functional 502 real(dp),intent(in) :: mu ! Parameter of PBE functional 503 real(dp),intent(in) :: kappa ! Parameter of PBE functional 504 integer, intent(in) :: iRel ! Relativistic exchange? 0=no, 1=yes 505 integer, intent(in) :: nSpin ! Number of spin components 506 real(dp),intent(in) :: Dens(nSpin) ! Local electron (spin) density 507 real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density 508 509 ! Output 510 real(dp),intent(out):: EX ! Exchange energy per electron 511 real(dp),intent(out):: EC ! Correlation energy per electron 512 real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 513 real(dp),intent(out):: dECdD(nSpin) ! dEc/dDens 514 real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 515 real(dp),intent(out):: dECdGD(3,nSpin) ! dEc/dGrad(Dens) 516 517 ! Internal variables 518 INTEGER :: IS, IX 519 real(dp) & 520 A, D(2), DADD, DECUDD, DENMIN, & 521 DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, & 522 DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), & 523 DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, & 524 DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), & 525 ECUNIF, EXUNIF, & 526 F, F1, F2, F3, F4, FC, FX, FOUTHD, & 527 GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 528 H, HALF, KF, KFS, KS, PHI, PI, RS, S, & 529 T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA 530 531 ! Lower bounds of density and its gradient to avoid divisions by zero 532 PARAMETER ( DENMIN = 1.D-12 ) 533 PARAMETER ( GDMIN = 1.D-12 ) 534 535 ! Fix some numerical parameters 536 PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, & 537 THD=1.D0/3.D0, THRHLF=1.5D0, & 538 TWO=2.D0, TWOTHD=2.D0/3.D0 ) 539 540 ! Fix some more numerical constants 541 PI = 4 * ATAN(1.D0) 542 GAMMA = (1 - LOG(TWO)) / PI**2 543 544 ! Translate density and its gradient to new variables 545 IF (nspin .EQ. 1) THEN 546 D(1) = HALF * Dens(1) 547 D(2) = D(1) 548 DT = MAX( DENMIN, Dens(1) ) 549 DO IX = 1,3 550 GD(IX,1) = HALF * GDens(IX,1) 551 GD(IX,2) = GD(IX,1) 552 GDT(IX) = GDens(IX,1) 553 end DO 554 ELSE 555 D(1) = Dens(1) 556 D(2) = Dens(2) 557 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 558 DO IX = 1,3 559 GD(IX,1) = GDens(IX,1) 560 GD(IX,2) = GDens(IX,2) 561 GDT(IX) = GDens(IX,1) + GDens(IX,2) 562 end DO 563 ENDIF 564 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 565 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 566 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 567 GDMT = MAX( GDMIN, GDMT ) 568 569 ! Find local correlation energy and potential 570 CALL PW92C( 2, D, ECUNIF, VCUNIF ) 571 572 ! Find total correlation energy 573 RS = ( 3 / (4*PI*DT) )**THD 574 KF = (3 * PI**2 * DT)**THD 575 KS = SQRT( 4 * KF / PI ) 576 ZETA = ( D(1) - D(2) ) / DT 577 ZETA = MAX( -1.D0+DENMIN, ZETA ) 578 ZETA = MIN( 1.D0-DENMIN, ZETA ) 579 PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD ) 580 T = GDMT / (2 * PHI * KS * DT) 581 F1 = ECUNIF / GAMMA / PHI**3 582 F2 = EXP(-F1) 583 A = BETA / GAMMA / (F2-1) 584 F3 = T**2 + A * T**4 585 F4 = BETA/GAMMA * F3 / (1 + A*F3) 586 H = GAMMA * PHI**3 * LOG( 1 + F4 ) 587 FC = ECUNIF + H 588 589 ! Find correlation energy derivatives 590 DRSDD = - (THD * RS / DT) 591 DKFDD = THD * KF / DT 592 DKSDD = HALF * KS * DKFDD / KF 593 DZDD(1) = 1 / DT - ZETA / DT 594 DZDD(2) = - (1 / DT) - ZETA / DT 595 DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD ) 596 DO IS = 1,2 597 DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT 598 DPDD = DPDZ * DZDD(IS) 599 DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT ) 600 DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI ) 601 DF2DD = (- F2) * DF1DD 602 DADD = (- A) * DF2DD / (F2-1) 603 DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4 604 DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) ) 605 DHDD = 3 * H * DPDD / PHI 606 DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4) 607 DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD 608 609 DO IX = 1,3 610 DTDGD = (T / GDMT) * GDT(IX) / GDMT 611 DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 ) 612 DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 613 DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4) 614 DFCDGD(IX,IS) = DT * DHDGD 615 end DO 616 end DO 617 618 ! Find exchange energy and potential 619 FX = 0 620 DO IS = 1,2 621 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 622 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 623 KFS = (3 * PI**2 * DS(IS))**THD 624 S = GDMS / (2 * KFS * DS(IS)) 625 F1 = 1 + MU * S**2 / KAPPA 626 F = 1 + KAPPA - KAPPA / F1 627 ! 628 ! Note nspin=1 in call to exchng... 629 ! 630 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 631 FX = FX + DS(IS) * EXUNIF * F 632 633 DKFDD = THD * KFS / DS(IS) 634 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 635 DF1DD = 2 * (F1-1) * DSDD / S 636 DFDD = KAPPA * DF1DD / F1**2 637 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 638 639 DO IX = 1,3 640 GDS = 2 * GD(IX,IS) 641 DSDGD = (S / GDMS) * GDS / GDMS 642 DF1DGD = 2 * MU * S * DSDGD / KAPPA 643 DFDGD = KAPPA * DF1DGD / F1**2 644 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 645 end DO 646 end DO 647 FX = HALF * FX / DT 648 649 ! Set output arguments 650 EX = FX 651 EC = FC 652 DO IS = 1,nspin 653 DEXDD(IS) = DFXDD(IS) 654 DECDD(IS) = DFCDD(IS) 655 DO IX = 1,3 656 DEXDGD(IX,IS) = DFXDGD(IX,IS) 657 DECDGD(IX,IS) = DFCDGD(IX,IS) 658 end DO 659 end DO 660 END SUBROUTINE PBEformXC 661 662 663 664 SUBROUTINE PBEXC( IREL, nspin, Dens, GDens, & 665 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 666 667 ! ********************************************************************* 668 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation. 669 ! Ref: J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996) 670 ! Modified to call PBEformXC by J.M.Soler. December 2009. 671 ! ******** INPUT ****************************************************** 672 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 673 ! INTEGER nspin : Number of spin polarizations (1 or 2) 674 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 675 ! spin electron density (if nspin=2) 676 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 677 ! ******** OUTPUT ***************************************************** 678 ! REAL*8 EX : Exchange energy density 679 ! REAL*8 EC : Correlation energy density 680 ! REAL*8 DEXDD(nspin) : Partial derivative 681 ! d(DensTot*Ex)/dDens(ispin), 682 ! where DensTot = Sum_ispin( Dens(ispin) ) 683 ! For a constant density, this is the 684 ! exchange potential 685 ! REAL*8 DECDD(nspin) : Partial derivative 686 ! d(DensTot*Ec)/dDens(ispin), 687 ! where DensTot = Sum_ispin( Dens(ispin) ) 688 ! For a constant density, this is the 689 ! correlation potential 690 ! REAL*8 DEXDGD(3,nspin): Partial derivative 691 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 692 ! REAL*8 DECDGD(3,nspin): Partial derivative 693 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 694 ! ********* UNITS **************************************************** 695 ! Lengths in Bohr 696 ! Densities in electrons per Bohr**3 697 ! Energies in Hartrees 698 ! Gradient vectors in cartesian coordinates 699 ! ********* ROUTINES CALLED ****************************************** 700 ! EXCHNG, PW92C 701 ! ******************************************************************** 702 703 ! Passed arguments 704 integer, intent(in) :: IREL, NSPIN 705 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 706 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 707 DEXDD(NSPIN), DEXDGD(3,NSPIN) 708 709 ! Internal variables 710 real(dp):: BETA, KAPPA, MU, PI 711 712 ! Fix values for PBE functional parameters 713 PI = 4 * ATAN(1._dp) 714 BETA = 0.066725_dp ! From grad. exp. for correl. rs->0 715 MU = BETA * PI**2 / 3 ! From Jell. response for x+c 716 KAPPA = 0.804_dp ! From general Lieb-Oxford bound 717 718 ! Call PBE routine with appropriate values for beta, mu, and kappa. 719 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 720 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 721 722 END SUBROUTINE PBEXC 723 724 725 726 SUBROUTINE REVPBEXC( IREL, nspin, Dens, GDens,& 727 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 728 729 ! ********************************************************************* 730 ! Implements revPBE: revised Perdew-Burke-Ernzerhof GGA. 731 ! Ref: Y. Zhang & W. Yang, Phys. Rev. Lett. 80, 890 (1998). 732 ! Same interface as PBEXC. 733 ! revPBE parameters introduced by E. Artacho in January 2006 734 ! Modified to call PBEformXC by J.M.Soler. December 2009. 735 ! ******************************************************************** 736 737 ! Passed arguments 738 integer, intent(in) :: IREL, NSPIN 739 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 740 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 741 DEXDD(NSPIN), DEXDGD(3,NSPIN) 742 743 ! Internal variables 744 real(dp):: BETA, KAPPA, MU, PI 745 746 ! Fix values for PBE functional parameters 747 PI = 4 * ATAN(1._dp) 748 BETA = 0.066725_dp ! From grad. exp. for correl. rs->0 749 MU = BETA * PI**2 / 3 ! From Jell. response for x+c 750 KAPPA = 1.245_dp ! From fit of molecular energies 751 752 ! Call PBE routine with appropriate values for beta, mu, and kappa. 753 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 754 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 755 756 END SUBROUTINE REVPBEXC 757 758 759 760 SUBROUTINE PBESOLXC( IREL, NSPIN, DENS, GDENS, & 761 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 762 763 ! ********************************************************************* 764 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation. 765 ! with the revised parameters for solids (PBEsol). 766 ! Ref: J.P.Perdew et al, PRL 100, 136406 (2008) 767 ! Same interface as PBEXC. 768 ! Modified by J.D. Gale for PBEsol. May 2009. 769 ! Modified to call PBEformXC by J.M.Soler. December 2009. 770 ! ******************************************************************** 771 772 ! Passed arguments 773 integer, intent(in) :: IREL, NSPIN 774 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 775 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 776 DEXDD(NSPIN), DEXDGD(3,NSPIN) 777 778 ! Internal variables 779 real(dp):: BETA, KAPPA, MU, PI 780 781 ! Fix values for PBE functional parameters 782 PI = 4 * ATAN(1._dp) 783 BETA = 0.046_dp ! From fit of Jell. surf. energies 784 MU = 10.0_dp / 81.0_dp ! From grad. exp. for exchange. 785 KAPPA = 0.804_dp ! From general Lieb-Oxford bound 786 787 ! Call PBE routine with appropriate values for beta, mu, and kappa. 788 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 789 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 790 791 END SUBROUTINE PBESOLXC 792 793 794 795 SUBROUTINE PBEJsJrLOxc( IREL, NSPIN, DENS, GDENS, & 796 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 797 798 ! ********************************************************************* 799 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation 800 ! functional form, with revised parameters of Capelle et al: 801 ! Js refers to Jellium surface energies, that fix parameter beta 802 ! Jr refers to Jellium response, that fixes parameter mu 803 ! LO refers to Lieb-Oxford bound, that fixes parameter kappa 804 ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009) 805 ! M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009) 806 ! Same interface as PBEXC. J.M.Soler. December 2009. 807 ! ******************************************************************** 808 809 ! Passed arguments 810 integer, intent(in) :: IREL, NSPIN 811 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 812 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 813 DEXDD(NSPIN), DEXDGD(3,NSPIN) 814 815 ! Internal variables 816 real(dp):: BETA, KAPPA, MU, PI 817 818 ! Fix values for PBE functional parameters 819 PI = 4 * ATAN(1._dp) 820 BETA = 0.046_dp ! From fit of Jell. surf. energies 821 MU = BETA * PI**2 / 3 ! From Jell. response for x+c 822 KAPPA = 0.804_dp ! From general Lieb-Oxford bound 823 824 ! Call PBE routine with appropriate values for beta, mu, and kappa. 825 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 826 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 827 828 END SUBROUTINE PBEJsJrLOxc 829 830 831 832 SUBROUTINE PBEJsJrHEGxc( IREL, NSPIN, DENS, GDENS, & 833 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 834 835 ! ********************************************************************* 836 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation 837 ! functional form, with revised parameters of Capelle et al: 838 ! Js refers to Jellium surface energies, that fix parameter beta 839 ! Jr refers to Jellium response, that fixes parameter mu 840 ! HGE refers to the Lieb-Oxford bound for the low-density limit of 841 ! the homogeneous electron gas, that fixes parameter kappa 842 ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009) 843 ! M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009) 844 ! Same interface as PBEXC. J.M.Soler. December 2009. 845 ! ******************************************************************** 846 847 ! Passed arguments 848 integer, intent(in) :: IREL, NSPIN 849 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 850 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 851 DEXDD(NSPIN), DEXDGD(3,NSPIN) 852 853 ! Internal variables 854 real(dp):: BETA, KAPPA, MU, PI 855 856 ! Fix values for PBE functional parameters 857 PI = 4 * ATAN(1._dp) 858 BETA = 0.046_dp ! From fit of Jell. surf. energies 859 MU = BETA * PI**2 / 3 ! From Jell. response for x+c 860 KAPPA = 0.552_dp ! From Lieb-Oxford bound for HEG 861 862 ! Call PBE routine with appropriate values for beta, mu, and kappa. 863 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 864 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 865 866 END SUBROUTINE PBEJsJrHEGxc 867 868 869 870 SUBROUTINE PBEGcGxLOxc( IREL, NSPIN, DENS, GDENS, & 871 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 872 873 ! ********************************************************************* 874 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation 875 ! functional form, with revised parameters of Capelle et al: 876 ! Gc refers to gradient exp. for correl., that fixes parameter beta 877 ! Gx refers to grad. expansion for exchange, that fixes parameter mu 878 ! LO refers to Lieb-Oxford bound, that fixes parameter kappa 879 ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009) 880 ! M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009) 881 ! Same interface as PBEXC. J.M.Soler. December 2009. 882 ! ******************************************************************** 883 884 ! Passed arguments 885 integer, intent(in) :: IREL, NSPIN 886 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 887 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 888 DEXDD(NSPIN), DEXDGD(3,NSPIN) 889 890 ! Internal variables 891 real(dp):: BETA, KAPPA, MU, PI 892 893 ! Fix values for PBE functional parameters 894 PI = 4 * ATAN(1._dp) 895 BETA = 0.066725_dp ! From grad. exp. for correl. rs->0 896 MU = 10._dp / 81._dp ! From grad. exp. for exchange. 897 KAPPA = 0.804_dp ! From general Lieb-Oxford bound 898 899 ! Call PBE routine with appropriate values for beta, mu, and kappa. 900 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 901 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 902 903 END SUBROUTINE PBEGcGxLOxc 904 905 906 907 SUBROUTINE PBEGcGxHEGxc( IREL, NSPIN, DENS, GDENS, & 908 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 909 910 ! ********************************************************************* 911 ! Implements Perdew-Burke-Ernzerhof Generalized-Gradient-Approximation 912 ! functional form, with revised parameters of Capelle et al: 913 ! Gc refers to gradient exp. for correl., that fixes parameter beta 914 ! Gx refers to grad. expansion for exchange, that fixes parameter mu 915 ! HGE refers to the Lieb-Oxford bound for the low-density limit of 916 ! the homogeneous electron gas, that fixes parameter kappa 917 ! Refs: L.S.Pedroza et al, PRB 79, 201106 (2009) 918 ! M.M.Odashima et al, J. Chem. Theory Comp. 5, 798 (2009) 919 ! Same interface as PBEXC. J.M.Soler. December 2009. 920 ! ******************************************************************** 921 922 ! Passed arguments 923 integer, intent(in) :: IREL, NSPIN 924 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 925 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 926 DEXDD(NSPIN), DEXDGD(3,NSPIN) 927 928 ! Internal variables 929 real(dp):: BETA, KAPPA, MU, PI 930 931 ! Fix values for PBE functional parameters 932 PI = 4 * ATAN(1._dp) 933 BETA = 0.066725_dp ! From grad. exp. for correl. rs->0 934 MU = 10._dp / 81._dp ! From grad. exp. for exchange. 935 KAPPA = 0.552_dp ! From Lieb-Oxford bound for HEG 936 937 ! Call PBE routine with appropriate values for beta, mu, and kappa. 938 CALL PBEformXC( BETA, MU, KAPPA, IREL, nspin, Dens, GDens, & 939 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 940 941 END SUBROUTINE PBEGcGxHEGxc 942 943 944 945 SUBROUTINE PW91XC( IREL, nspin, Dens, GDens, & 946 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 947 948 ! ********************************************************************* 949 ! Implements Perdew-Wang91 Generalized-Gradient-Approximation. 950 ! Ref: JCP 100, 1290 (1994) 951 ! Written by J.L. Martins August 2000 952 ! ******** INPUT ****************************************************** 953 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 954 ! INTEGER nspin : Number of spin polarizations (1 or 2) 955 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 956 ! spin electron density (if nspin=2) 957 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 958 ! ******** OUTPUT ***************************************************** 959 ! REAL*8 EX : Exchange energy density 960 ! REAL*8 EC : Correlation energy density 961 ! REAL*8 DEXDD(nspin) : Partial derivative 962 ! d(DensTot*Ex)/dDens(ispin), 963 ! where DensTot = Sum_ispin( Dens(ispin) ) 964 ! For a constant density, this is the 965 ! exchange potential 966 ! REAL*8 DECDD(nspin) : Partial derivative 967 ! d(DensTot*Ec)/dDens(ispin), 968 ! where DensTot = Sum_ispin( Dens(ispin) ) 969 ! For a constant density, this is the 970 ! correlation potential 971 ! REAL*8 DEXDGD(3,nspin): Partial derivative 972 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 973 ! REAL*8 DECDGD(3,nspin): Partial derivative 974 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 975 ! ********* UNITS **************************************************** 976 ! Lengths in Bohr 977 ! Densities in electrons per Bohr**3 978 ! Energies in Hartrees 979 ! Gradient vectors in cartesian coordinates 980 ! ********* ROUTINES CALLED ****************************************** 981 ! EXCHNG, PW92C 982 ! ******************************************************************** 983 984 INTEGER IREL, nspin 985 real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin), & 986 DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin) 987 988 ! Internal variables 989 INTEGER :: IS, IX 990 real(dp) & 991 A, BETA, D(2), DADD, DECUDD, DENMIN, & 992 DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, & 993 DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), & 994 DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, & 995 DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), & 996 EC, ECUNIF, EX, EXUNIF, & 997 F, F1, F2, F3, F4, FC, FX, FOUTHD, & 998 GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 999 H, HALF, KF, KFS, KS, PHI, PI, RS, S, & 1000 T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA 1001 1002 real(dp) F5, F6, F7, F8, ASINHS 1003 real(dp) DF5DD,DF6DD,DF7DD,DF8DD 1004 real(dp) DF1DS, DF2DS, DF3DS, DFDS, DF7DGD 1005 1006 ! Lower bounds of density and its gradient to avoid divisions by zero 1007 PARAMETER ( DENMIN = 1.D-12 ) 1008 PARAMETER ( GDMIN = 1.D-12 ) 1009 1010 ! Fix some numerical parameters 1011 PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, & 1012 THD=1.D0/3.D0, THRHLF=1.5D0, & 1013 TWO=2.D0, TWOTHD=2.D0/3.D0 ) 1014 1015 ! Fix some more numerical constants 1016 PI = 4.0_dp * ATAN(1.0_dp) 1017 BETA = 15.75592_dp * 0.004235_dp 1018 GAMMA = BETA**2 / (2.0_dp * 0.09_dp) 1019 1020 ! Translate density and its gradient to new variables 1021 IF (nspin .EQ. 1) THEN 1022 D(1) = HALF * Dens(1) 1023 D(2) = D(1) 1024 DT = MAX( DENMIN, Dens(1) ) 1025 DO IX = 1,3 1026 GD(IX,1) = HALF * GDens(IX,1) 1027 GD(IX,2) = GD(IX,1) 1028 GDT(IX) = GDens(IX,1) 1029 end DO 1030 ELSE 1031 D(1) = Dens(1) 1032 D(2) = Dens(2) 1033 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 1034 DO IX = 1,3 1035 GD(IX,1) = GDens(IX,1) 1036 GD(IX,2) = GDens(IX,2) 1037 GDT(IX) = GDens(IX,1) + GDens(IX,2) 1038 end DO 1039 ENDIF 1040 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 1041 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 1042 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 1043 GDMT = MAX( GDMIN, GDMT ) 1044 1045 ! Find local correlation energy and potential 1046 CALL PW92C( 2, D, ECUNIF, VCUNIF ) 1047 1048 ! Find total correlation energy 1049 RS = ( 3 / (4*PI*DT) )**THD 1050 KF = (3 * PI**2 * DT)**THD 1051 KS = SQRT( 4 * KF / PI ) 1052 S = GDMT / (2 * KF * DT) 1053 T = GDMT / (2 * KS * DT) 1054 ZETA = ( D(1) - D(2) ) / DT 1055 ZETA = MAX( -1.D0+DENMIN, ZETA ) 1056 ZETA = MIN( 1.D0-DENMIN, ZETA ) 1057 PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD ) 1058 F1 = ECUNIF / GAMMA / PHI**3 1059 F2 = EXP(-F1) 1060 A = BETA / GAMMA / (F2-1) 1061 F3 = T**2 + A * T**4 1062 F4 = BETA/GAMMA * F3 / (1 + A*F3) 1063 F5 = 0.002568D0 + 0.023266D0*RS + 7.389D-6*RS**2 1064 F6 = 1.0D0 + 8.723D0*RS + 0.472D0*RS**2 + 0.07389D0*RS**3 1065 F7 = EXP(-100.0D0 * S**2 * PHI**4) 1066 F8 = 15.75592D0*(0.001667212D0 + F5/F6 -0.004235D0 + & 1067 3.0D0*0.001667212D0/7.0D0) 1068 H = GAMMA * PHI**3 * LOG( 1 + F4 ) + F8 * T**2 * F7 1069 FC = ECUNIF + H 1070 1071 ! Find correlation energy derivatives 1072 DRSDD = - THD * RS / DT 1073 DKFDD = THD * KF / DT 1074 DKSDD = HALF * KS * DKFDD / KF 1075 DZDD(1) = 1 / DT - ZETA / DT 1076 DZDD(2) = - 1 / DT - ZETA / DT 1077 DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD ) 1078 DO IS = 1,2 1079 DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT 1080 DPDD = DPDZ * DZDD(IS) 1081 DSDD = -S*DKFDD/KF - S/DT ! JMS: corrected, May.2014 1082 DTDD = -T*DKSDD/KS - T/DT 1083 DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI ) 1084 DF2DD = - F2 * DF1DD 1085 DADD = - A * DF2DD / (F2-1) 1086 DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4 1087 DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) ) 1088 DF5DD = (0.023266D0 + 2.0D0*7.389D-6*RS)*DRSDD 1089 DF6DD = (8.723D0 + 2.0D0*0.472D0*RS & 1090 + 3.0D0*0.07389D0*RS**2)*DRSDD 1091 DF7DD = -200.0D0 * S * PHI**4 * DSDD * F7 & 1092 -100.0D0 * S**2 * 4.0D0* PHI**3 * DPDD * F7 1093 DF8DD = 15.75592D0 * DF5DD/F6 - 15.75592D0*F5*DF6DD / F6**2 1094 DHDD = 3 * H * DPDD / PHI 1095 DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4) 1096 DHDD = DHDD + DF8DD * T**2 * F7 1097 DHDD = DHDD + F8 * 2*T*DTDD *F7 1098 DHDD = DHDD + F8 * T**2 * DF7DD 1099 1100 DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD 1101 DO IX = 1,3 1102 DTDGD = (T / GDMT) * GDT(IX) / GDMT 1103 DSDGD = (S / GDMT) * GDT(IX) / GDMT 1104 DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 ) 1105 DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 1106 DF7DGD = -200.0D0 * S * PHI**4 * DSDGD * F7 1107 DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4) 1108 DHDGD = DHDGD + F8 * 2*T*DTDGD *F7 + F8 * T**2 *DF7DGD 1109 DFCDGD(IX,IS) = DT * DHDGD 1110 end DO 1111 end DO 1112 1113 ! Find exchange energy and potential 1114 FX = 0 1115 DO IS = 1,2 1116 DS(1) = MAX( DENMIN, 2 * D(IS) ) 1117 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 1118 KFS = (3 * PI**2 * DS(1))**THD 1119 S = GDMS / (2 * KFS * DS(1)) 1120 F4 = SQRT(1.0D0 + (7.7956D0*S)**2) 1121 ASINHS = LOG(7.7956D0*S + F4) 1122 F1 = 1.0D0 + 0.19645D0 * S * ASINHS 1123 F2 = 0.2743D0 - 0.15084D0*EXP(-100.0D0*S*S) 1124 F3 = 1.0D0 / (F1 + 0.004D0 * S*S*S*S) 1125 F = (F1 + F2 * S*S ) * F3 1126 CALL EXCHNG( IREL, 1, DS, EXUNIF, VXUNIF ) 1127 FX = FX + DS(1) * EXUNIF * F 1128 1129 DKFDD = THD * KFS / DS(1) 1130 DSDD = S * ( -DKFDD/KFS - 1/DS(1) ) 1131 DF1DS = 0.19645D0 * ASINHS + & 1132 0.19645D0 * S * 7.7956D0 / F4 1133 DF2DS = 0.15084D0*200.0D0*S*EXP(-100.0D0*S*S) 1134 DF3DS = - F3*F3 * (DF1DS + 4.0D0*0.004D0 * S*S*S) 1135 DFDS = DF1DS * F3 + DF2DS * S*S * F3 + 2.0D0 * S * F2 * F3 & 1136 + (F1 + F2 * S*S ) * DF3DS 1137 DFXDD(IS) = VXUNIF(1) * F + DS(1) * EXUNIF * DFDS * DSDD 1138 1139 DO IX = 1,3 1140 GDS = 2 * GD(IX,IS) 1141 DSDGD = (S / GDMS) * GDS / GDMS 1142 DFDGD = DFDS * DSDGD 1143 DFXDGD(IX,IS) = DS(1) * EXUNIF * DFDGD 1144 end DO 1145 end DO 1146 FX = HALF * FX / DT 1147 1148 ! Set output arguments 1149 EX = FX 1150 EC = FC 1151 DO IS = 1,nspin 1152 DEXDD(IS) = DFXDD(IS) 1153 DECDD(IS) = DFCDD(IS) 1154 DO IX = 1,3 1155 DEXDGD(IX,IS) = DFXDGD(IX,IS) 1156 DECDGD(IX,IS) = DFCDGD(IX,IS) 1157 end DO 1158 end DO 1159 1160 END SUBROUTINE PW91XC 1161 1162 1163 1164 subroutine blypxc(nspin,dens,gdens,EX,EC, & 1165 dEXdd,dECdd,dEXdgd,dECdgd) 1166 ! *************************************************************** 1167 ! Implements Becke gradient exchange functional (A.D. 1168 ! Becke, Phys. Rev. A 38, 3098 (1988)) and Lee, Yang, Parr 1169 ! correlation functional (C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 1170 ! 37, 785 (1988)), as modificated by Miehlich,Savin,Stoll and Preuss, 1171 ! Chem. Phys. Lett. 157,200 (1989). See also Johnson, Gill and Pople, 1172 ! J. Chem. Phys. 98, 5612 (1993). Some errors were detected in this 1173 ! last paper, so not all of the expressions correspond exactly to those 1174 ! implemented here. 1175 ! Written by Maider Machado. July 1998. 1176 ! **************** INPUT ******************************************** 1177 ! integer nspin : Number of spin polarizations (1 or 2) 1178 ! real*8 dens(nspin) : Total electron density (if nspin=1) or 1179 ! spin electron density (if nspin=2) 1180 ! real*8 gdens(3,nspin) : Total or spin density gradient 1181 ! ******** OUTPUT ***************************************************** 1182 ! real*8 ex : Exchange energy density 1183 ! real*8 ec : Correlation energy density 1184 ! real*8 dexdd(nspin) : Partial derivative 1185 ! d(DensTot*Ex)/dDens(ispin), 1186 ! where DensTot = Sum_ispin( Dens(ispin) ) 1187 ! For a constant density, this is the 1188 ! exchange potential 1189 ! real*8 decdd(nspin) : Partial derivative 1190 ! d(DensTot*Ec)/dDens(ispin), 1191 ! where DensTot = Sum_ispin( Dens(ispin) ) 1192 ! For a constant density, this is the 1193 ! correlation potential 1194 ! real*8 dexdgd(3,nspin): Partial derivative 1195 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 1196 ! real*8 decdgd(3,nspin): Partial derivative 1197 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 1198 ! ********* UNITS **************************************************** 1199 ! Lengths in Bohr 1200 ! Densities in electrons per Bohr**3 1201 ! Energies in Hartrees 1202 ! Gradient vectors in cartesian coordinates 1203 ! ******************************************************************** 1204 1205 integer nspin 1206 real(dp) dens(nspin), gdens(3,nspin), EX, EC, & 1207 dEXdd(nspin), dECdd(nspin), dEXdgd(3,nspin), & 1208 dECdgd(3,nspin) 1209 1210 ! Internal variables 1211 integer is,ix 1212 real(dp) pi, beta, thd, tthd, thrhlf, half, fothd, & 1213 d(2),gd(3,2),dmin, ash,gdm(2),denmin,dt, & 1214 g(2),x(2),a,b,c,dd,onzthd,gdmin, & 1215 ga, gb, gc,becke,dbecgd(3,2), & 1216 dgdx(2), dgdxa, dgdxb, dgdxc,dgdxd,dbecdd(2), & 1217 den,omega, domega, delta, ddelta,cf, & 1218 gam11, gam12, gam22, LYPa, LYPb1, & 1219 LYPb2,dLYP11,dLYP12,dLYP22,LYP, & 1220 dd1g11,dd1g12,dd1g22,dd2g12,dd2g11,dd2g22, & 1221 dLYPdd(2),dg11dd(3,2),dg22dd(3,2), & 1222 dg12dd(3,2),dLYPgd(3,2) 1223 1224 ! Lower bounds of density and its gradient to avoid divisions by zero 1225 parameter ( denmin=1.0e-8_dp ) 1226 parameter ( gdmin=1.0e-8_dp ) 1227 parameter ( dmin=1.0e-5_dp ) 1228 1229 ! Fix some numerical parameters 1230 parameter ( thd = 1.d0/3.d0, tthd=2.d0/3.d0 ) 1231 parameter ( thrhlf=1.5d0, half=0.5d0, & 1232 fothd=4.d0/3.d0, onzthd=11.d0/3.d0) 1233 1234 ! Empirical parameter for Becke exchange functional (a.u.) 1235 parameter(beta= 0.0042d0) 1236 1237 ! Constants for LYP functional (a.u.) 1238 parameter(a=0.04918d0, b=0.132d0, c=0.2533d0, dd=0.349d0) 1239 1240 pi= 4*atan(1.d0) 1241 1242 ! Translate density and its gradient to new variables 1243 if (nspin .eq. 1) then 1244 d(1) = half * dens(1) 1245 d(1) = max(denmin,d(1)) 1246 d(2) = d(1) 1247 dt = max( denmin, dens(1) ) 1248 do ix = 1,3 1249 gd(ix,1) = half * gdens(ix,1) 1250 gd(ix,2) = gd(ix,1) 1251 enddo 1252 else 1253 d(1) = dens(1) 1254 d(2) = dens(2) 1255 do is=1,2 1256 d(is) = max (denmin,d(is)) 1257 enddo 1258 dt = max( denmin, dens(1)+dens(2) ) 1259 do ix = 1,3 1260 gd(ix,1) = gdens(ix,1) 1261 gd(ix,2) = gdens(ix,2) 1262 enddo 1263 endif 1264 1265 gdm(1) = sqrt( gd(1,1)**2 + gd(2,1)**2 + gd(3,1)**2 ) 1266 gdm(2) = sqrt( gd(1,2)**2 + gd(2,2)**2 + gd(3,2)**2 ) 1267 1268 do is=1,2 1269 gdm(is)= max(gdm(is),gdmin) 1270 enddo 1271 1272 ! Find Becke exchange energy 1273 ga = -thrhlf*(3.d0/4.d0/pi)**thd 1274 do is=1,2 1275 if(d(is).lt.dmin) then 1276 g(is)=ga 1277 else 1278 x(is) = gdm(is)/d(is)**fothd 1279 gb = beta*x(is)**2 1280 ash=log(x(is)+sqrt(x(is)**2+1)) 1281 gc = 1+6*beta*x(is)*ash 1282 g(is) = ga-gb/gc 1283 endif 1284 enddo 1285 1286 ! Density of energy 1287 becke=(g(1)*d(1)**fothd+g(2)*d(2)**fothd)/dt 1288 1289 ! Exchange energy derivatives 1290 do is=1,2 1291 if(d(is).lt.dmin)then 1292 dbecdd(is)=0. 1293 do ix=1,3 1294 dbecgd(ix,is)=0. 1295 enddo 1296 else 1297 dgdxa=6*beta**2*x(is)**2 1298 ash=log(x(is)+sqrt(x(is)**2+1)) 1299 dgdxb=x(is)/sqrt(x(is)**2+1)-ash 1300 dgdxc=-2*beta*x(is) 1301 dgdxd=(1+6*beta*x(is)*ash)**2 1302 dgdx(is)=(dgdxa*dgdxb+dgdxc)/dgdxd 1303 dbecdd(is)=fothd*d(is)**thd*(g(is)-x(is)*dgdx(is)) 1304 do ix=1,3 1305 dbecgd(ix,is)=d(is)**(-fothd)*dgdx(is)*gd(ix,is)/x(is) 1306 enddo 1307 endif 1308 enddo 1309 1310 ! Lee-Yang-Parr correlation energy 1311 den=1+dd*dt**(-thd) 1312 omega=dt**(-onzthd)*exp(-c*dt**(-thd))/den 1313 delta=c*dt**(-thd)+dd*dt**(-thd)/den 1314 cf=3.*(3*pi**2)**tthd/10. 1315 gam11=gdm(1)**2 1316 gam12=gd(1,1)*gd(1,2)+gd(2,1)*gd(2,2)+gd(3,1)*gd(3,2) 1317 gam22=gdm(2)**2 1318 LYPa=-4*a*d(1)*d(2)/(den*dt) 1319 LYPb1=2**onzthd*cf*a*b*omega*d(1)*d(2) 1320 LYPb2=d(1)**(8./3.)+d(2)**(8./3.) 1321 dLYP11=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.) & 1322 *d(1)/dt)-d(2)**2) 1323 dLYP12=-a*b*omega*(d(1)*d(2)/9.*(47.-7.*delta) & 1324 -fothd*dt**2) 1325 dLYP22=-a*b*omega*(d(1)*d(2)/9.*(1.-3.*delta-(delta-11.)* & 1326 d(2)/dt)-d(1)**2) 1327 1328 ! Density of energy 1329 LYP=(LYPa-LYPb1*LYPb2+dLYP11*gam11+dLYP12*gam12 & 1330 +dLYP22*gam22)/dt 1331 1332 ! Correlation energy derivatives 1333 domega=-thd*dt**(-fothd)*omega*(11.*dt**thd-c-dd/den) 1334 ddelta=thd*(dd**2*dt**(-5./3.)/den**2-delta/dt) 1335 1336 ! Second derivatives with respect to the density 1337 dd1g11=domega/omega*dLYP11-a*b*omega*(d(2)/9.* & 1338 (1.-3.*delta-2*(delta-11.)*d(1)/dt)-d(1)*d(2)/9.* & 1339 ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)) 1340 1341 dd1g12=domega/omega*dLYP12-a*b*omega*(d(2)/9.* & 1342 (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt) 1343 1344 dd1g22=domega/omega*dLYP22-a*b*omega*(1./9.*d(2) & 1345 *(1.-3.*delta-(delta-11.)*d(2)/dt)-d(1)*d(2)/9.* & 1346 ((3.+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)-2*d(1)) 1347 1348 dd2g22=domega/omega*dLYP22-a*b*omega*(d(1)/9.* & 1349 (1.-3.*delta-2*(delta-11.)*d(2)/dt)-d(1)*d(2)/9.* & 1350 ((3+d(2)/dt)*ddelta-(delta-11.)*d(2)/dt**2)) 1351 1352 dd2g12=domega/omega*dLYP12-a*b*omega*(d(1)/9.* & 1353 (47.-7.*delta)-7./9.*d(1)*d(2)*ddelta-8./3.*dt) 1354 1355 dd2g11=domega/omega*dLYP11-a*b*omega*(1./9.*d(1) & 1356 *(1.-3.*delta-(delta-11.)*d(1)/dt)-d(1)*d(2)/9.* & 1357 ((3.+d(1)/dt)*ddelta-(delta-11.)*d(1)/dt**2)-2*d(2)) 1358 1359 dLYPdd(1)=-4*a/den*d(1)*d(2)/dt* & 1360 (thd*dd*dt**(-fothd)/den & 1361 +1./d(1)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)* & 1362 (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(2)*(onzthd* & 1363 d(1)**(8./3.)+d(2)**(8./3.)))+dd1g11*gam11+ & 1364 dd1g12*gam12+dd1g22*gam22 1365 1366 dLYPdd(2)=-4*a/den*d(1)*d(2)/dt*(thd*dd*dt**(-fothd)/den & 1367 +1./d(2)-1./dt)-2**onzthd*cf*a*b*(domega*d(1)*d(2)* & 1368 (d(1)**(8./3.)+d(2)**(8./3.))+omega*d(1)*(onzthd* & 1369 d(2)**(8./3.)+d(1)**(8./3.)))+dd2g22*gam22+ & 1370 dd2g12*gam12+dd2g11*gam11 1371 1372 ! second derivatives with respect to the density gradient 1373 do is=1,2 1374 do ix=1,3 1375 dg11dd(ix,is)=2*gd(ix,is) 1376 dg22dd(ix,is)=2*gd(ix,is) 1377 enddo 1378 enddo 1379 do ix=1,3 1380 dLYPgd(ix,1)=dLYP11*dg11dd(ix,1)+dLYP12*gd(ix,2) 1381 dLYPgd(ix,2)=dLYP22*dg22dd(ix,2)+dLYP12*gd(ix,1) 1382 enddo 1383 1384 EX=becke 1385 EC=LYP 1386 do is=1,nspin 1387 dEXdd(is)=dbecdd(is) 1388 dECdd(is)=dLYPdd(is) 1389 do ix=1,3 1390 dEXdgd(ix,is)=dbecgd(ix,is) 1391 dECdgd(ix,is)=dLYPgd(ix,is) 1392 enddo 1393 enddo 1394 end subroutine blypxc 1395 1396 1397 1398 SUBROUTINE RPBEXC( IREL, nspin, Dens, GDens, & 1399 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 1400 1401 ! ********************************************************************* 1402 ! Implements Hammer's RPBE Generalized-Gradient-Approximation (GGA). 1403 ! A revision of PBE (Perdew-Burke-Ernzerhof) 1404 ! Ref: Hammer, Hansen & Norskov, PRB 59, 7413 (1999) and 1405 ! J.P.Perdew, K.Burke & M.Ernzerhof, PRL 77, 3865 (1996) 1406 ! 1407 ! Written by M.V. Fernandez-Serra. March 2004. On the PBE routine of 1408 ! L.C.Balbas and J.M.Soler. December 1996. 1409 ! ******** INPUT ****************************************************** 1410 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 1411 ! INTEGER nspin : Number of spin polarizations (1 or 2) 1412 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 1413 ! spin electron density (if nspin=2) 1414 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 1415 ! ******** OUTPUT ***************************************************** 1416 ! REAL*8 EX : Exchange energy density 1417 ! REAL*8 EC : Correlation energy density 1418 ! REAL*8 DEXDD(nspin) : Partial derivative 1419 ! d(DensTot*Ex)/dDens(ispin), 1420 ! where DensTot = Sum_ispin( Dens(ispin) ) 1421 ! For a constant density, this is the 1422 ! exchange potential 1423 ! REAL*8 DECDD(nspin) : Partial derivative 1424 ! d(DensTot*Ec)/dDens(ispin), 1425 ! where DensTot = Sum_ispin( Dens(ispin) ) 1426 ! For a constant density, this is the 1427 ! correlation potential 1428 ! REAL*8 DEXDGD(3,nspin): Partial derivative 1429 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 1430 ! REAL*8 DECDGD(3,nspin): Partial derivative 1431 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 1432 ! ********* UNITS **************************************************** 1433 ! Lengths in Bohr 1434 ! Densities in electrons per Bohr**3 1435 ! Energies in Hartrees 1436 ! Gradient vectors in cartesian coordinates 1437 ! ********* ROUTINES CALLED ****************************************** 1438 ! EXCHNG, PW92C 1439 ! ******************************************************************** 1440 1441 INTEGER IREL, nspin 1442 real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin), & 1443 DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin) 1444 1445 ! Internal variables 1446 INTEGER :: IS, IX 1447 1448 real(dp) & 1449 A, BETA, D(2), DADD, DECUDD, DENMIN, & 1450 DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, & 1451 DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), & 1452 DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, & 1453 DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), & 1454 EC, ECUNIF, EX, EXUNIF, & 1455 F, F1, F2, F3, F4, FC, FX, FOUTHD, & 1456 GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 1457 H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S, & 1458 T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA 1459 1460 ! Lower bounds of density and its gradient to avoid divisions by zero 1461 PARAMETER ( DENMIN = 1.D-12 ) 1462 PARAMETER ( GDMIN = 1.D-12 ) 1463 1464 ! Fix some numerical parameters 1465 PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, & 1466 THD=1.D0/3.D0, THRHLF=1.5D0, & 1467 TWO=2.D0, TWOTHD=2.D0/3.D0 ) 1468 1469 ! Fix some more numerical constants 1470 PI = 4 * ATAN(1.D0) 1471 BETA = 0.066725D0 1472 GAMMA = (1 - LOG(TWO)) / PI**2 1473 MU = BETA * PI**2 / 3 1474 KAPPA = 0.804D0 1475 1476 ! Translate density and its gradient to new variables 1477 IF (nspin .EQ. 1) THEN 1478 D(1) = HALF * Dens(1) 1479 D(2) = D(1) 1480 DT = MAX( DENMIN, Dens(1) ) 1481 DO IX = 1,3 1482 GD(IX,1) = HALF * GDens(IX,1) 1483 GD(IX,2) = GD(IX,1) 1484 GDT(IX) = GDens(IX,1) 1485 end DO 1486 ELSE 1487 D(1) = Dens(1) 1488 D(2) = Dens(2) 1489 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 1490 DO IX = 1,3 1491 GD(IX,1) = GDens(IX,1) 1492 GD(IX,2) = GDens(IX,2) 1493 GDT(IX) = GDens(IX,1) + GDens(IX,2) 1494 end DO 1495 ENDIF 1496 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 1497 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 1498 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 1499 GDMT = MAX( GDMIN, GDMT ) 1500 1501 ! Find local correlation energy and potential 1502 CALL PW92C( 2, D, ECUNIF, VCUNIF ) 1503 1504 ! Find total correlation energy 1505 RS = ( 3 / (4*PI*DT) )**THD 1506 KF = (3 * PI**2 * DT)**THD 1507 KS = SQRT( 4 * KF / PI ) 1508 ZETA = ( D(1) - D(2) ) / DT 1509 ZETA = MAX( -1.D0+DENMIN, ZETA ) 1510 ZETA = MIN( 1.D0-DENMIN, ZETA ) 1511 PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD ) 1512 T = GDMT / (2 * PHI * KS * DT) 1513 F1 = ECUNIF / GAMMA / PHI**3 1514 F2 = EXP(-F1) 1515 A = BETA / GAMMA / (F2-1) 1516 F3 = T**2 + A * T**4 1517 F4 = BETA/GAMMA * F3 / (1 + A*F3) 1518 H = GAMMA * PHI**3 * LOG( 1 + F4 ) 1519 FC = ECUNIF + H 1520 1521 ! Find correlation energy derivatives 1522 DRSDD = - (THD * RS / DT) 1523 DKFDD = THD * KF / DT 1524 DKSDD = HALF * KS * DKFDD / KF 1525 DZDD(1) = 1 / DT - ZETA / DT 1526 DZDD(2) = - (1 / DT) - ZETA / DT 1527 DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD ) 1528 DO IS = 1,2 1529 DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT 1530 DPDD = DPDZ * DZDD(IS) 1531 DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT ) 1532 DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI ) 1533 DF2DD = (- F2) * DF1DD 1534 DADD = (- A) * DF2DD / (F2-1) 1535 DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4 1536 DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) ) 1537 DHDD = 3 * H * DPDD / PHI 1538 DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4) 1539 DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD 1540 1541 DO IX = 1,3 1542 DTDGD = (T / GDMT) * GDT(IX) / GDMT 1543 DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 ) 1544 DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 1545 DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4) 1546 DFCDGD(IX,IS) = DT * DHDGD 1547 end DO 1548 end DO 1549 1550 ! Find exchange energy and potential 1551 FX = 0 1552 DO IS = 1,2 1553 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 1554 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 1555 KFS = (3 * PI**2 * DS(IS))**THD 1556 S = GDMS / (2 * KFS * DS(IS)) 1557 !ea Hammer's RPBE (Hammer, Hansen & Norskov PRB 59 7413 (99) 1558 !ea F1 = DEXP( - MU * S**2 / KAPPA) 1559 !ea F = 1 + KAPPA * (1 - F1) 1560 !ea Following is standard PBE 1561 !ea F1 = 1 + MU * S**2 / KAPPA 1562 !ea F = 1 + KAPPA - KAPPA / F1 1563 !ea (If revPBE Zhang & Yang, PRL 80,890(1998),change PBE's KAPPA to 1.245) 1564 F1 = DEXP( - MU * S**2 / KAPPA) 1565 F = 1 + KAPPA * (1 - F1) 1566 1567 ! Note nspin=1 in call to exchng... 1568 1569 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 1570 FX = FX + DS(IS) * EXUNIF * F 1571 1572 !MVFS The derivatives of F also need to be changed for Hammer's RPBE. 1573 !MVFS DF1DD = 2 * F1 * DSDD * ( - MU * S / KAPPA) 1574 !MVFS DF1DGD= 2 * F1 * DSDGD * ( - MU * S / KAPPA) 1575 !MVFS DFDD = -1 * KAPPA * DF1DD 1576 !MVFS DFDGD = -1 * KAPPA * DFDGD 1577 1578 DKFDD = THD * KFS / DS(IS) 1579 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 1580 ! DF1DD = 2 * (F1-1) * DSDD / S 1581 ! DFDD = KAPPA * DF1DD / F1**2 1582 DF1DD = 2* F1 * DSDD * ( - MU * S / KAPPA) 1583 DFDD = -1 * KAPPA * DF1DD 1584 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 1585 1586 DO IX = 1,3 1587 GDS = 2 * GD(IX,IS) 1588 DSDGD = (S / GDMS) * GDS / GDMS 1589 ! DF1DGD = 2 * MU * S * DSDGD / KAPPA 1590 ! DFDGD = KAPPA * DF1DGD / F1**2 1591 DF1DGD =2*F1 * DSDGD * ( - MU * S / KAPPA) 1592 DFDGD = -1 * KAPPA * DF1DGD 1593 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 1594 end DO 1595 end DO 1596 FX = HALF * FX / DT 1597 1598 ! Set output arguments 1599 EX = FX 1600 EC = FC 1601 DO IS = 1,nspin 1602 DEXDD(IS) = DFXDD(IS) 1603 DECDD(IS) = DFCDD(IS) 1604 DO IX = 1,3 1605 DEXDGD(IX,IS) = DFXDGD(IX,IS) 1606 DECDGD(IX,IS) = DFCDGD(IX,IS) 1607 end DO 1608 end DO 1609 END SUBROUTINE RPBEXC 1610 1611 1612 1613 SUBROUTINE WCXC( IREL, nspin, Dens, GDens, & 1614 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 1615 1616 ! ********************************************************************* 1617 ! Implements Wu-Cohen Generalized-Gradient-Approximation. 1618 ! Ref: Z. Wu and R. E. Cohen PRB 73, 235116 (2006) 1619 ! Written by Marivi Fernandez-Serra, with contributions by 1620 ! Julian Gale and Alberto Garcia, 1621 ! over the PBEXC subroutine of L.C.Balbas and J.M.Soler. 1622 ! September, 2006. 1623 ! ******** INPUT ****************************************************** 1624 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 1625 ! INTEGER nspin : Number of spin polarizations (1 or 2) 1626 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 1627 ! spin electron density (if nspin=2) 1628 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 1629 ! ******** OUTPUT ***************************************************** 1630 ! REAL*8 EX : Exchange energy density 1631 ! REAL*8 EC : Correlation energy density 1632 ! REAL*8 DEXDD(nspin) : Partial derivative 1633 ! d(DensTot*Ex)/dDens(ispin), 1634 ! where DensTot = Sum_ispin( Dens(ispin) ) 1635 ! For a constant density, this is the 1636 ! exchange potential 1637 ! REAL*8 DECDD(nspin) : Partial derivative 1638 ! d(DensTot*Ec)/dDens(ispin), 1639 ! where DensTot = Sum_ispin( Dens(ispin) ) 1640 ! For a constant density, this is the 1641 ! correlation potential 1642 ! REAL*8 DEXDGD(3,nspin): Partial derivative 1643 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 1644 ! REAL*8 DECDGD(3,nspin): Partial derivative 1645 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 1646 ! ********* UNITS **************************************************** 1647 ! Lengths in Bohr 1648 ! Densities in electrons per Bohr**3 1649 ! Energies in Hartrees 1650 ! Gradient vectors in cartesian coordinates 1651 ! ********* ROUTINES CALLED ****************************************** 1652 ! EXCHNG, PW92C 1653 ! ******************************************************************** 1654 1655 INTEGER IREL, nspin 1656 real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin), & 1657 DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin) 1658 1659 ! Internal variables 1660 INTEGER :: IS, IX 1661 1662 real(dp) & 1663 A, BETA, D(2), DADD, DECUDD, DENMIN, & 1664 DF1DD, DF2DD, DF3DD, DF4DD, DF1DGD, DF3DGD, DF4DGD, & 1665 DFCDD(2), DFCDGD(3,2), DFDD, DFDGD, DFXDD(2), DFXDGD(3,2), & 1666 DHDD, DHDGD, DKFDD, DKSDD, DPDD, DPDZ, DRSDD, & 1667 DS(2), DSDD, DSDGD, DT, DTDD, DTDGD, DZDD(2), & 1668 XWC, DXWCDS, CWC, & 1669 EC, ECUNIF, EX, EXUNIF, & 1670 F, F1, F2, F3, F4, FC, FX, FOUTHD, & 1671 GAMMA, GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 1672 H, HALF, KAPPA, KF, KFS, KS, MU, PHI, PI, RS, S, & 1673 TEN81, & 1674 T, THD, THRHLF, TWO, TWOTHD, VCUNIF(2), VXUNIF(2), ZETA 1675 1676 ! Lower bounds of density and its gradient to avoid divisions by zero 1677 PARAMETER ( DENMIN = 1.D-12 ) 1678 PARAMETER ( GDMIN = 1.D-12 ) 1679 1680 ! Fix some numerical parameters 1681 PARAMETER ( FOUTHD=4.D0/3.D0, HALF=0.5D0, & 1682 THD=1.D0/3.D0, THRHLF=1.5D0, & 1683 TWO=2.D0, TWOTHD=2.D0/3.D0 ) 1684 PARAMETER ( TEN81 = 10.0d0/81.0d0 ) 1685 1686 ! Fix some more numerical constants 1687 PI = 4 * ATAN(1.D0) 1688 BETA = 0.066725D0 1689 GAMMA = (1 - LOG(TWO)) / PI**2 1690 MU = BETA * PI**2 / 3 1691 KAPPA = 0.804D0 1692 CWC = 0.0079325D0 1693 1694 ! Translate density and its gradient to new variables 1695 IF (nspin .EQ. 1) THEN 1696 D(1) = HALF * Dens(1) 1697 D(2) = D(1) 1698 DT = MAX( DENMIN, Dens(1) ) 1699 DO IX = 1,3 1700 GD(IX,1) = HALF * GDens(IX,1) 1701 GD(IX,2) = GD(IX,1) 1702 GDT(IX) = GDens(IX,1) 1703 end DO 1704 ELSE 1705 D(1) = Dens(1) 1706 D(2) = Dens(2) 1707 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 1708 DO IX = 1,3 1709 GD(IX,1) = GDens(IX,1) 1710 GD(IX,2) = GDens(IX,2) 1711 GDT(IX) = GDens(IX,1) + GDens(IX,2) 1712 end DO 1713 ENDIF 1714 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 1715 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 1716 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 1717 GDMT = MAX( GDMIN, GDMT ) 1718 1719 ! Find local correlation energy and potential 1720 CALL PW92C( 2, D, ECUNIF, VCUNIF ) 1721 1722 ! Find total correlation energy 1723 RS = ( 3 / (4*PI*DT) )**THD 1724 KF = (3 * PI**2 * DT)**THD 1725 KS = SQRT( 4 * KF / PI ) 1726 ZETA = ( D(1) - D(2) ) / DT 1727 ZETA = MAX( -1.D0+DENMIN, ZETA ) 1728 ZETA = MIN( 1.D0-DENMIN, ZETA ) 1729 PHI = HALF * ( (1+ZETA)**TWOTHD + (1-ZETA)**TWOTHD ) 1730 T = GDMT / (2 * PHI * KS * DT) 1731 F1 = ECUNIF / GAMMA / PHI**3 1732 F2 = EXP(-F1) 1733 A = BETA / GAMMA / (F2-1) 1734 F3 = T**2 + A * T**4 1735 F4 = BETA/GAMMA * F3 / (1 + A*F3) 1736 H = GAMMA * PHI**3 * LOG( 1 + F4 ) 1737 FC = ECUNIF + H 1738 1739 ! Find correlation energy derivatives 1740 DRSDD = - (THD * RS / DT) 1741 DKFDD = THD * KF / DT 1742 DKSDD = HALF * KS * DKFDD / KF 1743 DZDD(1) = 1 / DT - ZETA / DT 1744 DZDD(2) = - (1 / DT) - ZETA / DT 1745 DPDZ = HALF * TWOTHD * ( 1/(1+ZETA)**THD - 1/(1-ZETA)**THD ) 1746 DO IS = 1,2 1747 DECUDD = ( VCUNIF(IS) - ECUNIF ) / DT 1748 DPDD = DPDZ * DZDD(IS) 1749 DTDD = (- T) * ( DPDD/PHI + DKSDD/KS + 1/DT ) 1750 DF1DD = F1 * ( DECUDD/ECUNIF - 3*DPDD/PHI ) 1751 DF2DD = (- F2) * DF1DD 1752 DADD = (- A) * DF2DD / (F2-1) 1753 DF3DD = (2*T + 4*A*T**3) * DTDD + DADD * T**4 1754 DF4DD = F4 * ( DF3DD/F3 - (DADD*F3+A*DF3DD)/(1+A*F3) ) 1755 DHDD = 3 * H * DPDD / PHI 1756 DHDD = DHDD + GAMMA * PHI**3 * DF4DD / (1+F4) 1757 DFCDD(IS) = VCUNIF(IS) + H + DT * DHDD 1758 1759 DO IX = 1,3 1760 DTDGD = (T / GDMT) * GDT(IX) / GDMT 1761 DF3DGD = DTDGD * ( 2 * T + 4 * A * T**3 ) 1762 DF4DGD = F4 * DF3DGD * ( 1/F3 - A/(1+A*F3) ) 1763 DHDGD = GAMMA * PHI**3 * DF4DGD / (1+F4) 1764 DFCDGD(IX,IS) = DT * DHDGD 1765 end DO 1766 end DO 1767 1768 ! Find exchange energy and potential 1769 FX = 0 1770 DO IS = 1,2 1771 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 1772 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 1773 KFS = (3 * PI**2 * DS(IS))**THD 1774 S = GDMS / (2 * KFS * DS(IS)) 1775 ! 1776 ! For PBE: 1777 ! 1778 ! x = MU * S**2 1779 ! dxds = 2*MU*S 1780 ! 1781 ! Wu-Cohen form: 1782 ! 1783 XWC= TEN81 * s**2 + (MU- TEN81) * & 1784 S**2 * exp(-S**2) + log(1+ CWC * S**4) 1785 DXWCDS = 2 * TEN81 * S + (MU - TEN81) * exp(-S**2) * & 1786 2*S * (1 - S*S) + 4 * CWC * S**3 / (1 + CWC * S**4) 1787 1788 F1 = 1 + XWC / KAPPA 1789 F = 1 + KAPPA - KAPPA / F1 1790 ! 1791 ! Note nspin=1 in call to exchng... 1792 ! 1793 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 1794 FX = FX + DS(IS) * EXUNIF * F 1795 1796 DKFDD = THD * KFS / DS(IS) 1797 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 1798 DF1DD = DXWCDS * DSDD / KAPPA 1799 DFDD = KAPPA * DF1DD / F1**2 1800 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 1801 1802 DO IX = 1,3 1803 GDS = 2 * GD(IX,IS) 1804 DSDGD = (S / GDMS) * GDS / GDMS 1805 DF1DGD = DXWCDS * DSDGD / KAPPA 1806 DFDGD = KAPPA * DF1DGD / F1**2 1807 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 1808 end DO 1809 end DO 1810 FX = HALF * FX / DT 1811 1812 ! Set output arguments 1813 EX = FX 1814 EC = FC 1815 DO IS = 1,nspin 1816 DEXDD(IS) = DFXDD(IS) 1817 DECDD(IS) = DFCDD(IS) 1818 DO IX = 1,3 1819 DEXDGD(IX,IS) = DFXDGD(IX,IS) 1820 DECDGD(IX,IS) = DFCDGD(IX,IS) 1821 end DO 1822 end DO 1823 1824 END SUBROUTINE WCXC 1825 1826 1827 1828 SUBROUTINE AM05XC( IREL, nspin, Dens, GDens, & 1829 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 1830 1831 ! ********************************************************************* 1832 ! Implements the Armiento Mattsson AM05 GGA. 1833 ! Ref: R. Armiento and A. E. Mattsson, PRB 72, 085108 (2005) 1834 ! Written by L.C.Balbas and J.M.Soler originally for PBE. December 1996. 1835 ! Modified by J.D. Gale for AM05. May 2009. 1836 ! ******** INPUT ****************************************************** 1837 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 1838 ! INTEGER nspin : Number of spin polarizations (1 or 2) 1839 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 1840 ! spin electron density (if nspin=2) 1841 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 1842 ! ******** OUTPUT ***************************************************** 1843 ! REAL*8 EX : Exchange energy density 1844 ! REAL*8 EC : Correlation energy density 1845 ! REAL*8 DEXDD(nspin) : Partial derivative 1846 ! d(DensTot*Ex)/dDens(ispin), 1847 ! where DensTot = Sum_ispin( Dens(ispin) ) 1848 ! For a constant density, this is the 1849 ! exchange potential 1850 ! REAL*8 DECDD(nspin) : Partial derivative 1851 ! d(DensTot*Ec)/dDens(ispin), 1852 ! where DensTot = Sum_ispin( Dens(ispin) ) 1853 ! For a constant density, this is the 1854 ! correlation potential 1855 ! REAL*8 DEXDGD(3,nspin): Partial derivative 1856 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 1857 ! REAL*8 DECDGD(3,nspin): Partial derivative 1858 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 1859 ! ********* UNITS **************************************************** 1860 ! Lengths in Bohr 1861 ! Densities in electrons per Bohr**3 1862 ! Energies in Hartrees 1863 ! Gradient vectors in cartesian coordinates 1864 ! ********* ROUTINES CALLED ****************************************** 1865 ! am05wbs 1866 ! ******************************************************************** 1867 1868 use precision, only : dp 1869 use am05, only : am05wbs 1870 1871 integer irel, nspin 1872 real(dp) Dens(nspin), DECDD(nspin), DECDGD(3,nspin), & 1873 DEXDD(nspin), DEXDGD(3,nspin), GDens(3,nspin) 1874 1875 ! Internal variables 1876 integer :: is, ix 1877 1878 real(dp) & 1879 D(2), DENMIN, DFXDD(2), DFCDD(2), DFCDGD(3,2), & 1880 DFXDGD(3,2), DFXDG(2), DFCDG(2), & 1881 DS(2), DT, EC, EX, FX, FC, & 1882 GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3) 1883 1884 ! Lower bounds of density and its gradient to avoid divisions by zero 1885 parameter ( DENMIN = 1.D-12 ) 1886 parameter ( GDMIN = 1.D-12 ) 1887 1888 ! Translate density and its gradient to new variables 1889 if (nspin .eq. 1) then 1890 D(1) = 0.5_dp*Dens(1) 1891 D(2) = D(1) 1892 DT = max( DENMIN, Dens(1) ) 1893 do ix = 1,3 1894 GD(ix,1) = 0.5_dp*GDens(ix,1) 1895 GD(ix,2) = GD(ix,1) 1896 GDT(ix) = GDens(ix,1) 1897 enddo 1898 else 1899 D(1) = Dens(1) 1900 D(2) = Dens(2) 1901 DT = max( DENMIN, Dens(1)+Dens(2) ) 1902 do ix = 1,3 1903 GD(ix,1) = GDens(ix,1) 1904 GD(ix,2) = GDens(ix,2) 1905 GDT(ix) = GDens(ix,1) + GDens(ix,2) 1906 enddo 1907 endif 1908 GDM(1) = sqrt( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 1909 GDM(2) = sqrt( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 1910 GDMT = sqrt( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 1911 GDMT = max( GDMIN, GDMT ) 1912 1913 D(1) = max(D(1),denmin) 1914 D(2) = max(D(2),denmin) 1915 1916 ! Call AM05 subroutine 1917 call am05wbs(D(1), D(2), GDM(1), GDM(2), FX, FC, & 1918 DFXDD(1), DFXDD(2), DFCDD(1), DFCDD(2), & 1919 DFXDG(1), DFXDG(2), DFCDG(1), DFCDG(2)) 1920 1921 ! Convert gradient terms into vectors 1922 do is = 1,nspin 1923 do ix = 1,3 1924 DFXDGD(ix,is) = DFXDG(is)*GD(ix,is) 1925 DFCDGD(ix,is) = DFCDG(is)*GD(ix,is) 1926 enddo 1927 enddo 1928 1929 ! Convert FX to form required by SIESTA - note factor of 1/2 1930 ! is already applied in am05 code. 1931 FX = FX / DT 1932 1933 ! Set output arguments 1934 EX = FX 1935 EC = FC 1936 do is = 1,nspin 1937 DEXDD(is) = DFXDD(is) 1938 DECDD(is) = DFCDD(is) 1939 do ix = 1,3 1940 DEXDGD(ix,is) = DFXDGD(ix,is) 1941 DECDGD(ix,is) = DFCDGD(ix,is) 1942 enddo 1943 enddo 1944 1945 END SUBROUTINE AM05XC 1946 1947 1948 1949 SUBROUTINE PW86formX( a, b, c, iRel, nSpin, Dens, GDens, & 1950 EX, dEXdD, dEXdGD ) 1951 1952 ! ********************************************************************* 1953 ! Implements Perdew-Wang-86 Generalized-Gradient-Approximation exchange-only 1954 ! functional form, eps_x=eps_x_LDA*(1+15*a*s**2+b*s**4+c*s**6)**(1/15), 1955 ! but with variable values for parameters a, b, and c 1956 ! Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986) 1957 ! E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009) 1958 ! Written by J.M.Soler. March 2010. 1959 ! ******** INPUT ****************************************************** 1960 ! REAL*8 a, b, c : Parameter a, b, and c of the PW86 functional 1961 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 1962 ! INTEGER nspin : Number of spin polarizations (1 or 2) 1963 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 1964 ! spin electron density (if nspin=2) 1965 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 1966 ! ******** OUTPUT ***************************************************** 1967 ! REAL*8 EX : Exchange energy density 1968 ! REAL*8 DEXDD(nspin) : Partial derivative 1969 ! d(DensTot*Ex)/dDens(ispin), 1970 ! where DensTot = Sum_ispin( Dens(ispin) ) 1971 ! For a constant density, this is the 1972 ! exchange potential 1973 ! REAL*8 DEXDGD(3,nspin): Partial derivative 1974 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 1975 ! ********* UNITS **************************************************** 1976 ! Lengths in Bohr 1977 ! Densities in electrons per Bohr**3 1978 ! Energies in Hartrees 1979 ! Gradient vectors in cartesian coordinates 1980 ! ********* ROUTINES CALLED ****************************************** 1981 ! EXCHNG 1982 ! ******************************************************************** 1983 1984 ! Input 1985 real(dp),intent(in) :: a ! Parameter of PW86 functional 1986 real(dp),intent(in) :: b ! Parameter of PW86 functional 1987 real(dp),intent(in) :: c ! Parameter of PW86 functional 1988 integer, intent(in) :: iRel ! Relativistic exchange? 0=no, 1=yes 1989 integer, intent(in) :: nSpin ! Number of spin components 1990 real(dp),intent(in) :: Dens(nSpin) ! Local electron (spin) density 1991 real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density 1992 1993 ! Output 1994 real(dp),intent(out):: EX ! Exchange energy per electron 1995 real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 1996 real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 1997 1998 ! Internal variables 1999 INTEGER :: IS, IX 2000 real(dp) & 2001 D(2), DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFXDD(2), DFXDGD(3,2), & 2002 DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, EXUNIF, F, F1, FX, & 2003 GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 2004 KFS, PI, S, VXUNIF(2), ZETA 2005 2006 ! Lower bounds of density and its gradient to avoid divisions by zero 2007 PARAMETER ( DENMIN = 1.D-12 ) 2008 PARAMETER ( GDMIN = 1.D-12 ) 2009 2010 ! Fix some more numerical constants 2011 PI = 4 * ATAN(1.D0) 2012 2013 ! Translate density and its gradient to new variables 2014 IF (nspin .EQ. 1) THEN 2015 D(1) = Dens(1) / 2 2016 D(2) = D(1) 2017 DT = MAX( DENMIN, Dens(1) ) 2018 DO IX = 1,3 2019 GD(IX,1) = GDens(IX,1) / 2 2020 GD(IX,2) = GD(IX,1) 2021 GDT(IX) = GDens(IX,1) 2022 end DO 2023 ELSE 2024 D(1) = Dens(1) 2025 D(2) = Dens(2) 2026 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 2027 DO IX = 1,3 2028 GD(IX,1) = GDens(IX,1) 2029 GD(IX,2) = GDens(IX,2) 2030 GDT(IX) = GDens(IX,1) + GDens(IX,2) 2031 end DO 2032 ENDIF 2033 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 2034 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 2035 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 2036 GDMT = MAX( GDMIN, GDMT ) 2037 2038 ! Find exchange energy and potential 2039 FX = 0 2040 DO IS = 1,2 2041 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 2042 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 2043 KFS = (3 * PI**2 * DS(IS))**(1._dp/3) 2044 S = GDMS / (2 * KFS * DS(IS)) 2045 F1 = 1 + 15*a*S**2 + b*S**4 + c*S**6 2046 F = F1**(1._dp/15) 2047 ! 2048 ! Note nspin=1 in call to exchng... 2049 ! 2050 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 2051 FX = FX + DS(IS) * EXUNIF * F 2052 2053 DKFDD = KFS / DS(IS) / 3 2054 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 2055 DF1DS = 30*a*S + 4*b*S**3 + 6*c*S**5 2056 DFDS = F/F1/15 * DF1DS 2057 DFDD = DFDS * DSDD 2058 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 2059 2060 DO IX = 1,3 2061 GDS = 2 * GD(IX,IS) 2062 DSDGD = (S / GDMS) * GDS / GDMS 2063 DFDGD = DFDS * DSDGD 2064 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 2065 end DO 2066 end DO 2067 FX = FX / DT / 2 2068 2069 ! Set output arguments 2070 EX = FX 2071 DO IS = 1,nspin 2072 DEXDD(IS) = DFXDD(IS) 2073 DO IX = 1,3 2074 DEXDGD(IX,IS) = DFXDGD(IX,IS) 2075 end DO 2076 end DO 2077 2078 END SUBROUTINE PW86formX 2079 2080 2081 2082 SUBROUTINE PW86X( IREL, nspin, Dens, GDens, & 2083 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2084 2085 ! ********************************************************************* 2086 ! Implements Perdew-Wang-86 Generalized-Gradient-Approximation 2087 ! exchange-only functional. Correlation energy returns as zero. 2088 ! Ref: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986) 2089 ! Written by J.M.Soler. March 2010. 2090 ! ******** INPUT ****************************************************** 2091 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2092 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2093 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2094 ! spin electron density (if nspin=2) 2095 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2096 ! ******** OUTPUT ***************************************************** 2097 ! REAL*8 EX : Exchange energy density 2098 ! REAL*8 EC : Correlation energy density 2099 ! REAL*8 DEXDD(nspin) : Partial derivative 2100 ! d(DensTot*Ex)/dDens(ispin), 2101 ! where DensTot = Sum_ispin( Dens(ispin) ) 2102 ! For a constant density, this is the 2103 ! exchange potential 2104 ! REAL*8 DECDD(nspin) : Partial derivative 2105 ! d(DensTot*Ec)/dDens(ispin), 2106 ! where DensTot = Sum_ispin( Dens(ispin) ) 2107 ! For a constant density, this is the 2108 ! correlation potential 2109 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2110 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2111 ! REAL*8 DECDGD(3,nspin): Partial derivative 2112 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2113 ! ********* UNITS **************************************************** 2114 ! Lengths in Bohr 2115 ! Densities in electrons per Bohr**3 2116 ! Energies in Hartrees 2117 ! Gradient vectors in cartesian coordinates 2118 ! ********* ROUTINES CALLED ****************************************** 2119 ! PW86formX 2120 ! ******************************************************************** 2121 2122 ! Passed arguments 2123 integer, intent(in) :: IREL, NSPIN 2124 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2125 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2126 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2127 2128 ! Internal parameters of the PW86 exchange functional 2129 real(dp),parameter:: a = 0.0864_dp 2130 real(dp),parameter:: b = 14.0_dp 2131 real(dp),parameter:: c = 0.2_dp 2132 2133 ! Call PW86 routine with appropriate values for a, b, and c. 2134 CALL PW86formX( a, b, c, IREL, nspin, Dens, GDens, & 2135 EX, DEXDD, DEXDGD ) 2136 2137 ! Set correlation energy and derivatives to zero 2138 EC = 0 2139 DECDD(:) = 0 2140 DECDGD(:,:) = 0 2141 2142 END SUBROUTINE PW86X 2143 2144 2145 2146 SUBROUTINE PW86RX( IREL, nspin, Dens, GDens, & 2147 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2148 2149 ! ********************************************************************* 2150 ! Implements Perdew-Wang-86 Generalized-Gradient-Approximation 2151 ! exchange-only functional with the refitted parameters of 2152 ! Murray, Lee, and Langreth. Correlation energy returns as zero. 2153 ! Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986) 2154 ! E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009) 2155 ! Written by J.M.Soler. March 2010. 2156 ! ******** INPUT ****************************************************** 2157 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2158 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2159 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2160 ! spin electron density (if nspin=2) 2161 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2162 ! ******** OUTPUT ***************************************************** 2163 ! REAL*8 EX : Exchange energy density 2164 ! REAL*8 EC : Correlation energy density 2165 ! REAL*8 DEXDD(nspin) : Partial derivative 2166 ! d(DensTot*Ex)/dDens(ispin), 2167 ! where DensTot = Sum_ispin( Dens(ispin) ) 2168 ! For a constant density, this is the 2169 ! exchange potential 2170 ! REAL*8 DECDD(nspin) : Partial derivative 2171 ! d(DensTot*Ec)/dDens(ispin), 2172 ! where DensTot = Sum_ispin( Dens(ispin) ) 2173 ! For a constant density, this is the 2174 ! correlation potential 2175 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2176 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2177 ! REAL*8 DECDGD(3,nspin): Partial derivative 2178 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2179 ! ********* UNITS **************************************************** 2180 ! Lengths in Bohr 2181 ! Densities in electrons per Bohr**3 2182 ! Energies in Hartrees 2183 ! Gradient vectors in cartesian coordinates 2184 ! ********* ROUTINES CALLED ****************************************** 2185 ! PW86formX 2186 ! ******************************************************************** 2187 2188 ! Passed arguments 2189 integer, intent(in) :: IREL, NSPIN 2190 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2191 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2192 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2193 2194 ! Internal parameters of the refitted PW86 exchange functional 2195 real(dp),parameter:: a = 0.1234_dp 2196 real(dp),parameter:: b = 17.33_dp 2197 real(dp),parameter:: c = 0.163_dp 2198 2199 ! Call PW86 routine with appropriate values for a, b, and c. 2200 CALL PW86formX( a, b, c, IREL, nspin, Dens, GDens, & 2201 EX, DEXDD, DEXDGD ) 2202 2203 ! Set correlation energy and derivatives to zero 2204 EC = 0 2205 DECDD(:) = 0 2206 DECDGD(:,:) = 0 2207 2208 END SUBROUTINE PW86RX 2209 2210 2211 SUBROUTINE BHX( IREL, nspin, Dens, GDens, & 2212 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2213 2214 ! ********************************************************************* 2215 ! Implements the combination by Berland and Hyldgaard of Perdew-Wang-86 2216 ! Generalized-Gradient-Approximation exchange-only functional with the 2217 ! refitted parameters of Murray, Lee, and Langreth and Langreth-Vosko 2218 ! screened exchange. Correlation energy returns as zero. 2219 ! Refs: J.P.Perdew & Y.Wang, PRB 33, 8800 (1986) 2220 ! E.D.Murray, K.Lee & D.C.Langreth, JCTC 5, 2754 (2009) 2221 ! K.Berland & P.Hyldgaard, PRB 89, 035412 (2014) 2222 ! Written by Michelle Fritz Feb. 2014. 2223 ! ******** INPUT ****************************************************** 2224 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2225 ! INTEGER NSPIN : Number of spin polarizations (1 or 2) 2226 ! REAL*8 DENS(nspin) : Total electron density (if nspin=1) or 2227 ! spin electron density (if nspin=2) 2228 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2229 ! ******** OUTPUT ***************************************************** 2230 ! REAL*8 EX : Exchange energy density 2231 ! REAL*8 EC : Correlation energy density 2232 ! REAL*8 DEXDD(nspin) : Partial derivative 2233 ! d(DensTot*Ex)/dDens(ispin), 2234 ! where DensTot = Sum_ispin( Dens(ispin) ) 2235 ! For a constant density, this is the 2236 ! exchange potential 2237 ! REAL*8 DECDD(nspin) : Partial derivative 2238 ! d(DensTot*Ec)/dDens(ispin), 2239 ! where DensTot = Sum_ispin( Dens(ispin) ) 2240 ! For a constant density, this is the 2241 ! correlation potential 2242 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2243 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2244 ! REAL*8 DECDGD(3,nspin): Partial derivative 2245 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2246 ! ********* UNITS **************************************************** 2247 ! Lengths in Bohr 2248 ! Densities in electrons per Bohr**3 2249 ! Energies in Hartrees 2250 ! Gradient vectors in cartesian coordinates 2251 ! ********* ROUTINES CALLED ****************************************** 2252 ! EXCHNG 2253 ! ******************************************************************** 2254 2255 ! Passed arguments 2256 integer, intent(in) :: IREL, NSPIN 2257 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2258 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2259 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2260 2261 ! Internal variables 2262 INTEGER :: IS, IX 2263 real(dp) & 2264 D(2), DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFLVDS, DFPW86RDS, & 2265 DFXDD(2), DFXDGD(3,2), DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, & 2266 EXUNIF, FPW86R, FLV, F, F1, FX, GD(3,2), GDM(2), GDMIN, GDMS, & 2267 GDMT, GDS, GDT(3), KFS, PI, S, VXUNIF(2), ZETA, MLV 2268 2269 ! Internal parameters of the refitted PW86 exchange functional 2270 real(dp),parameter:: a = 0.1234_dp 2271 real(dp),parameter:: b = 17.33_dp 2272 real(dp),parameter:: c = 0.163_dp 2273 real(dp),parameter:: zab = -0.8491_dp 2274 real(dp),parameter:: alpha = 0.02178_dp 2275 real(dp),parameter:: beta = 1.15_dp 2276 2277 ! Lower bounds of density and its gradient to avoid divisions by zero 2278 PARAMETER ( DENMIN = 1.D-12 ) 2279 PARAMETER ( GDMIN = 1.D-12 ) 2280 2281 ! Fix some more numerical constants 2282 PI = 4 * ATAN(1.D0) 2283 2284 ! Translate density and its gradient to new variables 2285 IF (NSPIN .EQ. 1) THEN 2286 D(1) = DENS(1) / 2 2287 D(2) = D(1) 2288 DT = MAX( DENMIN, DENS(1) ) 2289 DO IX = 1,3 2290 GD(IX,1) = GDENS(IX,1) / 2 2291 GD(IX,2) = GD(IX,1) 2292 GDT(IX) = GDENS(IX,1) 2293 end DO 2294 ELSE 2295 D(1) = DENS(1) 2296 D(2) = DENS(2) 2297 DT = MAX( DENMIN, DENS(1)+DENS(2) ) 2298 DO IX = 1,3 2299 GD(IX,1) = GDENS(IX,1) 2300 GD(IX,2) = GDENS(IX,2) 2301 GDT(IX) = GDENS(IX,1) + GDENS(IX,2) 2302 end DO 2303 ENDIF 2304 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 2305 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 2306 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 2307 GDMT = MAX( GDMIN, GDMT ) 2308 2309 ! Find exchange energy and potential 2310 FX = 0 2311 DO IS = 1,2 2312 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 2313 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 2314 KFS = (3 * PI**2 * DS(IS))**(1._dp/3) 2315 S = GDMS / (2 * KFS * DS(IS)) 2316 F1 = 1 + 15*a*S**2 + b*S**4 + c*S**6 2317 FPW86R = F1**(1._dp/15) 2318 MLV = -zab/9._dp 2319 FLV = 1 + MLV*S**2 2320 F = (1 / (1 + alpha*S**6)) * FLV 2321 F = F + ((alpha*S**6) / (beta + alpha*S**6)) * FPW86R 2322 2323 ! 2324 ! Note nspin=1 in call to exchng... 2325 ! 2326 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 2327 FX = FX + DS(IS) * EXUNIF * F 2328 2329 DKFDD = KFS / DS(IS) / 3 2330 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 2331 DF1DS = 30*a*S + 4*b*S**3 + 6*c*S**5 2332 DFLVDS = 2*MLV*S 2333 DFPW86RDS = FPW86R/F1/15 * DF1DS 2334 DFDS = -((6*alpha*S**5) / (1 + alpha*S**6)**2) * FLV 2335 DFDS = DFDS + (1 / (1 + alpha*S**6)) * DFLVDS 2336 DFDS = DFDS + ((6*alpha*S**5) / (beta + alpha*S**6)) * FPW86R 2337 DFDS =DFDS-((6*alpha**2*S**11)/(beta + alpha*S**6)**2)*FPW86R 2338 DFDS = DFDS + ((alpha*S**6) / (beta + alpha*S**6)) * DFPW86RDS 2339 DFDD = DFDS * DSDD 2340 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 2341 2342 DO IX = 1,3 2343 GDS = 2 * GD(IX,IS) 2344 DSDGD = (S / GDMS) * GDS / GDMS 2345 DFDGD = DFDS * DSDGD 2346 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 2347 end DO 2348 end DO 2349 FX = FX / DT / 2 2350 2351 EX = FX 2352 DO IS = 1,NSPIN 2353 DEXDD(IS) = DFXDD(IS) 2354 DO IX = 1,3 2355 DEXDGD(IX,IS) = DFXDGD(IX,IS) 2356 end DO 2357 end DO 2358 2359 ! Set correlation energy and derivatives to zero 2360 EC = 0 2361 DECDD(:) = 0 2362 DECDGD(:,:) = 0 2363 2364 END SUBROUTINE BHX 2365 2366 2367 SUBROUTINE B88formX( beta, mu, c, iRel, nSpin, Dens, GDens, & 2368 EX, dEXdD, dEXdGD ) 2369 2370 ! ********************************************************************* 2371 ! Implements Becke-88 Generalized-Gradient-Approximation exchange-only 2372 ! functional form, eps_x=eps_x_LDA*(1+mu*s**2/(1+beta*s*asinh(c*s))), 2373 ! but with variable values for parameters beta, mu, and c 2374 ! Refs: A.D.Becke, PRA 38, 3098 (1988) 2375 ! J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009) 2376 ! Written by J.M.Soler. April 2010. 2377 ! ******** INPUT ****************************************************** 2378 ! REAL*8 beta, mu, c : Parameter beta mu, and c of the B88 functional, 2379 ! as formulated by Klimes et al 2380 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2381 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2382 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2383 ! spin electron density (if nspin=2) 2384 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2385 ! ******** OUTPUT ***************************************************** 2386 ! REAL*8 EX : Exchange energy density 2387 ! REAL*8 DEXDD(nspin) : Partial derivative 2388 ! d(DensTot*Ex)/dDens(ispin), 2389 ! where DensTot = Sum_ispin( Dens(ispin) ) 2390 ! For a constant density, this is the 2391 ! exchange potential 2392 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2393 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2394 ! ********* UNITS **************************************************** 2395 ! Lengths in Bohr 2396 ! Densities in electrons per Bohr**3 2397 ! Energies in Hartrees 2398 ! Gradient vectors in cartesian coordinates 2399 ! ********* ROUTINES CALLED ****************************************** 2400 ! EXCHNG 2401 ! ******************************************************************** 2402 2403 ! Input 2404 real(dp),intent(in) :: beta ! Parameter of B88 functional 2405 real(dp),intent(in) :: mu ! Parameter of B88 functional 2406 real(dp),intent(in) :: c ! Parameter of B88 functional 2407 integer, intent(in) :: iRel ! Relativistic exchange? 0=no, 1=yes 2408 integer, intent(in) :: nSpin ! Number of spin components 2409 real(dp),intent(in) :: Dens(nSpin) ! Local electron (spin) density 2410 real(dp),intent(in) :: GDens(3,nSpin) ! Gradient of electron density 2411 2412 ! Output 2413 real(dp),intent(out):: EX ! Exchange energy per electron 2414 real(dp),intent(out):: dEXdD(nSpin) ! dEx/dDens, Ex=Int(dens*epsX) 2415 real(dp),intent(out):: dEXdGD(3,nSpin) ! dEx/dGrad(Dens) 2416 2417 ! Internal variables 2418 INTEGER :: IS, IX 2419 real(dp) & 2420 ASINHCS, CS, D(2), DASINHDS, & 2421 DENMIN, DF1DS, DFDD, DFDGD, DFDS, DFXDD(2), DFXDGD(3,2), & 2422 DKFDD, DS(2), DSDD, DSDGD, DT, ECUNIF, EXUNIF, F, F1, FX, & 2423 GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 2424 KFS, PI, S, VXUNIF(2), ZETA 2425 2426 ! Lower bounds of density and its gradient to avoid divisions by zero 2427 PARAMETER ( DENMIN = 1.D-12 ) 2428 PARAMETER ( GDMIN = 1.D-12 ) 2429 2430 ! Fix some more numerical constants 2431 PI = 4 * ATAN(1.D0) 2432 2433 ! Translate density and its gradient to new variables 2434 IF (nspin .EQ. 1) THEN 2435 D(1) = Dens(1) / 2 2436 D(2) = D(1) 2437 DT = MAX( DENMIN, Dens(1) ) 2438 DO IX = 1,3 2439 GD(IX,1) = GDens(IX,1) / 2 2440 GD(IX,2) = GD(IX,1) 2441 GDT(IX) = GDens(IX,1) 2442 end DO 2443 ELSE 2444 D(1) = Dens(1) 2445 D(2) = Dens(2) 2446 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 2447 DO IX = 1,3 2448 GD(IX,1) = GDens(IX,1) 2449 GD(IX,2) = GDens(IX,2) 2450 GDT(IX) = GDens(IX,1) + GDens(IX,2) 2451 end DO 2452 ENDIF 2453 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 2454 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 2455 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 2456 GDMT = MAX( GDMIN, GDMT ) 2457 2458 ! Find exchange energy and potential 2459 FX = 0 2460 DO IS = 1,2 2461 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 2462 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 2463 KFS = (3 * PI**2 * DS(IS))**(1._dp/3) 2464 S = GDMS / (2 * KFS * DS(IS)) 2465 CS = C * S 2466 ! Next line is ASINH(CS). See Numerical Recipes 2467 ASINHCS = log(CS+sqrt(CS**2+1)) 2468 F1 = 1 + beta * S * ASINHCS 2469 F = 1 + mu * S**2 / F1 2470 ! 2471 ! Note nspin=1 in call to exchng... 2472 ! 2473 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 2474 FX = FX + DS(IS) * EXUNIF * F 2475 2476 DKFDD = KFS / DS(IS) / 3 2477 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 2478 DASINHDS = C * (1 + CS/sqrt(CS**2+1)) / (CS+sqrt(CS**2+1)) 2479 DF1DS = beta * ASINHCS + beta * S * DASINHDS 2480 DFDS = 2*mu*S/F1 - mu*S**2/F1**2 * DF1DS 2481 DFDD = DFDS * DSDD 2482 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 2483 2484 DO IX = 1,3 2485 GDS = 2 * GD(IX,IS) 2486 DSDGD = (S / GDMS) * GDS / GDMS 2487 DFDGD = DFDS * DSDGD 2488 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 2489 end DO 2490 end DO 2491 FX = FX / DT / 2 2492 2493 ! Set output arguments 2494 EX = FX 2495 DO IS = 1,nspin 2496 DEXDD(IS) = DFXDD(IS) 2497 DO IX = 1,3 2498 DEXDGD(IX,IS) = DFXDGD(IX,IS) 2499 end DO 2500 end DO 2501 2502 END SUBROUTINE B88formX 2503 2504 2505 2506 SUBROUTINE B88X( IREL, nspin, Dens, GDens, & 2507 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2508 2509 ! ********************************************************************* 2510 ! Implements Becke-88 Generalized-Gradient-Approximation 2511 ! exchange-only functional. Correlation energy returns as zero. 2512 ! Refs: A.D.Becke, PRA 38, 3098 (1988) 2513 ! J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009) 2514 ! Written by J.M.Soler. April 2010. 2515 ! ******** INPUT ****************************************************** 2516 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2517 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2518 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2519 ! spin electron density (if nspin=2) 2520 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2521 ! ******** OUTPUT ***************************************************** 2522 ! REAL*8 EX : Exchange energy density 2523 ! REAL*8 EC : Correlation energy density 2524 ! REAL*8 DEXDD(nspin) : Partial derivative 2525 ! d(DensTot*Ex)/dDens(ispin), 2526 ! where DensTot = Sum_ispin( Dens(ispin) ) 2527 ! For a constant density, this is the 2528 ! exchange potential 2529 ! REAL*8 DECDD(nspin) : Partial derivative 2530 ! d(DensTot*Ec)/dDens(ispin), 2531 ! where DensTot = Sum_ispin( Dens(ispin) ) 2532 ! For a constant density, this is the 2533 ! correlation potential 2534 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2535 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2536 ! REAL*8 DECDGD(3,nspin): Partial derivative 2537 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2538 ! ********* UNITS **************************************************** 2539 ! Lengths in Bohr 2540 ! Densities in electrons per Bohr**3 2541 ! Energies in Hartrees 2542 ! Gradient vectors in cartesian coordinates 2543 ! ********* ROUTINES CALLED ****************************************** 2544 ! PW86formX 2545 ! ******************************************************************** 2546 2547 ! Passed arguments 2548 integer, intent(in) :: IREL, NSPIN 2549 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2550 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2551 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2552 2553 ! Internal variables and arrays 2554 real(dp):: beta, c, mu, pi 2555 2556 ! Internal parameters of the B88 exchange functional, written as 2557 ! by Klimes et al 2558 pi = acos(-1._dp) 2559 c = 2*(6*pi**2)**(1._dp/3) 2560 mu = 0.2743_dp 2561 beta = 9 * mu * (6/pi)**(1._dp/3) / (2*c) 2562 2563 ! Call PW86 routine with appropriate values for a, b, and c. 2564 CALL B88formX( beta, mu, c, IREL, nspin, Dens, GDens, & 2565 EX, DEXDD, DEXDGD ) 2566 2567 ! Set correlation energy and derivatives to zero 2568 EC = 0 2569 DECDD(:) = 0 2570 DECDGD(:,:) = 0 2571 2572 END SUBROUTINE B88X 2573 2574 2575 2576 SUBROUTINE B88KBMX( IREL, nspin, Dens, GDens, & 2577 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2578 2579 ! ********************************************************************* 2580 ! Implements Becke-88 GGA exchange functional, reparametrized by Klimes 2581 ! et al for its combined use with vdW-DF nonlocal correlation. 2582 ! Correlation energy returns as zero. 2583 ! Refs: A.D.Becke, PRA 38, 3098 (1988) 2584 ! J.Klimes, D.R.Bowler, and A.Michaelides, JPCM 22, 022201 (2009) 2585 ! Written by J.M.Soler. April 2010. 2586 ! ******** INPUT ****************************************************** 2587 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2588 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2589 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2590 ! spin electron density (if nspin=2) 2591 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2592 ! ******** OUTPUT ***************************************************** 2593 ! REAL*8 EX : Exchange energy density 2594 ! REAL*8 EC : Correlation energy density 2595 ! REAL*8 DEXDD(nspin) : Partial derivative 2596 ! d(DensTot*Ex)/dDens(ispin), 2597 ! where DensTot = Sum_ispin( Dens(ispin) ) 2598 ! For a constant density, this is the 2599 ! exchange potential 2600 ! REAL*8 DECDD(nspin) : Partial derivative 2601 ! d(DensTot*Ec)/dDens(ispin), 2602 ! where DensTot = Sum_ispin( Dens(ispin) ) 2603 ! For a constant density, this is the 2604 ! correlation potential 2605 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2606 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2607 ! REAL*8 DECDGD(3,nspin): Partial derivative 2608 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2609 ! ********* UNITS **************************************************** 2610 ! Lengths in Bohr 2611 ! Densities in electrons per Bohr**3 2612 ! Energies in Hartrees 2613 ! Gradient vectors in cartesian coordinates 2614 ! ********* ROUTINES CALLED ****************************************** 2615 ! PW86formX 2616 ! ******************************************************************** 2617 2618 ! Passed arguments 2619 integer, intent(in) :: IREL, NSPIN 2620 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2621 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2622 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2623 2624 ! Internal variables and arrays 2625 real(dp):: beta, c, mu, pi 2626 2627 ! Klimes et al parameters for the B88 exchange functional 2628 pi = acos(-1._dp) 2629 c = 2*(6*pi**2)**(1._dp/3) 2630 mu = 0.22_dp 2631 beta = mu / 1.2_dp 2632 2633 ! Call PW86 routine with appropriate values for a, b, and c. 2634 CALL B88formX( beta, mu, c, IREL, nspin, Dens, GDens, & 2635 EX, DEXDD, DEXDGD ) 2636 2637 ! Set correlation energy and derivatives to zero 2638 EC = 0 2639 DECDD(:) = 0 2640 DECDGD(:,:) = 0 2641 2642 END SUBROUTINE B88KBMX 2643 2644 2645 SUBROUTINE C09X( IREL, nspin, Dens, GDens, & 2646 EX, EC, DEXDD, DECDD, DEXDGD, DECDGD ) 2647 2648 ! ********************************************************************* 2649 ! Implements Cooper-09 exchange-only GGA (for use with vdW-DF1 correlation). 2650 ! Correlation energy returns as zero. 2651 ! Ref: V.R.Cooper, PRB 81, 161104(R) (2010) 2652 ! Written by J.M.Soler. April 2010. 2653 ! ******** INPUT ****************************************************** 2654 ! INTEGER IREL : Relativistic-exchange switch (0=No, 1=Yes) 2655 ! INTEGER nspin : Number of spin polarizations (1 or 2) 2656 ! REAL*8 Dens(nspin) : Total electron density (if nspin=1) or 2657 ! spin electron density (if nspin=2) 2658 ! REAL*8 GDens(3,nspin) : Total or spin density gradient 2659 ! ******** OUTPUT ***************************************************** 2660 ! REAL*8 EX : Exchange energy density 2661 ! REAL*8 EC : Correlation energy density 2662 ! REAL*8 DEXDD(nspin) : Partial derivative 2663 ! d(DensTot*Ex)/dDens(ispin), 2664 ! where DensTot = Sum_ispin( Dens(ispin) ) 2665 ! For a constant density, this is the 2666 ! exchange potential 2667 ! REAL*8 DECDD(nspin) : Partial derivative 2668 ! d(DensTot*Ec)/dDens(ispin), 2669 ! where DensTot = Sum_ispin( Dens(ispin) ) 2670 ! For a constant density, this is the 2671 ! correlation potential 2672 ! REAL*8 DEXDGD(3,nspin): Partial derivative 2673 ! d(DensTot*Ex)/d(GradDens(i,ispin)) 2674 ! REAL*8 DECDGD(3,nspin): Partial derivative 2675 ! d(DensTot*Ec)/d(GradDens(i,ispin)) 2676 ! ********* UNITS **************************************************** 2677 ! Lengths in Bohr 2678 ! Densities in electrons per Bohr**3 2679 ! Energies in Hartrees 2680 ! Gradient vectors in cartesian coordinates 2681 ! ********* ROUTINES CALLED ****************************************** 2682 ! none 2683 ! ******************************************************************** 2684 2685 ! Passed arguments 2686 integer, intent(in) :: IREL, NSPIN 2687 real(dp),intent(in) :: DENS(NSPIN), GDENS(3,NSPIN) 2688 real(dp),intent(out):: EX, EC, DECDD(NSPIN), DECDGD(3,NSPIN), & 2689 DEXDD(NSPIN), DEXDGD(3,NSPIN) 2690 2691 ! Internal variables and arrays (V.R.Cooper, PRB 81, 161104(R) (2010)) 2692 real(dp),parameter:: alpha = 0.0483_dp 2693 real(dp),parameter:: kappa = 1.245_dp 2694 real(dp),parameter:: mu = 0.0617 2695 2696 ! Internal variables 2697 INTEGER :: IS, IX 2698 real(dp) & 2699 D(2), DENMIN, DES2DS, DF1DS, DFDD, DFDGD, DFDS, & 2700 DFXDD(2), DFXDGD(3,2), DKFDD, DS(2), DSDD, DSDGD, DT, & 2701 ECUNIF, ES2, EXUNIF, F, FX, & 2702 GD(3,2), GDM(2), GDMIN, GDMS, GDMT, GDS, GDT(3), & 2703 KFS, PI, S, S2, VXUNIF(2) 2704 2705 ! Lower bounds of density and its gradient to avoid divisions by zero 2706 PARAMETER ( DENMIN = 1.D-12 ) 2707 PARAMETER ( GDMIN = 1.D-12 ) 2708 2709 ! Fix some more numerical constants 2710 PI = 4 * ATAN(1.D0) 2711 2712 ! Translate density and its gradient to new variables 2713 IF (nspin .EQ. 1) THEN 2714 D(1) = Dens(1) / 2 2715 D(2) = D(1) 2716 DT = MAX( DENMIN, Dens(1) ) 2717 DO IX = 1,3 2718 GD(IX,1) = GDens(IX,1) / 2 2719 GD(IX,2) = GD(IX,1) 2720 GDT(IX) = GDens(IX,1) 2721 end DO 2722 ELSE 2723 D(1) = Dens(1) 2724 D(2) = Dens(2) 2725 DT = MAX( DENMIN, Dens(1)+Dens(2) ) 2726 DO IX = 1,3 2727 GD(IX,1) = GDens(IX,1) 2728 GD(IX,2) = GDens(IX,2) 2729 GDT(IX) = GDens(IX,1) + GDens(IX,2) 2730 end DO 2731 ENDIF 2732 GDM(1) = SQRT( GD(1,1)**2 + GD(2,1)**2 + GD(3,1)**2 ) 2733 GDM(2) = SQRT( GD(1,2)**2 + GD(2,2)**2 + GD(3,2)**2 ) 2734 GDMT = SQRT( GDT(1)**2 + GDT(2)**2 + GDT(3)**2 ) 2735 GDMT = MAX( GDMIN, GDMT ) 2736 2737 ! Find exchange energy and potential 2738 FX = 0 2739 DO IS = 1,2 2740 DS(IS) = MAX( DENMIN, 2 * D(IS) ) 2741 GDMS = MAX( GDMIN, 2 * GDM(IS) ) 2742 KFS = (3 * PI**2 * DS(IS))**(1._dp/3) 2743 S = GDMS / (2 * KFS * DS(IS)) 2744 S2 = S**2 2745 2746 ! Next lines are the core of the C09 functional 2747 ES2 = exp(-alpha*S2/2) 2748 F = 1 + mu*S2*ES2**2 + kappa*(1-ES2) 2749 DES2DS = -alpha*S*ES2 2750 DFDS = 2*mu*S*ES2**2 + 2*mu*S2*ES2*DES2DS - kappa*DES2DS 2751 2752 ! Note nspin=1 in call to exchng... 2753 CALL EXCHNG( IREL, 1, DS(IS), EXUNIF, VXUNIF(IS) ) 2754 FX = FX + DS(IS) * EXUNIF * F 2755 2756 DKFDD = KFS / DS(IS) / 3 2757 DSDD = S * ( -(DKFDD/KFS) - 1/DS(IS) ) 2758 DFDD = DFDS * DSDD 2759 DFXDD(IS) = VXUNIF(IS) * F + DS(IS) * EXUNIF * DFDD 2760 2761 DO IX = 1,3 2762 GDS = 2 * GD(IX,IS) 2763 DSDGD = (S / GDMS) * GDS / GDMS 2764 DFDGD = DFDS * DSDGD 2765 DFXDGD(IX,IS) = DS(IS) * EXUNIF * DFDGD 2766 end DO 2767 end DO 2768 FX = FX / DT / 2 2769 2770 ! Set output arguments 2771 EX = FX 2772 DO IS = 1,nspin 2773 DEXDD(IS) = DFXDD(IS) 2774 DO IX = 1,3 2775 DEXDGD(IX,IS) = DFXDGD(IX,IS) 2776 end DO 2777 end DO 2778 ! Set correlation energy and derivatives to zero 2779 EC = 0 2780 DECDD(:) = 0 2781 DECDGD(:,:) = 0 2782 2783 END SUBROUTINE C09X 2784 2785END MODULE m_ggaxc 2786