1 Subroutine xc_hcth(tol_rho, xfac, lxfac, nlxfac, 2 , cfac, lcfac, nlcfac,rho, delrho, 3 & Amat, Cmat, nq, ipol, Ex, Ec, qwght, 4 & ldew,func,funcname) 5c 6c$Id$ 7c 8 Implicit none 9c 10#include "dft2drv.fh" 11c 12 logical ldew ! [input] 13 logical lcfac, nlcfac, lxfac, nlxfac ! [input] 14 double precision func(*) ![input/output] 15 double precision cfac, xfac ![input] 16 character*4 funcname ! functional name [input] 17c 18 integer ipol ! no. of spin states [input] 19 integer nq ! no. of quadrature pts [input] 20 double precision tol_rho! [input]!threshold on density 21 double precision Ec ! Correlation energy [input/output] 22 double precision Ex ! Exchange energy [input/output] 23 double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input] 24 double precision delrho(nq,3,ipol) ! Charge Density Gradient[input] 25 double precision qwght(nq) ! Quadrature Weights [input] 26 double precision Amat(nq,ipol) !Sampling Matrices for the XC [output] 27 double precision Cmat(nq,*)!Potential & Energy [output] 28c 29c References: 30c F.A.Hamprecht, A.J.Cohen, D.J.Tozer and N.C.Handy, 31c J. Chem. Phys. 109, 6264-6271 (1998) 32c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch 33c 34 integer n 35 double precision gammaval 36c to hcth 37 double precision rhoa 38 double precision rhob 39 double precision za 40 double precision zb 41 double precision hE_x, hE_c 42 double precision dfdrax, dfdrac 43 double precision dfdrbx, dfdrbc 44 double precision dfdzac,dfdzax,dfdza 45 double precision dfdzbc,dfdzbx,dfdzb 46c 47 if(ipol.eq.1) then 48 do n=1,nq 49 if(rho(n,1).gt.tol_rho) then 50 rhoa=0.d0 51 rhoa=0.5d0*rho(n,1) 52 rhob=0.d0 53 rhob=rhoa 54 za=0.d0 55 gammaval=0.25d0*(delrho(n,1,1)*delrho(n,1,1) + 56 & delrho(n,2,1)*delrho(n,2,1) + 57 & delrho(n,3,1)*delrho(n,3,1)) 58 if(gammaval.gt.tol_rho**2) za=sqrt(gammaval) 59 zb=za 60 call hcth(ipol,funcname, 61 * dfdrax,dfdrac, dfdzax,dfdzac, 62 * dfdrbx,dfdrbc, dfdzbx,dfdzbc, 63 . rhoa, rhob, 64 & za, zb, hE_x, hE_c, tol_rho) 65 if(ldew) func(n)=func(n)+hE_c*cfac+hE_x*xfac 66 Ec=Ec+hE_c*qwght(n)*cfac 67 Ex=Ex+hE_x*qwght(n)*xfac 68 Amat(n,1) = Amat(n,1)+dfdrac*cfac+dfdrax*xfac 69 dfdza=(dfdzac*cfac+dfdzax*xfac)*0.5d0 70 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza 71 endif 72 enddo 73 else 74 do n=1,nq 75 if(rho(n,1).gt.tol_rho) then 76 rhoa=rho(n,2) 77 rhob=rho(n,3) 78 za=0.d0 79 gammaval=delrho(n,1,1)*delrho(n,1,1) + 80 & delrho(n,2,1)*delrho(n,2,1) + 81 & delrho(n,3,1)*delrho(n,3,1) 82 if(gammaval.gt.tol_rho**2) za=sqrt(gammaval) 83 zb=0.d0 84 gammaval=delrho(n,1,2)*delrho(n,1,2) + 85 & delrho(n,2,2)*delrho(n,2,2) + 86 & delrho(n,3,2)*delrho(n,3,2) 87 if(gammaval.gt.tol_rho**2) zb=sqrt(gammaval) 88 call hcth(ipol,funcname, 89 * dfdrax,dfdrac, dfdzax,dfdzac, 90 * dfdrbx,dfdrbc, dfdzbx,dfdzbc, 91 . rhoa, rhob, 92 & za, zb, hE_x, hE_c, tol_rho) 93 if(ldew) func(n)=func(n)+hE_c*cfac+hE_x*xfac 94 Ec=Ec+hE_c*qwght(n)*cfac 95 Ex=Ex+hE_x*qwght(n)*xfac 96 Amat(n,1) = Amat(n,1)+dfdrac*cfac+dfdrax*xfac 97 Amat(n,2) = Amat(n,2)+dfdrbc*cfac+dfdrbx*xfac 98 dfdza=dfdzac*cfac+dfdzax*xfac 99 dfdzb=dfdzbc*cfac+dfdzbx*xfac 100 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza*0.5d0 101 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dfdzb*0.5d0 102 endif 103 enddo 104 endif 105 return 106 end 107 108 Subroutine xc_hcth_d2(tol_rho, xfac, lxfac, nlxfac, 109 & cfac, lcfac, nlcfac,rho, delrho, 110 & Amat, Amat2, Cmat, Cmat2, nq, ipol, 111 & Ex, Ec, qwght, ldew,func,funcname) 112c 113c$Id$ 114c 115 Implicit none 116#include "errquit.fh" 117c 118#include "dft2drv.fh" 119#include "2ndDerivB97.h" 120#include "mafdecls.fh" 121c 122 logical ldew ! [input] 123 logical lcfac, nlcfac, lxfac, nlxfac ! [input] 124 double precision func(*) ![input/output] 125 double precision cfac, xfac ![input] 126 character*4 funcname ! functional name [input] 127c 128 integer ipol ! no. of spin states [input] 129 integer nq ! no. of quadrature pts [input] 130 double precision tol_rho! [input]!threshold on density 131 double precision Ec ! Correlation energy [input/output] 132 double precision Ex ! Exchange energy [input/output] 133 double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input] 134 double precision delrho(nq,3,ipol) ! Charge Density Gradient[input] 135 double precision qwght(nq) ! Quadrature Weights [input] 136 double precision Amat(nq,ipol) !Sampling Matrices for the XC [output] 137 double precision Cmat(nq,*)!Potential & Energy [output] 138 double precision Amat2(nq,NCOL_AMAT2) ! XC functional seconds [output] 139 double precision Cmat2(nq,NCOL_CMAT2) ! XC functional seconds [output] 140c 141 integer i,max_pow_u 142 double precision rho_a(0:3),rho_b(0:3) 143 double precision FX(0:_FXC_NUMDERI),FC(0:_FXC_NUMDERI) 144#include "xc_hcth.fh" 145 double precision sol((limpow+1)*3) 146 call xc_htch_loadcf(funcname,sol,max_pow_u) 147c if(max_pow_u.gt.2) call errquit(' not ready for maxpow=4',0,0) 148 do i=1,nq 149 if(ipol.eq.2) then 150 rho_a(0)=rho(i,2) 151 rho_b(0)=rho(i,3) 152 rho_a(1)=delrho(i,1,1)**2 153 rho_a(2)=delrho(i,2,1)**2 154 rho_a(3)=delrho(i,3,1)**2 155 rho_b(1)=delrho(i,1,2)**2 156 rho_b(2)=delrho(i,2,2)**2 157 rho_b(3)=delrho(i,3,2)**2 158 else 159 rho_a(0)=rho(i,1)*.5d0 160 rho_b(0)=rho_a(0) 161 rho_a(1)=(delrho(i,1,1)*.5d0)**2 162 rho_a(2)=(delrho(i,2,1)*.5d0)**2 163 rho_a(3)=(delrho(i,3,1)*.5d0)**2 164 rho_b(1)=rho_a(1) 165 rho_b(2)=rho_a(2) 166 rho_b(3)=rho_a(3) 167 endif 168 if(rho_a(0).gt.tol_rho.or.rho_b(0).gt.tol_rho) then 169 call dft_xckernel_xb97(rho_a, rho_b, 1d0, tol_rho, 170 F FX, max_pow_u, sol) 171 call dscal(_FXC_NUMDERI+1,xfac,FX,1) 172 Ex = Ex + FX(_FXC_E)*qwght(i) 173 if(ldew) func(i)=func(i) + FX(_FXC_E) 174 call xc_fxc2acmat(nq,i,ipol,FX, 175 A Amat,Cmat,Amat2,Cmat2) 176 call dft_xckernel_cb97(rho_a, rho_b, 1d0, tol_rho, 177 F FC, max_pow_u, sol) 178 call dscal(_FXC_NUMDERI+1,cfac,FC,1) 179 Ec = Ec + FC(_FXC_E)*qwght(i) 180 if(ldew) func(i)=func(i) + FC(_FXC_E) 181 call xc_fxc2acmat(nq,i,ipol,FC, 182 A Amat,Cmat,Amat2,Cmat2) 183 endif 184 enddo 185 186 return 187 end 188 Subroutine xc_hcth_d2num(tol_rho, xfac, lxfac, nlxfac, 189 & cfac, lcfac, nlcfac,rho, delrho, 190 & Amat, Amat2, Cmat, Cmat2, nq, ipol, 191 & Ex, Ec, qwght, ldew,func,funcname) 192c 193c$Id$ 194c 195 Implicit none 196#include "errquit.fh" 197c 198#include "dft2drv.fh" 199#include "mafdecls.fh" 200c 201 logical ldew ! [input] 202 logical lcfac, nlcfac, lxfac, nlxfac ! [input] 203 double precision func(*) ![input/output] 204 double precision cfac, xfac ![input] 205 character*4 funcname ! functional name [input] 206c 207 integer ipol ! no. of spin states [input] 208 integer nq ! no. of quadrature pts [input] 209 double precision tol_rho! [input]!threshold on density 210 double precision Ec ! Correlation energy [input/output] 211 double precision Ex ! Exchange energy [input/output] 212 double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input] 213 double precision delrho(nq,3,ipol) ! Charge Density Gradient[input] 214 double precision qwght(nq) ! Quadrature Weights [input] 215 double precision Amat(nq,ipol) !Sampling Matrices for the XC [output] 216 double precision Cmat(nq,*)!Potential & Energy [output] 217 double precision Amat2(nq,NCOL_AMAT2) ! XC functional seconds [output] 218 double precision Cmat2(nq,NCOL_CMAT2) ! XC functional seconds [output] 219c 220c Local variables 221c 222 integer l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc, 223 & i_qwght_copy, npert 224 double precision ExDum, EcDum 225c 226c First get the functional and first derivative values 227c 228 call xc_hcth(tol_rho, xfac, lxfac, nlxfac, cfac, lcfac, nlcfac, 229 & rho, delrho, Amat, Cmat, nq, ipol, Ex, Ec, qwght, ldew, func, 230 & funcname) 231c 232c Compute the second derivative values by finite difference 233c 234 call xc_setup_fd(tol_rho, rho, delrho, qwght, nq, ipol, .true., 235 & l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc, 236 & i_qwght_copy) 237c 238c Compute functional first derivatives at perturbed density parameter 239c values - note that the number of points is nq*2*npert and that the 240c routine is called as unrestricted 241c 242 npert = 5 243 call xc_hcth(tol_rho, xfac, lxfac, nlxfac, cfac, lcfac, nlcfac, 244 & dbl_mb(i_prho), dbl_mb(i_pdelrho), dbl_mb(i_pAmat), 245 & dbl_mb(i_pCmat), nq*2*npert, ipol, ExDum, EcDum, 246 & dbl_mb(i_qwght_copy), ldew, dbl_mb(i_pfunc), funcname) 247 call xc_make_fd(Amat2, Cmat2, nq, .true., dbl_mb(i_pAmat), 248 & dbl_mb(i_pCmat)) 249c 250c Free temporary storage allocated by xc_setup_fd 251c 252 if (.not.ma_free_heap(l_storage)) 253 & call errquit('xc_hcth_d2: cannot pop stack',0, MA_ERR) 254c 255 return 256 end 257 258 SUBROUTINE hcth(ipol,functional, 259 * dfdrax,dfdrac, dfdzax,dfdzac, 260 * dfdrbx,dfdrbc, dfdzbx,dfdzbc, 261 1 rhoa, rhob, za, zb, hE_x, hE_c, 262 2 tol_rho) 263 264c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch 265C SUPPLIED TO THE ROUTINE: 266C 267C rhoa -- value of rhoalpha at a given grid point 268C rhob -- value of rhobeta at a given grid point 269C za -- zeta_alpha, as defined in the TH1 paper (JCP 108 2545), 270C that is mod(grad(rhoalpha)), a scalar quantity. 271C zb -- mod(grad(rhobeta)) 272C zab -- zeta_{alpha beta} as defined in the TH1 paper, that is 273C grad(rhoalpha).grad(rhobeta) 274C energy -- a boolean variable deciding whether to compute the energy 275C contribution at the point in space (true) or the 276C appropriate derivatives (false) needed for the KS matrix 277C _and_ the energy contribution. 278 279C RETURNED FROM THE ROUTINE: 280 281C hE_x -- the contribution to the energy at this point in space. 282C hE_c -- the contribution to the energy at this point in space. 283C dfdra -- partial functional derivative of F_xc with respect to 284C rhoalpha 285C dfdrb -- partial functional derivative of F_xc with respect to 286C rhobeta 287C dfdza -- partial functional derivative of F_xc with respect to 288C mod(grad(rhoalpha)), divided by za !!!!!!!!! 289C i.e. 1 d f 290C --- * ---- 291C za d za 292C This is a consequence of the Cadpac implementation 293C dfdzb -- partial functional derivative of F_xc with respect to 294C mod(grad(rhobeta)), divided by zb !!!!!!!!! 295 296 297 implicit none 298#include "errquit.fh" 299 300 double precision rhoa ![input] 301 double precision rhob ![input] 302 double precision za ![input] 303 double precision zb ![input] 304 integer ipol ![input] 305 double precision hE_x ![output] 306 double precision hE_c ![output] 307 double precision tol_rho ![input] 308 double precision dfdrax ![output] 309 double precision dfdrac ![output] 310 double precision dfdrbx ![output] 311 double precision dfdrbc ![output] 312 double precision dfdzax ![output] 313 double precision dfdzac ![output] 314 double precision dfdzbx ![output] 315 double precision dfdzbc ![output] 316 character*4 functional 317 double precision pi 318 PARAMETER (PI=3.1415926535898D0) 319#include "xc_hcth.fh" 320c 321c variables passed to hcderiv 322c 323 integer nofunc,max_pow_u 324 double precision sol((limpow+1)*3), F((limpow+1)*3,4), 325C & FF((limpow+1)*3,5,4), 326 & F_xc((limpow+1)*3) 327 328Cfah sol -- contains the coefficients of the terms in F_xc 329Cfah convention: sol(1) = c_{x alpha, 0}, c_{x beta, 0} 330Cfah sol(2) = c_{c alpha alpha, 0}, c_{c beta beta, 0} 331Cfah sol(3) = c_{c alpha beta, 0} 332Cfah sol(4) = c_{x alpha, 1}, c_{x beta, 1} 333Cfah sol(5) = c_{c alpha alpha, 1}, c_{c beta beta, 1} 334Cfah sol(6) = c_{c alpha beta, 1} 335Cfah 336Cfah etc. 337Cfah 338Cfah f(5) -- contains the partial first functional derivatives of F_xc with 339Cfah respect to 340Cfah the four quantities (IN THIS ORDER): ra, rb, za, zb 341Cfah 342Cfah ff(5,5) contains the second derivatives with 343Cfah respect to the same five quantities 344 345Cfah F_xa -- contains the alpha exchange bit containing the various powers 346Cfah of u_{x alpha} (eq. (18) of Becke V paper) 347Cfah F_xb -- beta 348Cfah u_{x beta} 349Cfah F_caa -- contains the alpha parallel spin correlation bit with the powers 350Cfah of u_{c alpha alpha} 351Cfah F_cbb -- beta 352Cfah u_{c beta beta} 353Cfah F_cab -- contains the anti-parallel spin correlation bit with the powers 354Cfah of u_{c alpha beta} 355 356C Initialise 357 358 dfdrac = 0.D0 359 dfdrax = 0.D0 360 dfdrbc = 0.D0 361 dfdrbx = 0.D0 362 dfdzac = 0.D0 363 dfdzax = 0.D0 364 dfdzbc = 0.D0 365 dfdzbx = 0.D0 366 hE_c = 0.D0 367 hE_x = 0.D0 368 369 IF (rhoa .LT. tol_rho.and.rhob.lt.tol_rho) RETURN 370Cfah numerical cutoff: if the density is too low, its contribution is 371Cfah neglectable. 372 call xc_htch_loadcf(functional,sol,max_pow_u) 373 374 CALL hcderiv(max_pow_u,ipol, 375 & F, 376CFF, 377 & F_xc, 378 & rhoa,rhob,za,zb, 379 & sol,tol_rho) 380 381c if(ipol.eq.2) then 382c DO n = 1, (max_pow_u+1)*3 383c dfdra = dfdra + F(n,1) 384c dfdrb = dfdrb + F(n,2) 385c if(za.gt.tol_rho) dfdza = dfdza + F(n,3) / za 386c if(zb.gt.tol_rho) dfdzb = dfdzb + F(n,4) / zb 387c ENDDO 388c else 389c DO n = 1, (max_pow_u+1)*3 390c dfdra = dfdra + F(n,1) 391c enddo 392c if(za.gt.tol_rho) then 393c DO n = 1, (max_pow_u+1)*3 394c dfdza = dfdza + F(n,3) / za 395c enddo 396c endif 397c endif 398Cfah big thanks to NCH: cadpac requires df/(za * dza), NOT 399Cfah df/dza 400cDEC$ NOVECTOR 401 DO n = 0, max_pow_u 402 hE_x = hE_x + F_xc (n*3 + 1) 403 hE_c = hE_c + F_xc (n*3 + 2) + F_xc (n*3 + 3) 404 dfdrax = dfdrax + F(n*3+1,1) 405 dfdrac = dfdrac + F(n*3+2,1) + F(n*3+3,1) 406 if(za.gt.tol_rho) then 407 dfdzax = dfdzax + F(n*3+1,3) / za 408 dfdzac = dfdzac + (F(n*3+2,3)+F(n*3+3,3)) / za 409 endif 410 if(ipol.eq.2) then 411 dfdrbx = dfdrbx + F(n*3+1,2) 412 dfdrbc = dfdrbc + F(n*3+2,2) + F(n*3+3,2) 413 if(zb.gt.tol_rho) then 414 dfdzbx = dfdzbx + F(n*3+1,4) / zb 415 dfdzbc = dfdzbc + (F(n*3+2,4)+F(n*3+3,4)) / zb 416 endif 417 endif 418 ENDDO 419 RETURN 420 END 421 SUBROUTINE hcderiv(max_pow_u,ipol, 422 & F, 423CFF, 424 & F_xc, 425 & rhoa,rhob,za,zb, 426 & sol,tol_rho) 427 428c subroutine supplied by Fred Hamprecht fah@igc.phys.chem.ethz.ch 429 implicit none 430 INTEGER max_pow_u,ipol 431 integer limpow 432 parameter (limpow=4) 433 double precision f_xc((limpow + 1)*3) 434 double precision f((limpow+1)*3,4) 435C, ff((limpow+1)*3,5,4) 436 double precision rhoa, rhob, za, zb,tol_rho 437 438Cfah COMMON/special/h_atom 439Cfah LOGICAL h_atom 440 DOUBLE PRECISION sol((limpow+1)*3) 441 442 443 DOUBLE PRECISION dF_xa(4) 444 DOUBLE PRECISION dF_xb(4) 445 DOUBLE PRECISION dF_caa(4) 446 DOUBLE PRECISION dF_cbb(4) 447 DOUBLE PRECISION dF_cab(4) 448Cfah these are the first derivatives of the terms of F_xc with respect to 449Cfah the 4 quantities. the index 450Cfah runs over the particular partial derivatives of each term. 451Cfah More explicitly: these are the partial functional derivatives of 452Cfah F_XXX with respect to rhoa, rhob, za and zb. 453 454c DOUBLE PRECISION d2F_xa(4,4) 455c DOUBLE PRECISION d2F_xb(4,4) 456c DOUBLE PRECISION d2F_caa(4,4) 457c DOUBLE PRECISION d2F_cbb(4,4) 458c DOUBLE PRECISION d2F_cab(4,4) 459 460Cfah these are the first derivatives of the different transformed variables 461Cfah u with respect to rhoa, rhob, za and zb. These different derivatives 462Cfah with respect to these 4 quantities named above are stored in these 463Cfah arrays. 464 465 DOUBLE PRECISION Pi 466 PARAMETER (Pi = 3.1415926535898D0) 467 double precision rho 468 DOUBLE PRECISION s_a2, s_b2, s_avg2, u_caa, u_cbb, u_cab 469 DOUBLE PRECISION du_caa_by_drhoa, du_caa_by_dza, du_cbb_by_drhob 470 DOUBLE PRECISION du_cbb_by_dzb, du_cab_by_drhoa, du_cab_by_drhob 471 DOUBLE PRECISION du_cab_by_dza, du_cab_by_dzb 472C, du_caa_by_drhoa_dza 473C DOUBLE PRECISION du_caa_by_dza_dza, du_cbb_by_dzb_dzb 474C DOUBLE PRECISION du_cbb_by_drhob_dzb, du_cab_by_drhoa_dza 475C DOUBLE PRECISION du_cab_by_drhoa_dzb, du_cab_by_drhob_dza 476C DOUBLE PRECISION du_cab_by_drhob_dzb, du_cab_by_dza_dza, 477C ,du_cab_by_dza_dzb 478C DOUBLE PRECISION du_cab_by_dzb_dzb 479 DOUBLE PRECISION rsa, rsa12, rsa32, rsa21, rsb, 480 ,rsb12, rsb32, rsb21 481 DOUBLE PRECISION rsab, rsab12, rsab32, rsab21 482 DOUBLE PRECISION drsa_by_drhoa, drsb_by_drhob, drsab_by_drhoa 483 DOUBLE PRECISION drsab_by_drhob 484 DOUBLE PRECISION zeta, dzeta_by_drhoa, dzeta_by_drhob 485 DOUBLE PRECISION fzeta, dfzeta_by_dzeta, 486 , e_crsa1, e_crsb1 487 DOUBLE PRECISION e_crsab1, e_crsab0, a_crsab 488 DOUBLE PRECISION e_crsabzeta, de_crsa1_by_drsa, de_crsb1_by_drsb 489 DOUBLE PRECISION da_crsab_by_drsab, de_crsab0_by_drsab 490 DOUBLE PRECISION de_crsab1_by_drsab, de_crsabzeta_by_drsab 491 DOUBLE PRECISION de_crsabzeta_by_dzeta, e_caa, e_cbb, e_cab, 492 & de_caa_by_drhoa, de_cbb_by_drhob, de_cab_by_drhoa, 493 & de_cab_by_drhob, 494 & c_naa, c_nbb, c_nab 495 DOUBLE PRECISION F_xs,F_xs0 ! this is a function which is called. 496 DOUBLE PRECISION dF_xs_by_drhos, dF_xs_by_dzs, 497 ,dF_xs_by_drhos0,dF_xs_by_drhos1, dF_xs_by_dzs1 498 INTEGER i, j, n 499 integer n1 500 double precision x1,x2,x3,x4 501 double precision e_crs1 502 double precision drsbydrh 503 double precision d2ez 504 double precision decrsdrs 505C 506C F_xs computes HCTH contribution to exchange Energy 507C using Dirac functional as LDA part 508C usage of F_xs 509C F_xs(n, sol(), rhoa, za) 510C 511 F_xs(n1, x1, x2, x3) = 512 = (-3.D0*x1*(3.D0/Pi)**(1.D0/3.D0)*x2**(4.D0/3.D0)* 513 - ((0.004D0*x3**2.D0)/(0.004D0*x3**2.D0 + 514 - x2**(8.D0/3.D0)))**n1)/(2.D0*2.D0**(2.D0/3.D0)) 515 F_xs0( x1, x2, x3) = 516 = (-3.D0*x1*(3.D0/Pi)**(1.D0/3.D0)*x2**(4.D0/3.D0) 517 - )/(2.D0*2.D0**(2.D0/3.D0)) 518C 519C dF_xs_by_drhos computes dE_x/drho derivative 520Cfah computes the derivative of the term with u^n of the exchange part of 521Cfah F_xc with respect to rho of the same spin. 522Cfah n -- the power of u involved in this term 523Cfah c_xs -- the coefficient c_xs(n) of the term of spin s with the 524Cfah power n of u; is NOT passed over as an array. 525Cfah rhos -- rhosigma, that is, either rhoalpha or rhobeta 526Cfah zs -- mod(grad(rhosigma)), again for alpha or beta 527C usage dF_xs_by_drhos(n, c_xs, rhos, zs) 528C 529 dF_xs_by_drhos(n1, x1, x2, x3) = -(x1*(6.D0/Pi)** 530 * (1.D0/3.D0)*x2**(1.D0/3.D0)*((0.004D0*x3*x3)/ 531 / (x2**(8.D0/3.D0)+0.004D0*x3*x3))**n1)+ 532 + (x1*0.008D0*n1*(6.D0/Pi)**(1.D0/3.D0)* 533 * x2**3.D0*x3*x3*((0.004D0*x3*x3)/ 534 / (x2**(8.D0/3.D0)+0.004D0*x3*x3))**(-1 + n1))/ 535 / (x2**(8.D0/3.D0)+0.004D0*x3*x3)**2 536 dF_xs_by_drhos0(x1, x2, x3) = -(x1*(6.D0/Pi)** 537 * (1.D0/3.D0)*x2**(1.D0/3.D0)) 538 dF_xs_by_drhos1(x1, x2, x3) = -(x1*(6.D0/Pi)** 539 * (1.D0/3.D0)*x2**(1.D0/3.D0)*((0.004D0*x3*x3)/ 540 / (x2**(8.D0/3.D0)+0.004D0*x3*x3)))+ 541 + (x1*0.008D0*(6.D0/Pi)**(1.D0/3.D0)* 542 * x2**3.D0*x3*x3)/ 543 / (x2**(8.D0/3.D0)+0.004D0*x3*x3)**2 544C 545C F_xc with respect to zs 546Cfah see above (function dF_xs_by_drhos) for definition of the 547Cfah other variables 548C usage dF_xs_by_dzs (n, c_xs, rhos, zs) 549c 550 dF_xs_by_dzs(n1, x1, x2, x3) = 551 = (-3.d0*x1*n1*(3.D0/Pi)**(1.D0/3.D0)* 552 * x2**(4.D0/3.D0)*((0.004D0*x3*x3)/ 553 / (x2**(8.D0/3.D0)+ 0.004D0*x3*x3))**(-1+n1)* 554 * ((-0.000032D0*x3*x3*x3)/(x2**(8.D0/3.D0)+ 555 + 0.004D0*x3*x3)**2+(0.008D0*x3)/ 556 / (x2**(8.D0/3.D0)+0.004D0*x3*x3)))/ 557 / (2.D0**(5.D0/3.D0)) 558 dF_xs_by_dzs1(x1, x2, x3) = 559 = (-3.d0*x1*(3.D0/Pi)**(1.D0/3.D0)* 560 * x2**(4.D0/3.D0)* 561 * ((-0.000032D0*x3*x3*x3)/(x2**(8.D0/3.D0)+ 562 + 0.004D0*x3*x3)**2+(0.008D0*x3)/ 563 / (x2**(8.D0/3.D0)+0.004D0*x3*x3)))/ 564 / (2.D0**(5.D0/3.D0)) 565c dF_xs_by_dzs0(x1, x2, x3) = 0d0 566c 567c usage 568c e_crsa1 = e_crs1(rsa12,rsa,rsa32,rsa21) 569c 570 e_crs1(x1,x2,x3,x4) = -0.03108999999999999d0* 571 * dlog(1.d0 + 32.16468317787069d0/ 572 / (14.1189d0*x1+6.1977d0*x2 + 3.3662d0*x3 + 573 + 0.6251699999999999d0*x4))*(1.d0 + 0.20548d0*x2) 574 drsbydrh(x1) = -((1.d0/x1)**(4.D0/3.D0)/ 575 - (6.d0**(2.D0/3.D0)*Pi**(1.D0/3.D0))) 576c usage decrsdrs(rsa,rsa12,rsa21,rsa32) 577 decrsdrs(x1,x2,x3,x4) = ((1.d0 + 0.20548d0*x1)* 578 - (6.1977d0 + 7.05945d0/x2 + 1.25034d0*x1+5.0493d0*x2))/ 579 - ((6.1977d0*x1+14.1189d0*x2+0.6251699999999999d0*x3 + 580 + 3.3662d0*x4)**2d0*(1.d0 + 32.16468317787069d0/ 581 - (6.1977d0*x1+14.1189d0*x2+0.6251699999999999d0*x3+ 582 + 3.3662d0*x4))) - 0.006388373199999999d0* 583 - dlog(1.d0 + 32.16468317787069d0/(6.1977d0*x1 + 14.1189d0*x2 + 584 - 0.6251699999999999d0*x3 + 3.3662d0*x4)) 585c 586 DO j = 1, 4 587 DO n = 1, (max_pow_u+1)*3 588 F(n,j) = 0.D0 589Cfah later on, n has a different meaning: n as power of u, not 590Cfah as number of the coefficient. 591 ENDDO 592 dF_xa(j) = 0.D0 593 dF_xb(j) = 0.D0 594 dF_caa(j) = 0.D0 595 dF_cbb(j) = 0.D0 596 dF_cab(j) = 0.D0 597C DO k = 1, 4 598C DO n = 1, (max_pow_u+1)*3 599C FF(n,j,k) = 0.D0 600C ENDDO 601C d2F_xa(j,k) = 0.D0 602C d2F_xb(j,k) = 0.D0 603C d2F_caa(j,k) = 0.D0 604C d2F_cbb(j,k) = 0.D0 605C d2F_cab(j,k) = 0.D0 606C ENDDO 607 ENDDO 608 DO j = 1, (max_pow_u+1)*3 609 F_xc(j) = 0.D0 610 ENDDO 611 612Cfah -------------------------------------------------------------- 613 614Cfah call the expensive correlation parts here just once, and store their 615Cfah values in a temporary variable. Then compute the actual F_c derivatives 616Cfah with the various powers of u. 617 618 rho = rhoa + rhob 619 620 s_a2=0.d0 621 if(za.gt.tol_rho.and.rhoa.gt.tol_rho) 622 + s_a2 = za**2.D0 / rhoa**(8.D0/3.D0) 623 s_b2=0.d0 624 if(zb.gt.tol_rho.and.rhob.gt.tol_rho) 625 + s_b2 = zb**2.D0 / rhob**(8.D0/3.D0) 626 s_avg2 = 0.5D0*(s_a2 + s_b2) 627 628 u_caa = 0.2D0*s_a2/(1.D0+0.2D0*s_a2) 629 u_cbb = 0.2D0*s_b2/(1.D0+0.2D0*s_b2) 630 631 u_cab = 0.006D0*s_avg2/(1.d0+0.006D0*s_avg2) 632 if(rhoa.gt.tol_rho) then 633 rsa = ((3.d0/Pi)**(1.D0/3.D0)* 634 - (1.d0/rhoa)**(1.D0/3.D0))/ 635 - 2**(2.D0/3.D0) 636 rsa12 = rsa**(1.D0/2.D0) 637 rsa32 = rsa**(3.D0/2.D0) 638 rsa21 = rsa**2.D0 639 else 640 rsa=0d0 641 rsa12=0d0 642 rsa32=0d0 643 rsa21=0d0 644 endif 645 646 if(rhob.gt.tol_rho) then 647 rsb = ((3.d0/Pi)**(1.D0/3.D0)* 648 - (1.d0/rhob)**(1.D0/3.D0))/ 649 - 2**(2.D0/3.D0) 650 rsb12 = rsb**(1.D0/2.D0) 651 rsb32 = rsb**(3.D0/2.D0) 652 rsb21 = rsb**2.D0 653C 654C pw91 LDA Ecorr 655C 656 if(rhob.gt.tol_rho) then 657 call xc_pw91ldag(rsb, 2, e_crsb1 ,de_crsb1_by_drsb, d2ez) 658 else 659 e_crsb1 = 0d0 660 de_crsb1_by_drsb = 0d0 661 endif 662 663 du_cbb_by_drhob = (-1.6D0*zb**2*rhob**(5.D0/3.D0))/ 664 - (3.d0*(0.2D0*zb**2 + 665 - rhob**(8.D0/3.D0))**2) 666 du_cbb_by_dzb = (2*0.2D0*zb*rhob**(8.D0/3.D0))/ 667 - (0.2D0*zb**2 + rhob**(8.D0/3.D0))**2 668 if(rhoa.gt.0d0) then 669 du_cab_by_drhoa = (-16*0.006D0*za*za*rhoa**(5.D0/3.D0)* 670 - rhob**(16.D0/3.D0))/ 671 - (3.d0*(0.006D0*zb**2*rhoa**(8.D0/3.D0) + 672 - 0.006D0*za*za*rhob**(8.D0/3.D0) + 673 - 2.d0*rhoa**(8.D0/3.D0)* 674 - rhob**(8.D0/3.D0))**2) 675 du_cab_by_dza = (4*0.006D0*za*rhoa**(8.D0/3.D0)* 676 - rhob**(16.D0/3.D0))/ 677 - (0.006D0*zb**2*rhoa**(8.D0/3.D0) + 678 - 0.006D0*za*za*rhob**(8.D0/3.D0) + 679 - 2.d0*rhoa**(8.D0/3.D0)* 680 - rhob**(8.D0/3.D0))**2 681 du_cab_by_drhob = (-16*0.006D0*zb**2*rhob**(5.D0/3.D0)* 682 - rhoa**(16.D0/3.D0))/ 683 - (3.d0*(0.006D0*za*za*rhob**(8.D0/3.D0) + 684 - 0.006D0*zb**2*rhoa**(8.D0/3.D0) + 685 - 2.d0*rhob**(8.D0/3.D0)* 686 - rhoa**(8.D0/3.D0))**2) 687 du_cab_by_dzb = (4*0.006D0*zb*rhoa**(16.D0/3.D0)* 688 - rhob**(8.D0/3.D0))/ 689 - (0.006D0*zb**2*rhoa**(8.D0/3.D0) + 690 - 0.006D0*za*za*rhob**(8.D0/3.D0) + 691 - 2.d0*rhoa**(8.D0/3.D0)* 692 - rhob**(8.D0/3.D0))**2 693 else 694 du_cab_by_drhoa = 0d0 695 du_cab_by_dza = 0d0 696 du_cab_by_drhob = 0d0 697 du_cab_by_dzb = 0d0 698 endif 699 drsb_by_drhob = drsbydrh(rhob) 700 else 701 e_crsb1 = 0d0 702 de_crsb1_by_drsb = 0d0 703 du_cbb_by_drhob = 0d0 704 du_cbb_by_dzb = 0d0 705 du_cab_by_drhoa = 0d0 706 du_cab_by_drhob = 0d0 707 du_cab_by_dza = 0d0 708 du_cab_by_dzb = 0d0 709 drsb_by_drhob = 0d0 710 endif 711 712 rsab = ((3.d0/Pi)**(1.D0/3.D0)* 713 - (1.d0/rho)**(1.D0/3.D0))/ 714 - 2**(2.D0/3.D0) 715 rsab12 = rsab**(1.D0/2.D0) 716 rsab32 = rsab**(3.D0/2.D0) 717 rsab21 = rsab**2.D0 718 719 zeta = (rhoa-rhob)/rho 720 if(zeta.lt.-1d0) zeta=-1d0 721 if(zeta.gt.1d0) zeta=1d0 722 723 if(abs(1d0-zeta).gt.tol_rho) then 724 fzeta = (-2.d0 + sign(1d0,1.d0 - zeta)*(abs(1.d0 - zeta))** 725 - (4.D0/3.D0) + 726 - (1.d0 + zeta)**(4.D0/3.D0))/ 727 - (-2.d0 + 2.d0*2.d0**(1.D0/3.D0)) 728 else 729 fzeta = 1d0 730 endif 731 732C 733C pw91 LDA Ecorr 734C 735 if(rhoa.gt.tol_rho) then 736 call xc_pw91ldag(rsa, 2, e_crsa1, 737 D de_crsa1_by_drsa, d2ez) 738 else 739 e_crsa1 = 0d0 740 endif 741 if(rho.gt.tol_rho) then 742 call xc_pw91ldag(rsab, 2, e_crsab1, 743 D de_crsab1_by_drsab, d2ez) 744c eps_c(rs,0) 745 call xc_pw91ldag(rsab, 1, e_crsab0, 746 D de_crsab0_by_drsab, d2ez) 747c alpha(rs) 748 call xc_pw91ldag(rsab, 3, a_crsab, 749 D da_crsab_by_drsab, d2ez) 750 else 751 e_crsab1 = 0d0 752 e_crsab0 = 0d0 753 a_crsab=0d0 754 endif 755 756 e_crsabzeta = e_crsab0+a_crsab*fzeta*(1.d0-zeta**4)/1.709921D0+ 757 - (e_crsab1-e_crsab0)*fzeta*zeta**4 758 759 e_caa = rhoa*e_crsa1 760 e_cbb = rhob*e_crsb1 761 e_cab = rho*e_crsabzeta - rhoa*e_crsa1 - rhob*e_crsb1 762 if(rhoa.gt.tol_rho) then 763 du_caa_by_drhoa = (-1.6D0*za*za*rhoa**(5.D0/3.D0))/ 764 - (3.*(0.2D0*za*za + 765 - rhoa**(8.D0/3.D0))**2) 766 du_caa_by_dza = (2*0.2D0*za*rhoa**(8.D0/3.D0))/ 767 - (0.2D0*za*za + rhoa**(8.D0/3.D0))**2 768 else 769 du_caa_by_drhoa = 0d0 770 du_caa_by_dza = 0d0 771 endif 772 773 774 775Cfah Second derivatives are not required by cadpac. 776Cfah du_caa_by_drhoa_dza = (16*0.2D0*za*rhoa**(5.D0/3.D0)* 777Cfah - (0.2D0*za**2 - rhoa**(8.D0/3.D0)))/ 778Cfah - (3.*(0.2D0*za**2 + 779Cfah - rhoa**(8.D0/3.D0))**3) 780Cfah du_cbb_by_drhob_dzb = (16*0.2D0*zb*rhob**(5.D0/3.D0)* 781Cfah - (0.2D0*zb**2 - rhob**(8.D0/3.D0)))/ 782Cfah - (3.*(0.2D0*zb**2 + 783Cfah - rhob**(8.D0/3.D0))**3) 784Cfah 785Cfah du_caa_by_dza_dza = (2*0.2D0*rhoa**(8.D0/3.D0)* 786Cfah - (-3*0.2D0*za**2 + rhoa**(8.D0/3.D0)) 787Cfah - )/ 788Cfah - (0.2D0*za**2 + rhoa**(8.D0/3.D0))**3 789Cfah du_cbb_by_dzb_dzb = (2*0.2D0*rhob**(8.D0/3.D0)* 790Cfah - (-3*0.2D0*zb**2 + rhob**(8.D0/3.D0)) 791Cfah - )/ 792Cfah - (0.2D0*zb**2 + rhob**(8.D0/3.D0))**3 793Cfah 794Cfah du_cab_by_drhoa_dza = (-32*0.006D0*rhoa**(5.D0/3.D0)* 795Cfah - (0.006D0*za*zb**2* 796Cfah - rhoa**(8.D0/3.D0)* 797Cfah - rhob**(16.D0/3.D0) - 798Cfah - 0.006D0*za**3*rhob**8 + 799Cfah - 2*za*rhoa**(8.D0/3.D0)*rhob**8))/ 800Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) + 801Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 802Cfah - 2*rhoa**(8.D0/3.D0)* 803Cfah - rhob**(8.D0/3.D0))**3) 804Cfah du_cab_by_drhoa_dzb = (64*0.006D0**2*za**2*zb* 805Cfah - rhoa**(13.D0/3.D0)* 806Cfah - rhob**(16.D0/3.D0))/ 807Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) + 808Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 809Cfah - 2*rhoa**(8.D0/3.D0)* 810Cfah - rhob**(8.D0/3.D0))**3) 811Cfah du_cab_by_drhob_dza = (64*0.006D0**2*za*zb**2* 812Cfah - rhoa**(16.D0/3.D0)* 813Cfah - rhob**(13.D0/3.D0))/ 814Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) + 815Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 816Cfah - 2*rhoa**(8.D0/3.D0)* 817Cfah - rhob**(8.D0/3.D0))**3) 818Cfah du_cab_by_drhob_dzb = (-32*0.006D0*rhob**(5.D0/3.D0)* 819Cfah - (-(0.006D0*zb**3*rhoa**8) + 820Cfah - 0.006D0*za**2*zb* 821Cfah - rhoa**(16.D0/3.D0)* 822Cfah - rhob**(8.D0/3.D0) + 823Cfah - 2*zb*rhoa**8*rhob**(8.D0/3.D0)))/ 824Cfah - (3.*(0.006D0*zb**2*rhoa**(8.D0/3.D0) + 825Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 826Cfah - 2*rhoa**(8.D0/3.D0)* 827Cfah - rhob**(8.D0/3.D0))**3) 828Cfah du_cab_by_dza_dza = (4*0.006D0*(0.006D0*zb**2* 829Cfah - rhoa**(16.D0/3.D0)* 830Cfah - rhob**(16.D0/3.D0) - 831Cfah - 3*0.006D0*za**2*rhoa**(8.D0/3.D0)* 832Cfah - rhob**8 + 2*rhoa**(16.D0/3.D0)*rhob**8 833Cfah - ))/ 834Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) + 835Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 836Cfah - 2*rhoa**(8.D0/3.D0)* 837Cfah - rhob**(8.D0/3.D0))**3 838Cfah du_cab_by_dza_dzb = (-16*0.006D0**2*za*zb* 839Cfah - rhoa**(16.D0/3.D0)* 840Cfah - rhob**(16.D0/3.D0))/ 841Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) + 842Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 843Cfah - 2*rhoa**(8.D0/3.D0)* 844Cfah - rhob**(8.D0/3.D0))**3 845Cfah du_cab_by_dzb_dzb = (4*0.006D0*rhoa**(16.D0/3.D0)* 846Cfah - rhob**(8.D0/3.D0)* 847Cfah - (-3*0.006D0*zb**2*rhoa**(8.D0/3.D0) + 848Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 849Cfah - 2*rhoa**(8.D0/3.D0)* 850Cfah - rhob**(8.D0/3.D0)))/ 851Cfah - (0.006D0*zb**2*rhoa**(8.D0/3.D0) + 852Cfah - 0.006D0*za**2*rhob**(8.D0/3.D0) + 853Cfah - 2*rhoa**(8.D0/3.D0)* 854Cfah - rhob**(8.D0/3.D0))**3 855 856 if(rhoa.gt.tol_rho) then 857 drsa_by_drhoa = drsbydrh(rhoa) 858 else 859 drsa_by_drhoa =0d0 860 endif 861 if(rho.gt.tol_rho) then 862 drsab_by_drhoa = drsbydrh(rho) 863 drsab_by_drhob = drsab_by_drhoa 864 else 865 drsab_by_drhoa = 0d0 866 drsab_by_drhob = 0d0 867 endif 868 869 dzeta_by_drhoa = 2.d0*rhob/rho**2 870 dzeta_by_drhob = -2.d0*rhoa/rho**2 871 872 dfzeta_by_dzeta = ((-4.d0*sign(1d0,1.d0 - zeta)* 873 * (abs(1.d0 - zeta))**(1.D0/3.D0))/ 874 - 3.d0 + (4.d0*(1.d0 + zeta)** 875 - (1.D0/3.D0))/3.)/ 876 - (-2.d0 + 2d0*2**(1.D0/3.D0)) 877 878 if(rhoa.gt.tol_rho) then 879 call xc_pw91ldag(rsa, 2, e_crsa1, 880 D de_crsa1_by_drsa, d2ez) 881 else 882 de_crsa1_by_drsa = 0d0 883 de_crsab1_by_drsab = 0d0 884 endif 885 886 887 de_crsabzeta_by_drsab = 1.124999956683108D0*(1.d0 - zeta**4)* 888 - (-2.d0+sign(1d0,1.d0-zeta)*(abs(1.d0-zeta))**(4.D0/3.D0)+ 889 - (1.d0 + zeta)**(4.D0/3.D0))* 890 - da_crsab_by_drsab + 891 - de_crsab0_by_drsab + 892 - fzeta*zeta**4* 893 - (- de_crsab0_by_drsab + 894 - de_crsab1_by_drsab ) 895 896 de_crsabzeta_by_dzeta = 1.499999942244144D0*(-1.d0 + zeta**4)* 897 - (sign(1d0,1.d0 - zeta)*(abs(1.d0 - zeta))**(1.D0/3.D0) - 898 - (1.d0 + zeta)**(1.D0/3.D0))*a_crsab - 899 - 4.499999826732434D0*zeta**3* 900 - (-2.d0 + sign(1d0,1.d0-zeta)*abs(1.d0-zeta)**(4.D0/3.D0)+ 901 - (1.d0 + zeta)**(4.D0/3.D0))*a_crsab + 902 - (2*zeta**4*sign(1d0,1.d0-zeta)*(abs(1.d0-zeta)**(1.D0/3.D0)- 903 - (1.d0 + zeta)**(1.D0/3.D0))* 904 - (e_crsab0 - e_crsab1))/ 905 - (3.*(-1.d0 + 2**(1.D0/3.D0))) + 906 - 4*fzeta*(-e_crsab0 + 907 - e_crsab1)*zeta**3 908 909Cfah this is with application of the chain rule; I keep it that general 910Cfah because this way, I only have to define one "G". 911 de_caa_by_drhoa = e_crsa1 + rhoa*de_crsa1_by_drsa*drsa_by_drhoa 912 de_cbb_by_drhob = e_crsb1 + rhob*de_crsb1_by_drsb*drsb_by_drhob 913 914 de_cab_by_drhoa = -e_crsa1 + 915 - e_crsabzeta - 916 - rhoa*de_crsa1_by_drsa* 917 - drsa_by_drhoa + 918 - rho*(de_crsabzeta_by_drsab* 919 - drsab_by_drhoa + 920 - de_crsabzeta_by_dzeta* 921 - dzeta_by_drhoa) 922 923 de_cab_by_drhob = -e_crsb1 + 924 - e_crsabzeta - 925 - rhob*de_crsb1_by_drsb* 926 - drsb_by_drhob + 927 - rho*(de_crsabzeta_by_drsab* 928 - drsab_by_drhob + 929 - de_crsabzeta_by_dzeta* 930 - dzeta_by_drhob) 931 932 933 934Cfah Here starts the big outer loop over the powers u 935 DO n = 0, max_pow_u 936 c_naa = sol((n*3) + 2) 937 c_nbb = c_naa 938 c_nab = sol((n*3) + 3) 939 940Cfah construction of the F_xc itself 941Cfah ------------------------------- 942 IF (rhoa.GT.tol_rho) THEN 943 if(n.eq.0) then 944 F_xc(1) = F_xs0 (sol(1), rhoa, za) 945 else 946 F_xc(n*3+1) = F_xs (n, sol((n*3) + 1), rhoa, za) 947 endif 948 ENDIF 949 if(u_caa.gt.tol_rho)then 950 if(n.eq.0) then 951 F_xc(2) = e_caa*c_naa 952 else 953 F_xc(n*3+2) = e_caa*u_caa**n*c_naa 954 endif 955 endif 956 957 IF (rhob.GT.tol_rho) THEN 958 if(n.eq.0) then 959 F_xc(1) = F_xc(1)+F_xs0(sol(1), rhob, zb) 960 else 961 F_xc(n*3+1) = F_xc(n*3+1)+F_xs(n, sol((n*3) + 1), rhob, zb) 962 endif 963 ENDIF 964 if(u_cbb.gt.tol_rho) then 965 if(n.eq.0) then 966 F_xc(2) = F_xc(2)+e_cbb*c_nbb 967 else 968 F_xc(n*3+2) = F_xc(n*3+2)+e_cbb*u_cbb**n*c_nbb 969 endif 970 endif 971 972 if(u_cab.gt.tol_rho) then 973 if(n.eq.0) then 974 F_xc(3) = e_cab*c_nab 975 else 976 F_xc(n*3+3) = e_cab*u_cab**n*c_nab 977 endif 978 endif 979 980Cfah print*, 'in deriv:', e_cab, u_cab, c_nab 981 982Cfah First Derivatives 983Cfah --------------------- 984 985 if(za.gt.tol_rho)then 986 if(n.eq.0) then 987 dF_xa(1) = dF_xs_by_drhos0 ( sol(1), rhoa, za) 988 dF_xa(3) = 0d0 989 elseif(n.eq.1) then 990 dF_xa(1) = dF_xs_by_drhos1 (sol(4), rhoa, za) 991 dF_xa(3) = dF_xs_by_dzs1 ( sol(4), rhoa, za) 992 else 993 dF_xa(1) = dF_xs_by_drhos (n, sol((n*3) + 1), rhoa, za) 994 dF_xa(3) = dF_xs_by_dzs (n, sol((n*3) + 1), rhoa, za) 995 endif 996 endif 997 998 if(zb.gt.tol_rho) then 999 if(n.eq.0) then 1000 dF_xb(2) = dF_xs_by_drhos0(sol(1), rhob, zb) 1001 dF_xb(4) = 0d0 1002 elseif(n.eq.1) then 1003 dF_xb(2) = dF_xs_by_drhos1(sol(4), rhob, zb) 1004 dF_xb(4) = dF_xs_by_dzs1 ( sol(4), rhob, zb) 1005 else 1006 dF_xb(2) = dF_xs_by_drhos (n, sol((n*3) + 1), rhob, zb) 1007 dF_xb(4) = dF_xs_by_dzs (n, sol((n*3) + 1), rhob, zb) 1008 endif 1009 endif 1010 1011 if(u_caa.gt.tol_rho) then 1012 if(n.eq.0) then 1013 dF_caa(1) = c_naa*de_caa_by_drhoa 1014 dF_caa(3) = 0d0 1015 elseif(n.eq.1) then 1016 dF_caa(1) = c_naa*u_caa* 1017 * de_caa_by_drhoa+c_naa*e_caa*du_caa_by_drhoa 1018 dF_caa(3) = c_naa*e_caa*du_caa_by_dza 1019 else 1020 dF_caa(1) = c_naa*u_caa**n* 1021 * de_caa_by_drhoa+c_naa*n*e_caa*u_caa**(-1+n)* 1022 * du_caa_by_drhoa 1023 dF_caa(3) = c_naa*n*e_caa*u_caa**(-1+n)*du_caa_by_dza 1024 endif 1025 endif 1026 1027 if(u_cbb.gt.tol_rho) then 1028 if(n.eq.0) then 1029 dF_cbb(2) = c_nbb*de_cbb_by_drhob 1030 dF_cbb(4) = 0d0 1031 elseif(n.eq.1) then 1032 dF_cbb(2) = c_nbb*u_cbb*de_cbb_by_drhob + 1033 - c_nbb*e_cbb*du_cbb_by_drhob 1034 dF_cbb(4) = c_nbb*e_cbb*du_cbb_by_dzb 1035 else 1036 dF_cbb(2) = c_nbb*u_cbb**n*de_cbb_by_drhob + 1037 - c_nbb*n*e_cbb*u_cbb**(-1 + n)*du_cbb_by_drhob 1038 dF_cbb(4) = c_nbb*n*e_cbb*u_cbb**(-1+n)*du_cbb_by_dzb 1039 endif 1040 endif 1041 1042 1043 if(u_cab.gt.tol_rho) then 1044 if(n.eq.0) then 1045 dF_cab(1) = c_nab*de_cab_by_drhoa 1046 dF_cab(2) = c_nab*de_cab_by_drhob 1047 dF_cab(3) = 0d0 1048 dF_cab(4) = 0d0 1049 elseif(n.eq.1) then 1050 dF_cab(1) = c_nab*u_cab* 1051 * de_cab_by_drhoa+c_nab*n*e_cab*du_cab_by_drhoa 1052 dF_cab(2) = c_nab*u_cab* 1053 - de_cab_by_drhob+c_nab*n*e_cab*du_cab_by_drhob 1054 dF_cab(3) = c_nab*e_cab*du_cab_by_dza 1055 dF_cab(4) = c_nab*n*e_cab*du_cab_by_dzb 1056 else 1057 dF_cab(1) = c_nab*u_cab**n* 1058 * de_cab_by_drhoa+c_nab*n*e_cab*u_cab**(-1+n)* 1059 * du_cab_by_drhoa 1060 dF_cab(2) = c_nab*u_cab**n* 1061 - de_cab_by_drhob+c_nab*n*e_cab*u_cab**(-1+n)* 1062 - du_cab_by_drhob 1063 dF_cab(3) = c_nab*n*e_cab* 1064 - u_cab**(-1+n)*du_cab_by_dza 1065 dF_cab(4) = c_nab*n*e_cab* 1066 - u_cab**(-1+n)*du_cab_by_dzb 1067 endif 1068 endif 1069 1070Cfah Second Derivatives 1071Cfah ------------------ 1072 1073Cfah d2F_xa(1,1) = d2F_xs_by_drhos_drhos (n, sol((n*3) + 1), 1074Cfah & rhoa, za) 1075Cfah see comment below, for the (2,2) term. 1076Cfah d2F_xa(1,2) = 0 1077Cfah d2F_xa(1,3) = d2F_xs_by_drhos_dzs (n, sol((n*3) + 1), rhoa, za) 1078Cfah d2F_xa(1,4) = 0 1079Cfah d2F_xa(2,2) = 0 1080Cfah d2F_xa(2,3) = 0 1081Cfah d2F_xa(2,4) = 0 1082Cfah d2F_xa(3,3) = d2F_xs_by_dzs_dzs (n, sol((n*3) + 1), rhoa, za) 1083Cfah d2F_xa(3,4) = 0 1084Cfah d2F_xa(4,4) = 0 1085 1086Cfah for alpha spin, elements are non-zero when both indices are odd; 1087Cfah for beta spin, elements are non-zero when both indices are even. 1088Cfah the matrix is symmetric, and the upper triangle contains the 1089Cfah 10 elements given above and below. 1090 1091Cfah d2F_xb(1,1) = 0 1092Cfah d2F_xb(1,2) = 0 1093Cfah d2F_xb(1,3) = 0 1094Cfah d2F_xb(1,4) = 0 1095Cfah d2F_xb(2,2) = d2F_xs_by_drhos_drhos (n, sol((n*3) + 1), 1096Cfah & rhob, zb) 1097Cfah this term is NOT zero, but needs not be evaluated since we don't 1098Cfah need it for the construction of v (cf. routine "va" in the fit 1099Cfah program) 1100Cfah d2F_xb(2,3) = 0 1101Cfah d2F_xb(2,4) = d2F_xs_by_drhos_dzs (n, sol((n*3) + 1), rhob, zb) 1102Cfah d2F_xb(3,3) = 0 1103Cfah d2F_xb(3,4) = 0 1104Cfah d2F_xb(4,4) = d2F_xs_by_dzs_dzs (n, sol((n*3) + 1), rhob, zb) 1105 1106 1107Cfah d2F_caa(1,1) = !=0, but not needed 1108Cfah d2F_caa(1,2) = 0.D0 (not needed) 1109Cfah d2F_caa(1,3) = c_naa*n*u_caa**(-1 + n)* 1110Cfah - de_caa_by_drhoa* 1111Cfah - du_caa_by_dza + 1112Cfah - c_naa*(-1 + n)*n*e_caa* 1113Cfah - u_caa**(-2 + n)* 1114Cfah - du_caa_by_dza* 1115Cfah - du_caa_by_drhoa + 1116Cfah - c_naa*n*e_caa*u_caa**(-1 + n)* 1117Cfah - du_caa_by_drhoa_dza 1118Cfah d2F_caa(1,4) = 0.D0 1119Cfah d2F_caa(2,2) = 0.D0 (not needed) 1120Cfah d2F_caa(2,3) = 0.D0 1121Cfah d2F_caa(2,4) = 0.D0 1122Cfah d2F_caa(3,3) = c_naa*n*e_caa*u_caa**(-2 + n)* 1123Cfah - ((-1 + n)*du_caa_by_dza** 1124Cfah - 2 + u_caa* 1125Cfah - du_caa_by_dza_dza) 1126Cfah d2F_caa(3,4) = 0.D0 1127Cfah d2F_caa(4,4) = 0.D0 1128 1129 1130Cfah d2F_cbb(1,1) = 0.D0 (not needed) 1131Cfah d2F_cbb(1,2) = 0.D0 (not needed) 1132Cfah d2F_cbb(1,3) = 0.D0 1133Cfah d2F_cbb(1,4) = 0.D0 1134Cfah d2F_cbb(2,2) = !=0, but not needed 1135Cfah d2F_cbb(2,3) = 0.D0 1136Cfah d2F_cbb(2,4) = c_nbb*n*u_cbb**(-1 + n)* 1137Cfah - de_cbb_by_drhob* 1138Cfah - du_cbb_by_dzb + 1139Cfah - c_nbb*(-1 + n)*n*e_cbb* 1140Cfah - u_cbb**(-2 + n)* 1141Cfah - du_cbb_by_dzb* 1142Cfah - du_cbb_by_drhob + 1143Cfah - c_nbb*n*e_cbb*u_cbb**(-1 + n)* 1144Cfah - du_cbb_by_drhob_dzb 1145Cfah d2F_cbb(3,3) = 0.D0 1146Cfah d2F_cbb(3,4) = 0.D0 1147Cfah d2F_cbb(4,4) = c_nbb*n*e_cbb*u_cbb**(-2 + n)* 1148Cfah - ((-1 + n)*du_cbb_by_dzb** 1149Cfah - 2 + u_cbb* 1150Cfah - du_cbb_by_dzb_dzb) 1151 1152Cfah d2F_cab(1,1) = not needed 1153Cfah d2F_cab(1,2) = not needed 1154Cfah d2F_cab(1,3) = c_nab*n*u_cab**(-2 + n)* 1155Cfah - ((-1 + n)*e_cab* 1156Cfah - du_cab_by_dza 1157Cfah - *du_cab_by_drhoa + 1158Cfah - u_cab*(de_cab_by_drhoa* 1159Cfah - du_cab_by_dza + e_cab* 1160Cfah - du_cab_by_drhoa_dza 1161Cfah - )) 1162Cfah d2F_cab(1,4) = c_nab*n*u_cab**(-2 + n)* 1163Cfah - ( (-1 + n)*e_cab* 1164Cfah - du_cab_by_dzb 1165Cfah - *du_cab_by_drhoa + 1166Cfah - u_cab* 1167Cfah - (de_cab_by_drhoa* 1168Cfah - du_cab_by_dzb + 1169Cfah - e_cab* 1170Cfah - du_cab_by_drhoa_dzb)) 1171Cfah d2F_cab(2,2) = not needed 1172Cfah d2F_cab(2,3) = c_nab*n*u_cab**(-2 + n)* 1173Cfah - ((-1 + n)*e_cab* 1174Cfah - du_cab_by_dza 1175Cfah - *du_cab_by_drhob + 1176Cfah - u_cab* 1177Cfah - (de_cab_by_drhob* 1178Cfah - du_cab_by_dza + 1179Cfah - e_cab* 1180Cfah - du_cab_by_drhob_dza)) 1181Cfah d2F_cab(2,4) = c_nab*n*u_cab**(-2 + n)* 1182Cfah - ((-1 + n)*e_cab* 1183Cfah - du_cab_by_dzb 1184Cfah - *du_cab_by_drhob + 1185Cfah - u_cab*(de_cab_by_drhob* 1186Cfah - du_cab_by_dzb + e_cab* 1187Cfah - du_cab_by_drhob_dzb )) 1188Cfah d2F_cab(3,3) = c_nab*n*e_cab* 1189Cfah - u_cab**(-2 + n)* 1190Cfah - ((-1 + n)*du_cab_by_dza**2 + 1191Cfah - u_cab* 1192Cfah - du_cab_by_dza_dza) 1193Cfah d2F_cab(3,4) = c_nab*n*e_cab* 1194Cfah - u_cab**(-2 + n)* 1195Cfah - ((-1 + n)*du_cab_by_dzb* 1196Cfah - du_cab_by_dza 1197Cfah - + u_cab* 1198Cfah - du_cab_by_dza_dzb) 1199Cfah d2F_cab(4,4) = c_nab*n*e_cab* 1200Cfah - u_cab**(-2 + n)* 1201Cfah - ((-1 + n)*du_cab_by_dzb**2 + 1202Cfah - u_cab* 1203Cfah - du_cab_by_dzb_dzb) 1204Cfah 1205Cfah here, the second derivatives are completed (Schwartz's rule: 1206Cfah df/(dadb) = df/(dbda) 1207Cfah DO i = 1, 4 1208Cfah DO j = i, 4 1209Cfah d2F_xa(j,i) = d2F_xa(i,j) 1210Cfah d2F_xb(j,i) = d2F_xb(i,j) 1211Cfah d2F_caa(j,i) = d2F_caa(i,j) 1212Cfah d2F_cbb(j,i) = d2F_cbb(i,j) 1213Cfah d2F_cab(j,i) = d2F_cab(i,j) 1214Cfah ENDDO 1215Cfah ENDDO 1216 1217Cfah test for zero densities (as in beta part of H atom): 1218 IF (rhob.LT.tol_rho) THEN 1219 DO i = 1, 4 1220 dF_xb(i) = 0.D0 1221 dF_cbb(i) = 0.D0 1222 dF_cab(i) = 0.D0 1223Cfah DO j = 1, 4 1224Cfah d2F_xb(i,j) = 0.D0 1225Cfah d2F_cbb(i,j) = 0.D0 1226Cfah d2F_cab(i,j) = 0.D0 1227Cfah ENDDO 1228 ENDDO 1229 ENDIF 1230 1231 IF (rhoa.LT.tol_rho) THEN 1232 DO i = 1, 4 1233 dF_xa(i) = 0.D0 1234 dF_caa(i) = 0.D0 1235 dF_cab(i) = 0.D0 1236Cfah DO j = 1, 4 1237Cfah d2F_xa(i,j) = 0.D0 1238Cfah d2F_caa(i,j) = 0.D0 1239Cfah d2F_cab(i,j) = 0.D0 1240Cfah ENDDO 1241 ENDDO 1242 ENDIF 1243 1244 1245Cfah Sum up all the partial derivatives with respect to the same function 1246Cfah of terms containing different powers of u with the help of the big outer 1247Cfah loop: 1248 1249Cfah have the partial derivative 1250 1251 DO i = 1, 4 1252 F(n*3+1,i) = dF_xa(i) + dF_xb(i) 1253 F(n*3+2,i) = dF_caa(i) +dF_cbb(i) 1254 F(n*3+3,i) = dF_cab(i) 1255Cfah DO j = 1, 4 1256Cfah FF(n*3+1,i,j) = d2F_xa(i,j) + d2F_xb(i,j) 1257Cfah FF(n*3+2,i,j) = d2F_caa(i,j) + d2F_cbb(i,j) 1258Cfah FF(n*3+3,i,j) = d2F_cab(i,j) 1259Cfah ENDDO 1260 ENDDO 1261 1262Cfah these partial derivatives have not been computed because they are 1263Cfah zero since we don't have a gradrhoagradrhob term in the Becke V functional 1264C F(n*3+1,5) = 0 1265C F(n*3+2,5) = 0 1266C F(n*3+3,5) = 0 1267Cfah DO i = 1, 5 1268Cfah FF(n*3+1,i,5) = 0 1269Cfah FF(n*3+2,i,5) = 0 1270Cfah FF(n*3+3,i,5) = 0 1271Cfah FF(n*3+1,5,i) = 0 1272Cfah FF(n*3+2,5,i) = 0 1273Cfah FF(n*3+3,5,i) = 0 1274Cfah ENDDO 1275 1276 ENDDO 1277 1278 RETURN 1279 1280 END 1281 1282Cfah----------------------------------------------------------- 1283 1284Cfah----------------------------------------------------------- 1285 1286ccc DFT_XCKernel_PWLDA(ra, rb, FCLDA); 1287 subroutine DFT_XCKernel_PWLDA(ra, rb, FCLDA) 1288 implicit none 1289#include "2ndDerivB97.h" 1290 double precision ra, rb 1291 double precision FCLDA(0:_FXC_RARB) 1292c 1293 double precision Amat(2) 1294 double precision Amat2(3) 1295 double precision rho(3) 1296 double precision ec,qwght 1297c 1298 double precision func 1299c 1300 rho(1)=ra+rb 1301 rho(2)=ra 1302 rho(3)=rb 1303cinitialize 1304 ec = 0d0 1305 qwght=1d0 1306 Amat(1)=0d0 1307 Amat(2)=0d0 1308 Amat2(1)=0d0 1309 Amat2(2)=0d0 1310 Amat2(3)=0d0 1311 1312 call xc_pw91lda_d2(1d-20, 1d0, .true., .false., rho, 1313 & Amat, Amat2, 1, 2, Ec, qwght, .false., func) 1314 1315 FCLDA(_FXC_E)=ec 1316 FCLDA(_FXC_RA)=Amat(1) 1317 FCLDA(_FXC_RB)=Amat(2) 1318c 1319 FCLDA(_FXC_RARA)=Amat2(D2_RA_RA) 1320 FCLDA(_FXC_RBRB)=Amat2(D2_RB_RB) 1321 FCLDA(_FXC_RARB)=Amat2(D2_RA_RB) 1322c rubbish 1323 return 1324 end 1325 subroutine xc_fxc2acmat(nq,i,ipol,FX, 1326 A Amat,Cmat,Amat2,Cmat2) 1327 implicit none 1328#include "2ndDerivB97.h" 1329#include "mafdecls.fh" 1330 integer nq,ipol 1331 integer i 1332 double precision Amat(nq,ipol), Cmat(nq,*) 1333 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 1334 double precision FX(0:_FXC_NUMDERI) 1335 Amat(i,1) = Amat(i,1) + FX(_FXC_RA) 1336 Cmat(i,D1_GAA) = Cmat(i,D1_GAA) + FX(_FXC_GAA) 1337 Amat2(i,D2_RA_RA) = Amat2(i,D2_RA_RA) + FX(_FXC_RARA) 1338 Amat2(i,D2_RA_RB) = Amat2(i,D2_RA_RB) + FX(_FXC_RARB) 1339 Cmat2(i,D2_GAA_GAA) = Cmat2(i,D2_GAA_GAA) + 1340 + FX(_FXC_GAAGAA) 1341 Cmat2(i,D2_GAA_GBB) = Cmat2(i,D2_GAA_GBB) + 1342 + FX(_FXC_GAAGBB) 1343 Cmat2(i,D2_RA_GAA) = Cmat2(i,D2_RA_GAA) + 1344 + FX(_FXC_RAGAA) 1345 Cmat2(i,D2_RA_GBB) = Cmat2(i,D2_RA_GBB) + 1346 + FX(_FXC_RAGBB) 1347 Cmat2(i,D2_RB_GAA) = Cmat2(i,D2_RB_GAA) + 1348 + FX(_FXC_RBGAA) 1349 if(ipol.eq.2) then 1350 Amat(i,2) = Amat(i,2) + FX(_FXC_RB) 1351 Cmat(i,D1_GBB) = Cmat(i,D1_GBB) + FX(_FXC_GBB) 1352 Amat2(i,D2_RB_RB) = Amat2(i,D2_RB_RB) + FX(_FXC_RBRB) 1353 Cmat2(i,D2_GBB_GBB) = Cmat2(i,D2_GBB_GBB) + 1354 + FX(_FXC_GBBGBB) 1355 Cmat2(i,D2_RB_GBB) = Cmat2(i,D2_RB_GBB) + 1356 + FX(_FXC_RBGBB) 1357 endif 1358 return 1359 end 1360 subroutine xc_htch_loadcf(functional,sol,max_pow_u) 1361 implicit none 1362#include "errquit.fh" 1363 character*4 functional ! functional name [input] 1364 integer max_pow_u ! functional name [input] 1365#include "xc_hcth.fh" 1366 double precision sol((limpow+1)*3) ! [out] 1367c 1368 integer nofunc 1369c 1370 nofunc = -1 ! take care of compiler warnings 1371 do n=1,numfunc 1372 if(functional.eq.funcnam(n)) nofunc=n 1373 enddo 1374 if(nofunc.eq.-1) call errquit('xchcth: cant pair funcname ',0, 1375 & UNKNOWN_ERR) 1376 max_pow_u=maxpow(nofunc) 1377 do n=1,3*(limpow+1) 1378 sol(n)=coeffs(n,nofunc) 1379 enddo 1380 1381C please refer to these coeffs as THCH1/iterate-e750-g500-v1-m4-n4 1382 1383c sol( 1) = 0.109320D+01 1384c sol( 2) = 0.222601D+00 1385c sol( 3) = 0.729974D+00 1386c sol( 4) = -0.744056D+00 1387c sol( 5) = -0.338622D-01 1388c sol( 6) = 0.335287D+01 1389c sol( 7) = 0.559920D+01 1390c sol( 8) = -0.125170D-01 1391c sol( 9) = -0.115430D+02 1392c sol(10) = -0.678549D+01 1393c sol(11) = -0.802496D+00 1394c sol(12) = 0.808564D+01 1395c sol(13) = 0.449357D+01 1396c sol(14) = 0.155396D+01 1397c sol(15) = -0.447857D+01 1398 return 1399 end 1400