1c 2c Wrapper routine for evaluating XC functional and derivatives, 3c and combining them with quadrature weights. 4c Copied from xc_quadv0_a. 5c 6c BGJ (8/98) 7c 8C> \ingroup nwdft_xc 9C> @{ 10C> 11C> \brief the functional evaluation routine 12C> 13C> This routine iterates over the components of the current density 14C> functional and evaluates and accumulates all density dependent 15C> terms. 16C> At present it will evaluate 17C> 18C> - the exchange-correlation energy 19C> 20C> - the exchange-correlation 1st order partial derivatives 21C> 22C> - the exchange-correlation 2nd order partial derivatives 23C> 24C> as desired. 25C> 26 Subroutine xc_eval_fnl(rho, delrho, Amat, Amat2, Cmat, Cmat2, 27 & nq, Ex, Ec, qwght, grad, ldew, func, 28 & do_2nd, ttau, kske, Mmat, Mmat2, 29 & laprho, kslap, Lmat, 30 & StericEnergy, 31 & do_3rd, Amat3, Cmat3, ldew2) 32c 33 implicit none 34c 35#include "cdft.fh" 36c xc-second derivative header file 37#include "dft2drv.fh" 38c xc-third derivative header file 39#include "dft3drv.fh" 40#include "stdio.fh" 41#include "steric.fh" 42cHvD 43#include "errquit.fh" 44#include "mafdecls.fh" 45#include "util.fh" 46 integer l_f, k_f 47 integer l_g, k_g 48 integer ii,iq 49cHvD 50c 51c List of functionals with implemented third derivatives. 52c 53c Exchange: Slater/Dirac 54c PBE 96 (original, RPBE, and revPBE) 55c Becke 88 56c CAM-Becke 88 57c CAM-PBE 96 58c BNL 59c LC-wPBE 60c Correlation: VWN (1-5 and 1 RPA) 61c LYP 62c PBE 96 63c PW91 LDA (needed for the default local part of PBE) 64c Perdew 86 (restricted only at the moment) 65c Perdew 81 LDA (needed for the default local part of Perdew 86) 66c 67 integer nq, is12x 68 double precision rho(*) !< [Input] The electron density 69 !< \f$\rho^\alpha\f$ and 70 !< \f$\rho^\beta\f$. 71 double precision delrho(*) !< [Input] The electron density gradient 72 !< \f$\partial\rho^\alpha/\partial r\f$ 73 !< and 74 !< \f$\partial\rho^\beta/\partial r\f$ 75 double precision laprho(*) 76c 77 double precision Amat(*), Cmat(*) 78 double precision Amat2(*),Cmat2(*) 79 double precision Amat3(*),Cmat3(*) 80 double precision ttau(*) !< [Input] The kinetic energy density 81 double precision Mmat(*), Mmat2(*), Mmat3(2) 82 double precision Lmat(*) 83c 84 double precision Ex !< [Output] The exchange energy 85 double precision Ec !< [Output] The correlation energy 86 double precision StericEnergy 87 double precision qwght(nq), func(nq) 88 logical grad, kske, dohcth, use_nwxc, kslap 89c 90 character*4 whichf 91 92 logical do_2nd, do_3rd 93 logical ldew, ldew2, ldew3 94cc AJL/Unused 95c logical out1 96cc AJL/End 97c 98 double precision eps,dumd 99 integer nx,nc,dumi 100 parameter (eps=1.e-8) 101c 102c Initialize the XC potential and energy sampling matrices. 103c 104 call dcopy(ipol*nq, 0.0d0, 0, Amat, 1) 105 if(grad) call dcopy(3*nq*ipol, 0.0d0, 0, Cmat, 1) 106 if(kske) call dcopy(nq*ipol, 0.0d0, 0, Mmat, 1) 107 if(kslap) call dcopy(nq*ipol, 0.0d0, 0, Lmat, 1) 108 if (do_2nd) then 109 call dcopy(nq*NCOL_AMAT2, 0.0d0, 0, Amat2, 1) 110 if (grad) call dcopy(nq*NCOL_CMAT2, 0.0d0, 0, Cmat2, 1) 111 if (kske) call dcopy(nq*NCOL_MMAT2, 0.0d0, 0, Mmat2, 1) 112 endif 113c 114c Initialize the 3rd derivatives. Note that we may need 115c the 2nd derivative matrices in some cases as well, so 116c those are initialized also. 117 if (do_3rd) then 118 call dcopy(nq*NCOL_AMAT2, 0.0d0, 0, Amat2, 1) 119 call dcopy(nq*NCOL_AMAT3, 0.0d0, 0, Amat3, 1) 120 if (grad) call dcopy(nq*NCOL_CMAT2, 0.0d0, 0, Cmat2, 1) 121 if (grad) call dcopy(nq*NCOL_CMAT3, 0.0d0, 0, Cmat3, 1) 122 if (kske) call dcopy(nq*NCOL_MMAT2, 0.0d0, 0, Mmat2, 1) 123 if (kske) call errquit("xc_eval_fnc: do_3rd and kske",0, 124 + CAPMIS_ERR) 125 endif 126c 127c This prevents the XC-kernel from being multiplied 128c with the quadrature weights in TDDFT gradients. 129c 130 ldew3 = (do_3rd.and.ldew) 131 if (ldew.or.ldew2.or.ldew3) call dcopy(nq, 0.d0, 0, func, 1) 132c 133 use_nwxc = util_module_avail("nwxc") 134 if (use_nwxc) then 135 call nwxc_getvals("nwxc_is_on",use_nwxc) 136 endif 137 if (use_nwxc) then 138c 139c Use the nwxc functional library 140c 141 if (.not.ma_push_get(mt_dbl,nq,"ffunc",l_f,k_f)) 142 + call errquit("xc_eval_fnc: no memory",nq,MA_ERR) 143 if (.not.ma_push_get(mt_dbl,3*nq,"rgamma",l_g,k_g)) 144 + call errquit("xc_eval_fnc: no memory",3*nq,MA_ERR) 145 if (grad) call calc_rgamma(ipol,nq,delrho,dbl_mb(k_g)) 146 if (ipol.eq.1) then 147 if (do_3rd) then 148 call nwxc_eval_df3(ipol,nq,rho(1),dbl_mb(k_g),ttau, 149 + dbl_mb(k_f), 150 + Amat,Amat2,Amat3,Cmat,Cmat2,Cmat3, 151 + Mmat,Mmat2,Mmat3) 152 else if (do_2nd) then 153 call nwxc_eval_df2(ipol,nq,rho(1),dbl_mb(k_g),ttau, 154 + dbl_mb(k_f), 155 + Amat,Amat2,Cmat,Cmat2,Mmat,Mmat2) 156 else 157 call nwxc_eval_df(ipol,nq,rho(1),dbl_mb(k_g),ttau, 158 + dbl_mb(k_f), 159 + Amat,Cmat,Mmat) 160 endif 161 else 162 if (do_3rd) then 163 call nwxc_eval_df3(ipol,nq,rho(nq+1),dbl_mb(k_g),ttau, 164 + dbl_mb(k_f), 165 + Amat,Amat2,Amat3,Cmat,Cmat2,Cmat3, 166 + Mmat,Mmat2,Mmat3) 167 else if (do_2nd) then 168 call nwxc_eval_df2(ipol,nq,rho(nq+1),dbl_mb(k_g),ttau, 169 + dbl_mb(k_f), 170 + Amat,Amat2,Cmat,Cmat2,Mmat,Mmat2) 171 else 172 call nwxc_eval_df(ipol,nq,rho(nq+1),dbl_mb(k_g),ttau, 173 + dbl_mb(k_f), 174 + Amat,Cmat,Mmat) 175 endif 176 endif 177c 178c The NWXC module calculates derivatives wrt the proper 179c kinetic energy density. It seems that the original subroutines 180c calculate derivatives wrt. twice the kinetic energy density 181c (in many cases people choose to absorb the factor half in the 182c kinetic energy density in the functional expression). As a 183c result there is factor half between what NWXC calculates and 184c what NWChem expects. Hence we need to scale Mmat here. 185c 186 if (kske) then 187 call dscal(nq*ipol,0.5d0,mmat,1) 188 endif 189c 190 Ex = Ex + ddot(nq,dbl_mb(k_f),1,qwght,1) 191 if (ldew) then 192 call dcopy(nq,dbl_mb(k_f),1,func,1) 193 endif 194 if (.not.ma_pop_stack(l_g)) 195 + call errquit("xc_eval_fnc: no free",3*nq,MA_ERR) 196 if (.not.ma_pop_stack(l_f)) 197 + call errquit("xc_eval_fnc: no free",nq,MA_ERR) 198c 199c Combine with quadrature weights 200c 201c Daniel (1-11-13): Added XC-third derivative stuff. This currently 202c doesn't include any functionality for meta-GGAs. 203c Daniel (1-17-13): Because we need derivatives of the quadrature 204c weights for TDDFT gradients, we don't want to multiply the A 205c and C matrices with the quadrature weights here (which would be 206c double counting). This prevents us from multiplying the 207c quadrature weights into the matrices and then dividing them out 208c later, which is prone to numerical problems (and slow). 209c Daniel (2-4-13): ldew3 is for the dfxc*(X+Y)*(X+Y) term in the TDDFT 210c gradients. ldew2 is for the dVxc*P term in the TDDFT gradients. 211 if (.not.(ldew2.or.ldew3)) then 212 213 if (.not. do_2nd .and. .not. do_3rd) then 214 call setACmat(delrho, Amat, Cmat, qwght, ipol, nq, grad, 215 & (.not. do_2nd), kske, Mmat, kslap, Lmat) 216 else if (.not. do_3rd) then 217 call setACmat_d2(delrho, Amat, Amat2, Cmat, Cmat2, qwght, 218 & ipol, nq, grad, (.not. do_2nd), kske, Mmat, Mmat2, 219 & .false.) 220 else 221 call setACmat_d3(delrho, Amat, Amat2, Amat3, Cmat, Cmat2, 222 & Cmat3, qwght, ipol, nq, grad, (.not. do_3rd), 223 & .false.,.false.) 224 endif 225 endif 226 return 227 endif 228c if (ldew.or.ldew3) call dfill(nq, 0.d0, func, 1) 229c if (ldew) call dfill(nq, 0.d0, func, 1) 230c 231c warning!! xc_dirac has to be called before all the other 232c XC routines 233c 234 if (abs(xfac(2)).gt.eps)then 235 if (.not. do_2nd .and. .not. do_3rd) then 236 call xc_dirac(tol_rho, xfac(2), lxfac(2), nlxfac(2), rho, 237 & Amat, nq, ipol, Ex, qwght, 238 & ldew, func) 239 else if (.not. do_3rd) then 240 call xc_dirac_d2(tol_rho, xfac(2), lxfac(2), nlxfac(2), rho, 241 & Amat, Amat2, nq, ipol, Ex, qwght, 242 & .false., func) 243 else 244 call xc_dirac_d3(tol_rho, xfac(2), lxfac(2), nlxfac(2), rho, 245 & Amat, Amat2, Amat3, nq, ipol, Ex, qwght, 246 & .false., func) 247 endif 248 endif 249c 250c LB94 potential and energy (not variational!) 251c 252 if (abs(xfac(70)).gt.eps)then 253 if (.not. do_2nd) then 254 call xc_lb94_e(tol_rho, xfac(70), rho, delrho, Amat, nq, 255 I ipol, Ex, qwght,ldew,func,ttau) 256 else 257 call xc_lb94_e_d2(tol_rho, xfac(70), rho, delrho, Amat, 258 I Amat2, nq, ipol, Ex, qwght,.false.,func,ttau) 259c call errquit( ' lb94 2nds not ready yet',0,0) 260 endif 261 endif 262c 263 if (abs(xfac(3)).gt.eps)then 264 if (.not. do_2nd .and. .not. do_3rd) then 265 call xc_becke88(tol_rho, xfac(3), lxfac(3), nlxfac(3), 266 & rho, delrho, Amat, Cmat, nq, ipol, 267 & Ex, qwght,ldew,func) 268 else if (.not. do_3rd) then 269 call xc_becke88_d2(tol_rho, xfac(3), lxfac(3), nlxfac(3), 270 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 271 & Ex, qwght,ldew,func) 272 else 273 call xc_becke88_d3(tol_rho, xfac(3), lxfac(3), nlxfac(3), 274 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 275 & nq, ipol, Ex, qwght,ldew,func) 276 endif 277 endif 278c 279 if (abs(xfac(16)).gt.eps)then 280 call xc_optx(tol_rho, xfac(16), 2, 281 & rho, delrho, Amat, Cmat, nq, ipol, 282 & Ex, qwght,ldew,func) 283 endif 284C 285C hcth , becke97s functionals 286C 287 nx = 4 ! take care of compiler warnings 288 nc = 13 289c 290 dohcth=.false. 291 if (abs(xfac(4)).gt.eps.or.abs(cfac(13)).gt.eps)then 292 whichf='hcth' 293 dohcth=.true. 294 nx=4 295 nc=13 296 elseif (abs(xfac(10)).gt.eps.or.abs(cfac(16)).gt.eps)then 297 whichf='h120' 298 dohcth=.true. 299 nx=10 300 nc=16 301 elseif (abs(xfac(11)).gt.eps.or.abs(cfac(17)).gt.eps)then 302 whichf='h147' 303 dohcth=.true. 304 nx=11 305 nc=17 306 elseif (abs(xfac(5)).gt.eps.or.abs(cfac(14)).gt.eps)then 307 whichf='b970' 308 dohcth=.true. 309 nx=5 310 nc=14 311 elseif (abs(xfac(6)).gt.eps.or.abs(cfac(15)).gt.eps)then 312 whichf='b971' 313 dohcth=.true. 314 nx=6 315 nc=15 316 elseif (abs(xfac(12)).gt.eps.or.abs(cfac(18)).gt.eps)then 317 whichf='b980' 318 dohcth=.true. 319 nx=12 320 nc=18 321 elseif (abs(xfac(13)).gt.eps.or.abs(cfac(19)).gt.eps)then 322 whichf='b97g' 323 dohcth=.true. 324 nx=13 325 nc=19 326 elseif (abs(xfac(14)).gt.eps.or.abs(cfac(20)).gt.eps)then 327 whichf='h407' 328 dohcth=.true. 329 nx=14 330 nc=20 331 elseif (abs(xfac(15)).gt.eps.or.abs(cfac(21)).gt.eps)then 332 whichf='hp14' 333 dohcth=.true. 334 nx=15 335 nc=21 336 elseif (abs(xfac(17)).gt.eps.or.abs(cfac(23)).gt.eps)then 337 whichf='b972' 338 dohcth=.true. 339 nx=17 340 nc=23 341 elseif (abs(xfac(20)).gt.eps.or.abs(cfac(26)).gt.eps)then 342 whichf='407p' 343 dohcth=.true. 344 nx=20 345 nc=26 346 elseif (abs(xfac(22)).gt.eps.or.abs(cfac(28)).gt.eps)then 347 whichf='b973' 348 dohcth=.true. 349 nx=22 350 nc=28 351 elseif (abs(xfac(39)).gt.eps.or.abs(cfac(41)).gt.eps)then 352 whichf='b97d' 353 dohcth=.true. 354 nx=39 355 nc=41 356 elseif (abs(cfac(45)).gt.eps)then 357 whichf='n120' 358 dohcth=.true. 359 nx=45 360 nc=45 361 endif 362 if(dohcth) then 363 if (.not. do_2nd) then 364 call xc_hcth(tol_rho, xfac(nx), lxfac(nx), nlxfac(nx), 365 , cfac(nc), lcfac(nc), nlxfac(nc), rho, 366 & delrho, Amat, Cmat, nq, ipol, Ex, Ec, qwght, 367 & ldew, func,whichf) 368 else 369 call xc_hcth_d2(tol_rho, xfac(nx), lxfac(nx), nlxfac(nx), 370 , cfac(nc), lcfac(nc), nlxfac(nc), rho, 371 & delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, Ec, 372 & qwght, ldew, func,whichf) 373 endif ! do_2nd 374 endif ! dohctch 375c 376c compute partial derivatives of the correlation energy functional. 377c 378 if (abs(cfac(1)).gt.eps)then 379 if (.not. do_2nd .and..not. do_3rd) then 380 call xc_vwn_5(tol_rho, cfac(1), rho, 381 & Amat, nq, ipol, Ec, qwght, 382 & ldew, func) 383 else if (.not. do_3rd) then 384 call xc_vwn_5_d2(tol_rho, cfac(1), rho, 385 & Amat, Amat2, nq, ipol, Ec, qwght, 386 & ldew, func) 387 else 388 call xc_vwn_5_d3(tol_rho, cfac(1), rho, 389 & Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 390 & ldew, func) 391 endif 392 endif 393c 394 if (abs(cfac(7)).gt.eps)then 395 if (.not. do_2nd .and..not. do_3rd) then 396 call xc_vwn_1_rpa(tol_rho, cfac(7), 397 & rho, Amat, nq, ipol, Ec, 398 & qwght, ldew, func) 399 else if (.not. do_3rd) then 400 call xc_vwn_1_rpa_d2(tol_rho, cfac(7), 401 & rho, Amat, Amat2, nq, ipol, Ec, 402 & qwght, ldew, func) 403 else 404 call xc_vwn_1_rpa_d3(tol_rho, cfac(7), 405 & rho, Amat, Amat2, Amat3, nq, ipol, Ec, 406 & qwght, ldew, func) 407 endif 408 endif 409c 410 if (abs(cfac(8)).gt.eps)then 411 if (.not. do_2nd .and..not. do_3rd) then 412 call xc_vwn_1(tol_rho, cfac(8), rho, 413 & Amat, nq, ipol, Ec, qwght, 414 & ldew, func) 415 else if (.not. do_3rd) then 416 call xc_vwn_1_d2(tol_rho, cfac(8), rho, 417 & Amat, Amat2, nq, ipol, Ec, qwght, 418 & ldew, func) 419 else 420 call xc_vwn_1_d3(tol_rho, cfac(8), rho, 421 & Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 422 & ldew, func) 423 endif 424 endif 425c 426 if (abs(cfac(9)).gt.eps)then 427 if (.not. do_2nd .and..not. do_3rd) then 428 call xc_vwn_2(tol_rho, cfac(9), rho, 429 & Amat, nq, ipol, Ec, qwght, 430 & ldew, func) 431 else if (.not. do_3rd) then 432 call xc_vwn_2_d2(tol_rho, cfac(9), rho, 433 & Amat, Amat2, nq, ipol, Ec, qwght, 434 & ldew, func) 435 else 436 call xc_vwn_2_d3(tol_rho, cfac(9), rho, 437 & Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 438 & ldew, func) 439 endif 440 endif 441c 442 if (abs(cfac(10)).gt.eps)then 443 if (.not. do_2nd .and..not. do_3rd) then 444 call xc_vwn_3(tol_rho, cfac(10), 445 & rho, Amat, nq, ipol, Ec, qwght, 446 & ldew, func) 447 else if (.not. do_3rd) then 448 call xc_vwn_3_d2(tol_rho, cfac(10), 449 & rho, Amat, Amat2, nq, ipol, Ec, qwght, 450 & ldew, func) 451 else 452 call xc_vwn_3_d3(tol_rho, cfac(10), 453 & rho, Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 454 & ldew, func) 455 endif 456 endif 457c 458 if (abs(cfac(11)).gt.eps)then 459 if (.not. do_2nd .and..not. do_3rd) then 460 call xc_vwn_4(tol_rho, cfac(11), 461 & rho, Amat, nq, ipol, Ec, qwght, 462 & ldew, func) 463 else if (.not. do_3rd) then 464 call xc_vwn_4_d2(tol_rho, cfac(11), 465 & rho, Amat, Amat2, nq, ipol, Ec, qwght, 466 & ldew, func) 467 else 468 call xc_vwn_4_d3(tol_rho, cfac(11), 469 & rho, Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 470 & ldew, func) 471 endif 472 endif 473c 474 if (abs(cfac(6)).gt.eps)then 475 if (.not. do_2nd .and..not. do_3rd) then 476 call xc_pw91lda(tol_rho, cfac(6), lcfac(6), nlcfac(6), 477 & rho, Amat, nq, ipol, Ec, 478 & qwght, ldew, func) 479 else if (.not. do_3rd) then 480 call xc_pw91lda_d2(tol_rho, cfac(6), lcfac(6), nlcfac(6), 481 & rho, Amat, Amat2, nq, ipol, Ec, qwght, 482 & ldew, func) 483 else 484 call xc_pw91lda_d3(tol_rho, cfac(6), lcfac(6), nlcfac(6), 485 & rho, Amat, Amat2, Amat3, nq, ipol, Ec, qwght, 486 & ldew, func) 487 endif 488 endif 489c 490 if (abs(cfac(2)).gt.eps)then 491 if (.not. do_2nd .and. .not. do_3rd) then 492 call xc_lyp88(tol_rho, cfac(2), 493 & rho, delrho, Amat, Cmat, nq, ipol, Ec, 494 & qwght, ldew, func) 495 else if (.not. do_3rd) then 496 call xc_lyp88_d2(tol_rho, cfac(2), 497 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, Ec, 498 & qwght, ldew, func) 499 else 500 call xc_lyp88_d3(tol_rho, cfac(2), 501 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 502 & nq, ipol, Ec, qwght, ldew, func) 503 endif 504 endif 505c 506 if (abs(cfac(3)).gt.eps)then 507 if (.not. do_2nd .and..not. do_3rd) then 508 call xc_p81(tol_rho, cfac(3), lcfac(3), nlcfac(3), rho, 509 & Amat, nq, ipol, Ec, qwght, ldew, func) 510 else if (.not. do_3rd) then 511 call xc_p81_d2(tol_rho, cfac(3), lcfac(3), nlcfac(3), rho, 512 & Amat, Amat2, nq, ipol, Ec, qwght, ldew, func) 513 else 514 call xc_p81_d3(tol_rho, cfac(3), lcfac(3), nlcfac(3), rho, 515 & Amat, Amat2, Amat3, nq, ipol, Ec, qwght, ldew, func) 516 endif 517 endif 518c 519 if (abs(cfac(4)).gt.eps)then 520 if (.not. do_2nd .and..not. do_3rd) then 521 call xc_perdew86(tol_rho, cfac(4), lcfac(4), nlcfac(4), 522 & rho, delrho, Amat, Cmat, nq, ipol, 523 & Ec, qwght, ldew, func) 524 else if (.not. do_3rd) then 525 call xc_perdew86_d2(tol_rho, cfac(4), lcfac(4), nlcfac(4), 526 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 527 & Ec, qwght, ldew, func) 528 else 529 call xc_perdew86_d3(tol_rho, cfac(4), lcfac(4), nlcfac(4), 530 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 531 & nq, ipol, Ec, qwght, ldew, func) 532 endif 533 endif 534c 535c PW91 is special in that the GGA part is dependent on 536c the E(LDA) ... so more info has to be passed in. 537c 538 if (abs(cfac(5)).gt.eps)then 539 if (.not. do_2nd) then 540 call xc_perdew91(tol_rho, cfac, lcfac, nlcfac, rho, 541 & delrho, Amat, Cmat, nq, ipol, 542 & Ec, qwght, ldew, func) 543 else 544 call xc_perdew91_d2(tol_rho, cfac, lcfac, nlcfac, rho, 545 & delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 546 & Ec, qwght, ldew, func) 547 endif 548 endif 549c 550c PBE96 is special in that the GGA part is dependent on 551c the E(LDA) ... so more info has to be passed in. 552c 553c same is true for revTPSS-variant of it 554c 555 if (abs(cfac(12)).gt.eps.or.abs(cfac(64)).gt.eps)then 556 if (.not. do_2nd .and..not. do_3rd) then 557 call xc_cpbe96(tol_rho, cfac, lcfac, nlcfac, rho, 558 & delrho, Amat, Cmat, nq, ipol, 559 & Ec, qwght, ldew, func) 560 else if (.not. do_3rd) then 561 call xc_cpbe96_d2(tol_rho, cfac, lcfac, nlcfac, rho, 562 & delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 563 & Ec, qwght, ldew, func) 564 else 565 call xc_cpbe96_d3(tol_rho, cfac, lcfac, nlcfac, rho, 566 & delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 567 & nq, ipol, Ec, qwght, ldew, func) 568 endif 569 endif 570 if (abs(xfac(7)).gt.eps)then 571 if (.not. do_2nd .and..not. do_3rd) then 572 call xc_xpbe96('orig', 573 T tol_rho, xfac(7), lxfac(7), nlxfac(7), 574 & rho, delrho, Amat, Cmat, nq, ipol, 575 & Ex, qwght,ldew,func) 576 else if (.not. do_3rd) then 577 call xc_xpbe96_d2('orig', 578 T tol_rho, xfac(7), lxfac(7), nlxfac(7), 579 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 580 & Ex, qwght,ldew,func) 581 else 582 call xc_xpbe96_d3('orig', 583 T tol_rho, xfac(7), lxfac(7), nlxfac(7), 584 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 585 & nq, ipol, Ex, qwght,ldew,func) 586 endif 587 endif 588 if (abs(xfac(30)).gt.eps)then 589 if (.not. do_2nd .and..not. do_3rd) then 590 call xc_xpbe96('rpbe', 591 T tol_rho, xfac(30), lxfac(30), nlxfac(30), 592 & rho, delrho, Amat, Cmat, nq, ipol, 593 & Ex, qwght,ldew,func) 594 else if (.not. do_3rd) then 595 call xc_xpbe96_d2('rpbe', 596 T tol_rho, xfac(30), lxfac(30), nlxfac(30), 597 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 598 & Ex, qwght,ldew,func) 599 else 600 call xc_xpbe96_d3('rpbe', 601 T tol_rho, xfac(30), lxfac(30), nlxfac(30), 602 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 603 & nq, ipol, Ex, qwght,ldew,func) 604 endif 605 endif 606 if (abs(xfac(31)).gt.eps)then 607 if (.not. do_2nd .and..not. do_3rd) then 608 call xc_xpbe96('revp', 609 T tol_rho, xfac(31), lxfac(31), nlxfac(31), 610 & rho, delrho, Amat, Cmat, nq, ipol, 611 & Ex, qwght,ldew,func) 612 else if (.not. do_3rd) then 613 call xc_xpbe96_d2('revp', 614 T tol_rho, xfac(31), lxfac(31), nlxfac(31), 615 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 616 & Ex, qwght,ldew,func) 617 else 618 call xc_xpbe96_d3('revp', 619 T tol_rho, xfac(31), lxfac(31), nlxfac(31), 620 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 621 & nq, ipol, Ex, qwght,ldew,func) 622 endif 623 endif 624c 625c SSB-D 626c consists of three parts, SSB-1 (correction to PBEx), 627c portion of KT1 gradient correction, and sPBEc 628c (it also includes a portion of Grimme's dispersion correction) 629c see: Swart, Sola, Bickelhaupt JCP 2009, 131, 094103 630c 631c sPBEc is special in that the GGA part is dependent on 632c the E(LDA) ... so more info has to be passed in. 633c 634 if (abs(cfac(46)).gt.eps)then 635 if (.not. do_2nd) then 636 call xc_spbe96(tol_rho, cfac, lcfac, nlcfac, 637 & rho, delrho, Amat, Cmat, nq, ipol, 638 & Ec, qwght, ldew, func) 639 else 640 call xc_spbe96_d2(tol_rho, cfac, lcfac, nlcfac, 641 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 642 & Ec, qwght, ldew, func) 643 endif 644 endif 645 if (abs(xfac(46)).gt.eps)then 646 if (.not. do_2nd) then 647c 648c first the part that depends on s (correction to PBEx) 649c 650 call xc_ssbD_1(tol_rho, xfac(46), lxfac(46), nlxfac(46), 651 & rho, delrho, Amat, Cmat, nq, ipol, 652 & Ex, qwght, ldew, func) 653 else 654c 655c first the part that depends on s (correction to PBEx) 656c 657 call xc_ssbD_1_d2(tol_rho, xfac(46), lxfac(46), nlxfac(46), 658 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 659 & Ex, qwght, ldew, func) 660 endif 661 endif 662c 663c kt1 664c 665 if (abs(xfac(47)).gt.eps)then 666 if (.not. do_2nd) then 667 call xc_kt1(tol_rho, xfac(47), rho, delrho, 668 & Amat, Cmat, nq, ipol, Ex, qwght, ldew, func) 669 670 else 671 call xc_kt1_d2(tol_rho, xfac(47), rho, delrho, 672 & Amat, Amat2, Cmat, Cmat2, nq, ipol, 673 & Ex, qwght, ldew, func) 674 endif 675 endif 676c 677c s12g 678c 679 if (abs(xfac(60)).gt.eps) then 680 is12x = 1 681 if (.not. do_2nd) then 682 call xc_s12x(tol_rho, xfac(60), lxfac(60), nlxfac(60), 683 & rho, delrho, Amat, Cmat, nq, ipol, 684 & Ex, qwght, ldew, func, is12x) 685 else 686 call xc_s12x_d2(tol_rho, xfac(60), lxfac(60), nlxfac(60), 687 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 688 & Ex, qwght, ldew, func, is12x) 689 endif 690 endif 691c 692c s12h 693c 694 if (abs(xfac(61)).gt.eps) then 695 is12x = 2 696 if (.not. do_2nd) then 697 call xc_s12x(tol_rho, xfac(61), lxfac(61), nlxfac(61), 698 & rho, delrho, Amat, Cmat, nq, ipol, 699 & Ex, qwght, ldew, func, is12x) 700 else 701 call xc_s12x_d2(tol_rho, xfac(61), lxfac(61), nlxfac(61), 702 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 703 & Ex, qwght, ldew, func, is12x) 704 endif 705 endif 706c 707c cam-s12g 708c 709 if (abs(xfac(62)).gt.eps) then 710 is12x = 1 711 if (.not. do_2nd) then 712 call xc_cams12x(tol_rho, xfac(62), lxfac(62), nlxfac(62), 713 & rho, delrho, Amat, Cmat, nq, ipol, 714 & Ex, qwght, ldew, func, is12x) 715 else 716 call xc_cams12x_d2(tol_rho, xfac(62), lxfac(62), nlxfac(62), 717 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 718 & Ex, qwght, ldew, func, is12x) 719 endif 720 endif 721c 722c cam-s12h 723c 724 if (abs(xfac(63)).gt.eps) then 725 is12x = 2 726 if (.not. do_2nd) then 727 call xc_cams12x(tol_rho, xfac(63), lxfac(63), nlxfac(63), 728 & rho, delrho, Amat, Cmat, nq, ipol, 729 & Ex, qwght, ldew, func, is12x) 730 else 731 call xc_cams12x_d2(tol_rho, xfac(63), lxfac(63), nlxfac(63), 732 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 733 & Ex, qwght, ldew, func, is12x) 734 endif 735 endif 736c 737 if (abs(xfac(8)).gt.eps)then 738 if (.not. do_2nd) then 739 call xc_gill96(tol_rho, xfac(8), lxfac(8), nlxfac(8), 740 & rho, delrho, Amat, Cmat, nq, ipol, 741 & Ex, qwght,ldew,func) 742 else 743 call xc_gill96_d2(tol_rho, xfac(8), lxfac(8), nlxfac(8), 744 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 745 & Ex, qwght,ldew,func) 746 endif 747 endif 748 if (abs(xfac(9)).gt.eps)then 749 if (.not. do_2nd) then 750 call xc_xpw91(tol_rho, xfac(9), lxfac(9), nlxfac(9), 751 & rho, delrho, Amat, Cmat, nq, ipol, 752 & Ex, qwght,ldew,func) 753 else 754 call xc_xpw91_d2(tol_rho, xfac(9), lxfac(9), nlxfac(9), 755 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 756 & Ex, qwght,ldew,func) 757 endif 758 endif 759 if (abs(cfac(22)).gt.eps)then 760 call xc_optc(rho, delrho, 761 & Amat, Cmat, nq, Ec, qwght,ldew,func) 762 endif 763 if (abs(xfac(19)).gt.eps)then 764 if (.not. do_2nd) then 765 call xc_xmpw91(tol_rho,xfac(19),lxfac(19),nlxfac(19), 766 & rho, delrho, Amat, Cmat, nq, ipol, 767 & Ex, qwght,ldew,func) 768 else 769 call xc_xmpw91_d2(tol_rho,xfac(19),lxfac(19),nlxfac(19), 770 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 771 & Ex, qwght,ldew,func) 772 endif 773 endif 774c 775 if (abs(xfac(25)).gt.eps.or.abs(cfac(24)).gt.eps)then 776 if (.not. do_2nd) then 777 call xc_ft97(tol_rho,xfac(25),lxfac(25),nlxfac(25), 778 . cfac(24),lcfac(24),nlcfac(24), 779 & rho, delrho, Amat, Cmat, nq, ipol, 780 & Ex, Ec, qwght,ldew,func) 781 else 782 call errquit('2nd derivative not available 783 &for this xc functional',0,0) 784! call xc_ft97_d2() 785 endif 786 endif 787c 788 if ((abs(cfac(36)).gt.eps).or.(abs(cfac(37)).gt.eps))then 789 if(abs(cfac(36)).gt.eps) then 790 nc=36 791 whichf='be88' 792 endif 793 if(abs(cfac(37)).gt.eps) then 794 nc=37 795 whichf='pb96' 796 endif 797 if (.not. do_2nd) then 798 call xc_op(tol_rho,whichf, 799 . cfac(nc),lcfac(nc),nlcfac(nc), 800 & rho, delrho, Amat, Cmat, nq, ipol, 801 & Ec, qwght,ldew,func) 802 else 803 call errquit('2nd derivative not available 804 &for this xc functional',0,0) 805! call xc_op_d2() 806 endif 807 endif 808c 809c meta GGA 810c 811 if (abs(xfac(18)).gt.eps)then 812 if (.not. do_2nd) then 813 call xc_xpkzb99(tol_rho, xfac(18), lxfac(18), nlxfac(18), 814 & rho, delrho, Amat, Cmat, nq, ipol, 815 & Ex, qwght,ldew,func,ttau, Mmat) 816 else 817 call xc_xpkzb99_d2() 818 endif 819 endif 820c 821c LB94 or CS00 correction is added to xc potential only 822c (xc functional and functional 2nds are unchanged) 823c 824 if (cs00) then 825 call xc_cs00(tol_rho, xfac(1), rho, delrho, Amat, nq, ipol, 826 & delta_ac, e_homo) 827 else if (lb94) then 828 call xc_lb94(tol_rho, xfac(1), rho, delrho, Amat, nq, ipol) 829 else if (ncap) then 830 if (delta_ac.eq.1.0d99) then 831 Amat(1:nq) = Amat(1:nq) - 832 & 0.053805222d0*(1d0 + dsqrt(1d0 - 2d0*e_homo/0.053805222d0)) 833 if (ipol.eq.2) then 834 Amat(nq+1:2*nq) = Amat(nq+1:2*nq) - 835 & 0.053805222d0*(1d0 + dsqrt(1d0 - 2d0*e_homo/0.053805222d0)) 836 endif 837 else 838 Amat(1:nq) = Amat(1:nq) - delta_ac 839 if (ipol.eq.2) Amat(nq+1:2*nq) = Amat(nq+1:2*nq) - delta_ac 840 endif 841 endif 842c 843c PKZB99-COR is special in that the GGA part is 844c defined to be PBE COR GGA and also is dependent on 845c the E(LDA) ... 846c the decision has been made to use the PW91-LDA as the 847c LDA-correlation. at present, this LDA cannot be 848c set by the user 849c 850 if (abs(cfac(25)).gt.eps)then 851 if (.not. do_2nd) then 852 call xc_cpkzb99(tol_rho, cfac(25), lcfac(25), nlcfac(25), 853 & rho, delrho, nq, ipol, 854 & Ec, qwght, ldew, func, ttau,Amat,Cmat,Mmat) 855 else 856 call xc_cpkzb99_d2() 857 endif 858 endif 859c 860c TPSS meta GGA 861c 862 if (abs(xfac(21)).gt.eps)then 863 if (.not. do_2nd) then 864 call xc_xtpss03(tol_rho, xfac(21), 865 & rho, delrho, Amat, Cmat, nq, ipol, 866 & Ex, qwght,ldew,func,ttau,Mmat) 867 else 868 call xc_xtpss03_d2() 869 endif 870 endif 871c 872c MVS meta GGA 873c 874 if (abs(xfac(64)).gt.eps)then 875 if (.not. do_2nd) then 876 call xc_xmvs15(tol_rho, xfac(64), 877 & rho, delrho, Amat, Cmat, nq, ipol, 878 & Ex, qwght,ldew,func,ttau,Mmat) 879 else 880 call xc_xmvs15_d2() 881 endif 882 endif 883c 884c TPSS03-COR is special in that the GGA part is 885c defined to be PBE COR GGA and also is dependent on 886c the E(LDA) ... 887c the decision has been made to use the PW91-LDA as the 888c LDA-correlation. at present, this LDA cannot be 889c set by the user 890c 891 if (abs(cfac(27)).gt.eps)then 892 if (.not. do_2nd) then 893 call xc_ctpss03(tol_rho, cfac(27), lcfac(27), nlcfac(27), 894 & rho, delrho, nq, ipol, 895 & Ec, qwght, ldew, func, ttau,Amat,Cmat,Mmat) 896 897 else 898 call xc_ctpss03_d2() 899 endif 900 endif 901c 902c Bc95-COR is special in that the GGA part is 903c defined to be dependent on 904c the E(LDA) ... 905c the decision has been made to use the PW91-LDA as the 906c LDA-correlation. at present, this LDA cannot be 907c set by the user 908c 909c Note that bc95, cpw6b95, cpwb6k use the same subroutine xc_bc95() 910c 911 if (abs(cfac(31)).gt.eps)then 912 if (.not. do_2nd) then 913 call xc_bc95(tol_rho, cfac(31), lcfac(31), nlcfac(31), 914 & rho, delrho, nq, ipol, 915 & Ec, qwght, ldew,func, ttau,Amat,Cmat,Mmat,0) 916 917 else 918 call xc_bc95_d2() 919 endif 920 endif 921c 922c PW6B95 Exchange part 923c 924 if (abs(xfac(26)).gt.eps)then 925 if (.not. do_2nd) then 926 call xc_xpw6(tol_rho,xfac(26),lxfac(26),nlxfac(26), 927 & rho, delrho, Amat, Cmat, nq, ipol, 928 & Ex, qwght,ldew,func,1) 929 else 930 call xc_xpw6_d2(tol_rho,xfac(26),lxfac(26),nlxfac(26), 931 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 932 & Ex, qwght,ldew,func,1) 933 endif 934 endif 935c 936c PWB6K Exchange part 937c 938 if (abs(xfac(27)).gt.eps)then 939 if (.not. do_2nd) then 940 call xc_xpw6(tol_rho,xfac(27),lxfac(27),nlxfac(27), 941 & rho, delrho, Amat, Cmat, nq, ipol, 942 & Ex, qwght,ldew,func,2) 943 else 944 call xc_xpw6_d2(tol_rho,xfac(27),lxfac(27),nlxfac(27), 945 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 946 & Ex, qwght,ldew,func,2) 947 endif 948 endif 949c 950c M05 meta GGA Exchange 951c 952 if (abs(xfac(28)).gt.eps)then 953 if (.not. do_2nd) then 954 call xc_xm05(tol_rho, xfac(28), lxfac(28), nlxfac(28), 955 & rho, delrho, Amat, Cmat, nq, ipol, 956 & Ex, qwght,ldew,func,ttau,Mmat,1) 957 else 958 call xc_xm05_d2() 959 endif 960 endif 961c 962c M05-2X meta GGA Exchange 963c 964 if (abs(xfac(29)).gt.eps)then 965 if (.not. do_2nd) then 966 call xc_xm05(tol_rho, xfac(29), lxfac(29), nlxfac(29), 967 & rho, delrho, Amat, Cmat, nq, ipol, 968 & Ex, qwght,ldew, func,ttau,Mmat,2) 969 else 970 call xc_xm05_d2() 971 endif 972 endif 973c 974c dlDF meta GGA Exchange 975c 976 if (abs(xfac(32)).gt.eps)then 977 if (.not. do_2nd) then 978 call xc_xdldf(tol_rho, xfac(32), lxfac(32), nlxfac(32), 979 & rho, delrho, Amat, Cmat, nq, ipol, 980 & Ex, qwght,ldew, func,ttau,Mmat) 981 else 982 call xc_xdldf_d2() 983 endif 984 endif 985c 986c VSXC meta GGA Exchange 987c 988 if (abs(xfac(33)).gt.eps)then 989 if (.not. do_2nd) then 990 call xc_xvs98(tol_rho, xfac(33), lxfac(33), nlxfac(33), 991 & rho, delrho, Amat, Cmat, nq, ipol, 992 & Ex, qwght,ldew, func,ttau,Mmat,1) 993 else 994 call xc_xvs98_d2() 995 endif 996 endif 997c 998c M06-L meta GGA Exchange 999c 1000 if (abs(xfac(34)).gt.eps)then 1001 if (.not. do_2nd) then 1002 call xc_xm06(tol_rho, xfac(34), lxfac(34), nlxfac(34), 1003 & rho, delrho, Amat, Cmat, nq, ipol, 1004 & Ex, qwght,ldew, func,ttau,Mmat,1) 1005 else 1006 call xc_xm06_d2() 1007 endif 1008 endif 1009c 1010c revM06 meta GGA Exchange 1011c 1012 if (abs(xfac(68)).gt.eps)then 1013 if (.not. do_2nd) then 1014 call xc_xm06(tol_rho, xfac(68), lxfac(68), nlxfac(68), 1015 & rho, delrho, Amat, Cmat, nq, ipol, 1016 & Ex, qwght,ldew, func,ttau,Mmat,6) 1017 else 1018 call xc_xm06_d2() 1019 endif 1020 endif 1021c 1022c revM06-L meta GGA Exchange 1023c 1024 if (abs(xfac(69)).gt.eps)then 1025 if (.not. do_2nd) then 1026 call xc_xm06(tol_rho, xfac(69), lxfac(69), nlxfac(69), 1027 & rho, delrho, Amat, Cmat, nq, ipol, 1028 & Ex, qwght,ldew, func,ttau,Mmat,5) 1029 else 1030 call xc_xm06_d2() 1031 endif 1032 endif 1033c 1034c M06-HF meta GGA Exchange 1035c 1036 if (abs(xfac(35)).gt.eps)then 1037 if (.not. do_2nd) then 1038 call xc_xm06(tol_rho, xfac(35), lxfac(35), nlxfac(35), 1039 & rho, delrho, Amat, Cmat, nq, ipol, 1040 & Ex, qwght,ldew, func,ttau,Mmat,2) 1041 else 1042 call xc_xm06_d2() 1043 endif 1044 endif 1045c 1046c M06 meta GGA Exchange 1047c 1048 if (abs(xfac(36)).gt.eps)then 1049 if (.not. do_2nd) then 1050 call xc_xm06(tol_rho, xfac(36), lxfac(36), nlxfac(36), 1051 & rho, delrho, Amat, Cmat, nq, ipol, 1052 & Ex, qwght,ldew, func,ttau,Mmat,3) 1053 else 1054 call xc_xm06_d2() 1055 endif 1056 endif 1057c 1058c M06-2X meta GGA Exchange 1059c 1060 if (abs(xfac(37)).gt.eps)then 1061 if (.not. do_2nd) then 1062 call xc_xm06(tol_rho, xfac(37), lxfac(37), nlxfac(37), 1063 & rho, delrho, Amat, Cmat, nq, ipol, 1064 & Ex, qwght,ldew, func,ttau,Mmat,4) 1065 else 1066 call xc_xm06_d2() 1067 endif 1068 endif 1069c 1070c M08-HX meta GGA Exchange 1071c 1072 if (abs(xfac(48)).gt.eps)then 1073 if (.not. do_2nd) then 1074 call xc_xm11(tol_rho, xfac(48), lxfac(48), nlxfac(48), 1075 & rho, delrho, Amat, Cmat, nq, ipol, 1076 & Ex, qwght,ldew, func,ttau,Mmat,1) 1077 else 1078 call xc_xm11_d2() 1079 endif 1080 endif 1081c 1082c M08-SO meta GGA Exchange 1083c 1084 if (abs(xfac(49)).gt.eps)then 1085 if (.not. do_2nd) then 1086 call xc_xm11(tol_rho, xfac(49), lxfac(49), nlxfac(49), 1087 & rho, delrho, Amat, Cmat, nq, ipol, 1088 & Ex, qwght,ldew, func,ttau,Mmat,2) 1089 else 1090 call xc_xm11_d2() 1091 endif 1092 endif 1093c 1094c M11 meta GGA Exchange 1095c 1096 if (abs(xfac(50)).gt.eps)then 1097 if (.not. do_2nd) then 1098 call xc_xm11(tol_rho, xfac(50), lxfac(50), nlxfac(50), 1099 & rho, delrho, Amat, Cmat, nq, ipol, 1100 & Ex, qwght,ldew, func,ttau,Mmat,3) 1101 else 1102 call xc_xm11_d2() 1103 endif 1104 endif 1105c 1106c M11-L meta GGA Exchange 1107c 1108 if (abs(xfac(51)).gt.eps)then 1109 if (.not. do_2nd) then 1110 call xc_xm11(tol_rho, xfac(51), lxfac(51), nlxfac(51), 1111 & rho, delrho, Amat, Cmat, nq, ipol, 1112 & Ex, qwght,ldew, func,ttau,Mmat,4) 1113 else 1114 call xc_xm11_d2() 1115 endif 1116 endif 1117c 1118c SOGGA GGA Exchange 1119c 1120 if (abs(xfac(52)).gt.eps)then 1121 if (.not. do_2nd) then 1122 call xc_xsogga(tol_rho, xfac(48), lxfac(48), nlxfac(48), 1123 & rho, delrho, Amat, Cmat, nq, ipol, 1124 & Ex, qwght,ldew, func,1) 1125 else 1126 call xc_xsogga_d2() 1127 endif 1128 endif 1129c 1130c SOGGA11 GGA Exchange 1131c 1132 if (abs(xfac(53)).gt.eps)then 1133 if (.not. do_2nd) then 1134 call xc_xsogga(tol_rho, xfac(49), lxfac(49), nlxfac(49), 1135 & rho, delrho, Amat, Cmat, nq, ipol, 1136 & Ex, qwght,ldew, func,2) 1137 else 1138 call xc_xsogga_d2() 1139 endif 1140 endif 1141c 1142c SOGGA11-X GGA Exchange 1143c 1144 if (abs(xfac(54)).gt.eps)then 1145 if (.not. do_2nd) then 1146 call xc_xsogga(tol_rho, xfac(50), lxfac(50), nlxfac(50), 1147 & rho, delrho, Amat, Cmat, nq, ipol, 1148 & Ex, qwght,ldew, func,3) 1149 else 1150 call xc_xsogga_d2() 1151 endif 1152 endif 1153c 1154 if (abs(xfac(55)).gt.eps)then 1155 if (.not. do_2nd) then 1156 call xc_xb86b(tol_rho, xfac(55), lxfac(55), nlxfac(55), 1157 & rho, delrho, Amat, Cmat, nq, ipol, 1158 & Ex, qwght,ldew,func) 1159 else 1160 call xc_xb86b_d2(tol_rho, xfac(7), lxfac(7), nlxfac(7), 1161 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 1162 & Ex, qwght,ldew,func) 1163 endif 1164 endif 1165 if (abs(xfac(56)).gt.eps)then 1166 if (.not. do_2nd) then 1167 call xc_xpw86(tol_rho, xfac(56), lxfac(56), nlxfac(56), 1168 & rho, delrho, Amat, Cmat, nq, ipol, 1169 & Ex, qwght,ldew,func) 1170 else 1171 call xc_xpw86(tol_rho, xfac(56), lxfac(56), nlxfac(56), 1172 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 1173 & Ex, qwght,ldew,func) 1174 endif 1175 endif 1176c 1177c cm08-hx is special in that the GGA part is 1178c defined to be dependent on 1179c the E(LDA) ... 1180c the decision has been made to use the PW91-LDA as the 1181c LDA-correlation. at present, this LDA cannot be 1182c set by the user 1183c 1184 if (abs(cfac(48)).gt.eps)then 1185 if (.not. do_2nd) then 1186 call xc_cm11(tol_rho, cfac(48), lcfac(48), nlcfac(48), 1187 & rho, delrho, nq, ipol, 1188 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,1) 1189 1190 else 1191 call xc_cm11_d2() 1192 endif 1193 endif 1194c 1195c cm08-so is special in that the GGA part is 1196c defined to be dependent on 1197c the E(LDA) ... 1198c the decision has been made to use the PW91-LDA as the 1199c LDA-correlation. at present, this LDA cannot be 1200c set by the user 1201c 1202 if (abs(cfac(49)).gt.eps)then 1203 if (.not. do_2nd) then 1204 call xc_cm11(tol_rho, cfac(49), lcfac(49), nlcfac(49), 1205 & rho, delrho, nq, ipol, 1206 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,2) 1207 1208 else 1209 call xc_cm11_d2() 1210 endif 1211 endif 1212c 1213c cm11 is special in that the GGA part is 1214c defined to be dependent on 1215c the E(LDA) ... 1216c the decision has been made to use the PW91-LDA as the 1217c LDA-correlation. at present, this LDA cannot be 1218c set by the user 1219c 1220 if (abs(cfac(50)).gt.eps)then 1221 if (.not. do_2nd) then 1222 call xc_cm11(tol_rho, cfac(50), lcfac(50), nlcfac(50), 1223 & rho, delrho, nq, ipol, 1224 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,3) 1225 1226 else 1227 call xc_cm11_d2() 1228 endif 1229 endif 1230c 1231c cm11-l is special in that the GGA part is 1232c defined to be dependent on 1233c the E(LDA) ... 1234c the decision has been made to use the PW91-LDA as the 1235c LDA-correlation. at present, this LDA cannot be 1236c set by the user 1237c 1238 if (abs(cfac(51)).gt.eps)then 1239 if (.not. do_2nd) then 1240 call xc_cm11(tol_rho, cfac(51), lcfac(51), nlcfac(51), 1241 & rho, delrho, nq, ipol, 1242 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,4) 1243 1244 else 1245 call xc_cm11_d2() 1246 endif 1247 endif 1248c 1249c csogga is special in that the GGA part is 1250c defined to be dependent on 1251c the E(LDA) ... 1252c the decision has been made to use the PW91-LDA as the 1253c LDA-correlation. at present, this LDA cannot be 1254c set by the user 1255c 1256 if (abs(cfac(52)).gt.eps)then 1257 if (.not. do_2nd) then 1258 call xc_cpbe96(tol_rho, cfac, lcfac, nlcfac, rho, 1259 & delrho, Amat, Cmat, nq, ipol, 1260 & Ec, qwght, ldew, func) 1261 else 1262 call xc_csogga_d2() 1263 endif 1264 endif 1265c 1266c csogga11 is special in that the GGA part is 1267c defined to be dependent on 1268c the E(LDA) ... 1269c the decision has been made to use the PW91-LDA as the 1270c LDA-correlation. at present, this LDA cannot be 1271c set by the user 1272c 1273 if (abs(cfac(53)).gt.eps)then 1274 if (.not. do_2nd) then 1275 call xc_csogga(tol_rho, cfac(49), lcfac(49), nlcfac(49), 1276 & rho, delrho, nq, ipol, 1277 & Ec, qwght, ldew,func,Amat,Cmat,1) 1278 1279 else 1280 call xc_csogga_d2() 1281 endif 1282 endif 1283c 1284c csogga11-x is special in that the GGA part is 1285c defined to be dependent on 1286c the E(LDA) ... 1287c the decision has been made to use the PW91-LDA as the 1288c LDA-correlation. at present, this LDA cannot be 1289c set by the user 1290c 1291 if (abs(cfac(54)).gt.eps)then 1292 if (.not. do_2nd) then 1293 call xc_csogga(tol_rho, cfac(50), lcfac(50), nlcfac(50), 1294 & rho, delrho, nq, ipol, 1295 & Ec, qwght, ldew,func,Amat,Cmat,2) 1296 1297 else 1298 call xc_csogga_d2() 1299 endif 1300 endif 1301c 1302c LC-BNL 2007 Exchange 1303c 1304 if (abs(xfac(38)).gt.eps)then 1305 if (.not. do_2nd .and..not. do_3rd) then 1306 call xc_bnl(tol_rho, xfac(38), lxfac(38), nlxfac(38), 1307 & rho, Amat, nq, ipol, Ex, qwght, ldew, func) 1308 else if (.not. do_3rd) then 1309 call xc_bnl_d2(tol_rho, xfac(38), lxfac(38), nlxfac(38), 1310 & rho, Amat, Amat2, nq, ipol, Ex, qwght, ldew, func) 1311 else 1312 call xc_bnl_d3(tol_rho, xfac(38), lxfac(38), nlxfac(38), 1313 & rho, Amat, Amat2, Amat3, nq, ipol, Ex, qwght, ldew, func) 1314 endif 1315 endif 1316c 1317c CAM-B88 Exchange 1318c 1319 if (abs(xfac(40)).gt.eps)then 1320 if (.not. do_2nd .and..not. do_3rd) then 1321 call xc_camb88(tol_rho, xfac(40), lxfac(40), nlxfac(40), 1322 & rho, delrho, Amat, Cmat, nq, ipol, 1323 & Ex, qwght,ldew,func) 1324 else if (.not. do_3rd) then 1325 call xc_camb88_d2(tol_rho, xfac(40), lxfac(40), nlxfac(40), 1326 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 1327 & Ex, qwght,ldew,func) 1328 else 1329 call xc_camb88_d3(tol_rho, xfac(40), lxfac(40), nlxfac(40), 1330 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 1331 & nq, ipol, Ex, qwght,ldew,func) 1332 endif 1333 endif 1334c 1335c CAM-PBE96 Exchange 1336c 1337 if (abs(xfac(41)).gt.eps)then 1338 if (.not. do_2nd .and..not. do_3rd) then 1339 call xc_camxpbe96('orig', 1340 T tol_rho, xfac(41), lxfac(41), nlxfac(41), 1341 & rho, delrho, Amat, Cmat, nq, ipol, 1342 & Ex, qwght,ldew,func) 1343 else if (.not. do_3rd) then 1344 call xc_camxpbe96_d2('orig', 1345 T tol_rho, xfac(41), lxfac(41), nlxfac(41), 1346 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 1347 & Ex, qwght,ldew,func) 1348 else 1349 call xc_camxpbe96_d3('orig', 1350 T tol_rho, xfac(41), lxfac(41), nlxfac(41), 1351 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 1352 & nq, ipol, Ex, qwght,ldew,func) 1353 endif 1354 endif 1355c 1356c CAM-LSD Exchange 1357c 1358 if (abs(xfac(42)).gt.eps)then 1359 if (.not. do_2nd .and..not. do_3rd) then 1360 call xc_camxlsd(tol_rho, xfac(42), lxfac(42), nlxfac(42), 1361 & rho, Amat, nq, ipol, Ex, qwght, ldew, func) 1362 else if (.not. do_3rd) then 1363 call xc_camxlsd_d2(tol_rho, xfac(42), lxfac(42), nlxfac(42), 1364 & rho, Amat, Amat2, nq, ipol, Ex, qwght, .false., func) 1365 else 1366 call xc_camxlsd_d3(tol_rho, xfac(42), lxfac(42), nlxfac(42), 1367 & rho, Amat, Amat2, Amat3, nq, ipol, Ex, qwght, .false., 1368 & func) 1369 endif 1370 endif 1371c 1372c xwpbe exchange: HSE screened exchange 1373c 1374 if (abs(xfac(43)).gt.eps)then 1375 if (.not. do_2nd .and..not. do_3rd) then 1376 call xc_xwpbe(tol_rho, xfac(43), lxfac(43), nlxfac(43), 1377 & rho, delrho, Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 1378 else if (.not. do_3rd) then 1379 call xc_xwpbe_d2(tol_rho, xfac(43), lxfac(43), nlxfac(43), 1380 & rho, delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 1381 & qwght,ldew,func) 1382 else 1383 call xc_xwpbe_d3(tol_rho, xfac(43), lxfac(43), nlxfac(43), 1384 & rho, delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 1385 & nq, ipol, Ex, qwght, ldew, func) 1386 endif 1387 endif 1388c 1389c cpw6b95 is special in that the GGA part is 1390c defined to be dependent on 1391c the E(LDA) ... 1392c the decision has been made to use the PW91-LDA as the 1393c LDA-correlation. at present, this LDA cannot be 1394c set by the user 1395c Note that bc95, cpw6b95, cpwb6k use the same subroutine xc_bc95() 1396c 1397 if (abs(cfac(32)).gt.eps)then 1398 if (.not. do_2nd) then 1399 call xc_bc95(tol_rho, cfac(32), lcfac(32), nlcfac(32), 1400 & rho, delrho, nq, ipol, 1401 & Ec, qwght, ldew,func, ttau,Amat,Cmat,Mmat,1) 1402 1403 else 1404 call xc_bc95_d2() 1405 endif 1406 endif 1407c 1408c cpwb6k is special in that the GGA part is 1409c defined to be dependent on 1410c the E(LDA) ... 1411c the decision has been made to use the PW91-LDA as the 1412c LDA-correlation. at present, this LDA cannot be 1413c set by the user 1414c 1415 if (abs(cfac(33)).gt.eps)then 1416 if (.not. do_2nd) then 1417 1418 call xc_bc95(tol_rho, cfac(33), lcfac(33), nlcfac(33), 1419 & rho, delrho, nq, ipol, 1420 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,2) 1421 1422 else 1423 call xc_bc95_d2() 1424 endif 1425 endif 1426c 1427c cm05 is special in that the GGA part is 1428c defined to be dependent on 1429c the E(LDA) ... 1430c the decision has been made to use the PW91-LDA as the 1431c LDA-correlation. at present, this LDA cannot be 1432c set by the user 1433c 1434 if (abs(cfac(34)).gt.eps)then 1435 if (.not. do_2nd) then 1436 call xc_cm05(tol_rho, cfac(34), lcfac(34), nlcfac(34), 1437 & rho, delrho, nq, ipol, 1438 & Ec, qwght, ldew, func,ttau,Amat,Cmat,Mmat,1) 1439 1440 else 1441 call xc_cm05_d2() 1442 endif 1443 endif 1444c 1445c cm05-2x is special in that the GGA part is 1446c defined to be dependent on 1447c the E(LDA) ... 1448c the decision has been made to use the PW91-LDA as the 1449c LDA-correlation. at present, this LDA cannot be 1450c set by the user 1451c 1452 if (abs(cfac(35)).gt.eps)then 1453 if (.not. do_2nd) then 1454 call xc_cm05(tol_rho, cfac(35), lcfac(35), nlcfac(35), 1455 & rho, delrho, nq, ipol, 1456 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,2) 1457 1458 else 1459 call xc_cm05_d2() 1460 endif 1461 endif 1462c 1463c dlDF Correlation 1464c 1465c cdldf is special in that the GGA part is 1466c defined to be dependent on 1467c the E(LDA) ... 1468c the decision has been made to use the PW91-LDA as the 1469c LDA-correlation. at present, this LDA cannot be 1470c set by the user 1471c 1472 if (abs(cfac(42)).gt.eps)then 1473 if (.not. do_2nd) then 1474 call xc_cdldf(tol_rho, cfac(42), lcfac(42), nlcfac(42), 1475 & rho, delrho, nq, ipol, 1476 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat) 1477 1478 else 1479 call xc_cdldf_d2() 1480 endif 1481 endif 1482c 1483c cvs98 is special in that the GGA part is 1484c defined to be dependent on 1485c the E(LDA) ... 1486c the decision has been made to use the PW91-LDA as the 1487c LDA-correlation. at present, this LDA cannot be 1488c set by the user 1489c 1490 if (abs(cfac(29)).gt.eps)then 1491 if (.not. do_2nd) then 1492 call xc_cvs98(tol_rho, cfac(29), lcfac(29), nlcfac(29), 1493 & rho, delrho, nq, ipol, 1494 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,1) 1495 1496 else 1497 call xc_cvs98_d2() 1498 endif 1499 endif 1500c 1501c cm06-L is special in that the GGA part is 1502c defined to be dependent on 1503c the E(LDA) ... 1504c the decision has been made to use the PW91-LDA as the 1505c LDA-correlation. at present, this LDA cannot be 1506c set by the user 1507c 1508 if (abs(cfac(30)).gt.eps)then 1509 if (.not. do_2nd) then 1510 call xc_cm06(tol_rho, cfac(30), lcfac(30), nlcfac(30), 1511 & rho, delrho, nq, ipol, 1512 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,1) 1513 1514 else 1515 call xc_cm06_d2() 1516 endif 1517 endif 1518c 1519c cm06-hf is special in that the GGA part is 1520c defined to be dependent on 1521c the E(LDA) ... 1522c the decision has been made to use the PW91-LDA as the 1523c LDA-correlation. at present, this LDA cannot be 1524c set by the user 1525c 1526 if (abs(cfac(38)).gt.eps)then 1527 if (.not. do_2nd) then 1528 call xc_cm06(tol_rho, cfac(38), lcfac(38), nlcfac(38), 1529 & rho, delrho, nq, ipol, 1530 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,2) 1531 1532 else 1533 call xc_cm06_d2() 1534 endif 1535 endif 1536c 1537c cm06 is special in that the GGA part is 1538c defined to be dependent on 1539c the E(LDA) ... 1540c the decision has been made to use the PW91-LDA as the 1541c LDA-correlation. at present, this LDA cannot be 1542c set by the user 1543c 1544 if (abs(cfac(39)).gt.eps)then 1545 if (.not. do_2nd) then 1546 call xc_cm06(tol_rho, cfac(39), lcfac(39), nlcfac(39), 1547 & rho, delrho, nq, ipol, 1548 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,3) 1549 1550 else 1551 call xc_cm06_d2() 1552 endif 1553 endif 1554c 1555c cm06-2x is special in that the GGA part is 1556c defined to be dependent on 1557c the E(LDA) ... 1558c the decision has been made to use the PW91-LDA as the 1559c LDA-correlation. at present, this LDA cannot be 1560c set by the user 1561c 1562 if (abs(cfac(40)).gt.eps)then 1563 if (.not. do_2nd) then 1564 call xc_cm06(tol_rho, cfac(40), lcfac(40), nlcfac(40), 1565 & rho, delrho, nq, ipol, 1566 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,4) 1567 1568 else 1569 call xc_cm06_d2() 1570 endif 1571 endif 1572c 1573c N12x 1574c 1575 if (abs(xfac(45)).gt.eps)then 1576 if (.not. do_2nd) then 1577 call xc_xn12(tol_rho, xfac(45), lxfac(45), nlxfac(45), 1578 & rho, delrho, Amat, Cmat, nq, ipol, 1579 & Ex, qwght, ldew,func,1) 1580 1581 else 1582 call xc_xn12_d2() 1583 endif 1584 endif 1585c 1586c revM06 correlation 1587c 1588 if (abs(cfac(68)).gt.eps)then 1589 if (.not. do_2nd) then 1590 call xc_cm06(tol_rho, cfac(68), lcfac(68), nlcfac(68), 1591 & rho, delrho, nq, ipol, 1592 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,6) 1593 1594 else 1595 call xc_cm06_d2() 1596 endif 1597 endif 1598c 1599c revM06-L correlation 1600c 1601 if (abs(cfac(69)).gt.eps)then 1602 if (.not. do_2nd) then 1603 call xc_cm06(tol_rho, cfac(69), lcfac(69), nlcfac(69), 1604 & rho, delrho, nq, ipol, 1605 & Ec, qwght, ldew,func,ttau,Amat,Cmat,Mmat,5) 1606 1607 else 1608 call xc_cm06_d2() 1609 endif 1610 endif 1611c 1612c SCAN meta GGA 1613c 1614 if (abs(xfac(66)).gt.eps) then 1615 if(.not. do_2nd) then 1616 call xc_xscan('orig', tol_rho, xfac(66), rho, delrho, Amat, 1617 & Cmat, nq, ipol, Ex, qwght, ldew, func, ttau, 1618 & Mmat) 1619 else 1620 call xc_xscan_d2() 1621 endif 1622 endif 1623c 1624 if (abs(cfac(66)).gt.eps) then 1625 if(.not. do_2nd) then 1626 call xc_cscan('orig', tol_rho, cfac(66), rho, delrho, Amat, 1627 & Cmat, nq, ipol, Ec, qwght, ldew, func, ttau, 1628 & Mmat) 1629 else 1630 call xc_cscan_d2() 1631 endif 1632 endif 1633c 1634c regularized SCAN 1635c 1636 if (abs(xfac(71)).gt.eps) then 1637 if(.not. do_2nd) then 1638 call xc_xscan('regu', tol_rho, xfac(71), rho, delrho, Amat, 1639 & Cmat, nq, ipol, Ex, qwght, ldew, func, ttau, 1640 & Mmat) 1641 else 1642 call xc_xscan_d2() 1643 endif 1644 endif 1645c 1646 if (abs(cfac(71)).gt.eps) then 1647 if(.not. do_2nd) then 1648 call xc_cscan('regu', tol_rho, cfac(71), rho, delrho, Amat, 1649 & Cmat, nq, ipol, Ec, qwght, ldew, func, ttau, 1650 & Mmat) 1651 else 1652 call xc_cscan_d2() 1653 endif 1654 endif 1655c 1656 if (abs(xfac(72)).gt.eps)then 1657 if (.not. do_2nd .and..not. do_3rd) then 1658 call xc_xncap(tol_rho, xfac(72), rho, 1659 & delrho, Amat, Cmat, nq, ipol, 1660 & Ex, qwght, ldew, func) 1661 else if (.not. do_3rd) then 1662 call xc_xncap_d2(tol_rho, xfac(72), rho, 1663 & delrho, Amat, Amat2, Cmat, Cmat2, nq, ipol, 1664 & Ex, qwght, ldew, func) 1665 else 1666 call xc_xncap_d3(tol_rho, xfac(72), rho, 1667 & delrho, Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 1668 & nq, ipol, Ex, qwght, ldew, func) 1669 endif 1670 endif 1671c 1672c SCAN-L Laplacian meta GGA 1673c 1674 if (abs(xfac(67)).gt.eps) then 1675 if(.not. do_2nd) then 1676 call xc_xscanl(tol_rho, xfac(67), rho, delrho, laprho, 1677 & Amat, Cmat, Lmat, nq, ipol, Ex, qwght, ldew, 1678 & func) 1679 else 1680 call xc_xscanl_d2() 1681 endif 1682 endif 1683c 1684 if (abs(cfac(67)).gt.eps) then 1685 if(.not. do_2nd) then 1686 call xc_cscanl(tol_rho, cfac(67), rho, delrho, laprho, 1687 & Amat, Cmat, Lmat, nq, ipol, Ec, qwght, ldew, 1688 & func) 1689 else 1690 call xc_cscanl_d2() 1691 endif 1692 endif 1693c 1694c Calculate the steric energy 1695c 1696 if (lsteric) then 1697 StericEnergy = 0.d0 1698 call steric_energy(tol_rho,xfac(1),rho,delrho,nq, 1699 & qwght,ipol,StericEnergy) 1700 endif 1701c 1702c Combine with quadrature weights 1703c 1704c Daniel (1-11-13): Added XC-third derivative stuff. This currently 1705c doesn't include any functionality for meta-GGAs. 1706c Daniel (1-17-13): Because we need derivatives of the quadrature 1707c weights for TDDFT gradients, we don't want to multiply the A 1708c and C matrices with the quadrature weights here (which would be 1709c double counting). This prevents us from multiplying the 1710c quadrature weights into the matrices and then dividing them out 1711c later, which is prone to numerical problems (and slow). 1712c Daniel (2-4-13): ldew3 is for the dfxc*(X+Y)*(X+Y) term in the TDDFT 1713c gradients. ldew2 is for the dVxc*P term in the TDDFT gradients. 1714 if (.not.(ldew2.or.ldew3)) then 1715 1716 if (.not. do_2nd .and. .not. do_3rd) then 1717 call setACmat(delrho, Amat, Cmat, qwght, ipol, nq, grad, 1718 & (.not. do_2nd), kske, Mmat, kslap, Lmat) 1719 else if (.not. do_3rd) then 1720 call setACmat_d2(delrho, Amat, Amat2, Cmat, Cmat2, qwght, 1721 & ipol, nq, grad, (.not. do_2nd), kske, Mmat, Mmat2, 1722 & .false.) 1723 else 1724 call setACmat_d3(delrho, Amat, Amat2, Amat3, Cmat, Cmat2, 1725 & Cmat3, qwght, ipol, nq, grad, (.not. do_3rd), .false., 1726 & .false.) 1727 endif 1728 endif 1729c 1730 return 1731 end 1732 1733cc AJL/Begin/FDE 1734c 1735 Subroutine xc_eval_fnl_fde(rho, delrho, Amat, Amat2, Cmat, Cmat2, 1736 & nq, Ex, Ec, qwght, grad, ldew, func, 1737 & do_2nd, ttau, kske, Mmat, Mmat2, 1738 & StericEnergy, do_3rd, Amat3, Cmat3, ldew2) 1739c 1740 implicit none 1741c 1742#include "cdft.fh" 1743c xc-second derivative header file 1744c#include "dft2drv.fh" 1745c xc-third derivative header file 1746c#include "dft3drv.fh" 1747c#include "stdio.fh" 1748#include "steric.fh" 1749c 1750 integer nq 1751 double precision rho(*) !< [Input] The electron density 1752 !< \f$\rho^\alpha\f$ and 1753 !< \f$\rho^\beta\f$. 1754 double precision delrho(*) !< [Input] The electron density gradient 1755 !< \f$\partial\rho^\alpha/\partial r\f$ 1756 !< and 1757 !< \f$\partial\rho^\beta/\partial r\f$ 1758c 1759 double precision Amat(*), Cmat(*) 1760 double precision Amat2(*),Cmat2(*) 1761 double precision Amat3(*),Cmat3(*) 1762 double precision ttau(*) !< [Input] The kinetic energy density 1763 double precision Mmat(*), Mmat2(*) 1764c 1765 double precision Ex !< [Output] The exchange energy 1766 double precision Ec !< [Output] The correlation energy 1767 double precision StericEnergy 1768 double precision qwght(nq), func(nq) 1769 logical grad, kske 1770c 1771 logical do_2nd, do_3rd 1772 logical ldew, ldew2 1773c Local 1774 double precision cfac_temp(numfunc), xfac_temp(numfunc) 1775 logical lcfac_temp(numfunc), nlcfac_temp(numfunc) 1776 logical lxfac_temp(numfunc), nlxfac_temp(numfunc) 1777 integer i 1778c 1779 do i=1,numfunc 1780 cfac_temp(i) = cfac(i) 1781 xfac_temp(i) = xfac(i) 1782 lcfac_temp(i) = lcfac(i) 1783 nlcfac_temp(i) = nlcfac(i) 1784 lxfac_temp(i) = lxfac(i) 1785 nlxfac_temp(i) = nlxfac(i) 1786 1787 cfac(i) = cfac_fde(i) 1788 xfac(i) = xfac_fde(i) 1789 lcfac(i) = lcfac_fde(i) 1790 nlcfac(i) = nlcfac_fde(i) 1791 lxfac(i) = lxfac_fde(i) 1792 nlxfac(i) = nlxfac_fde(i) 1793 enddo 1794 1795 call xc_eval_fnl(rho, delrho, Amat, Amat2, Cmat, Cmat2, 1796 & nq, Ex, Ec, qwght, grad, ldew, func, 1797 & do_2nd, ttau, kske, Mmat, Mmat2, 1798 & 0d0, .false., 0d0, 1799 & StericEnergy, do_3rd, Amat3, Cmat3, ldew2) 1800 1801 do i=1,numfunc 1802 cfac(i) = cfac_temp(i) 1803 xfac(i) = xfac_temp(i) 1804 lcfac(i) = lcfac_temp(i) 1805 nlcfac(i) = nlcfac_temp(i) 1806 lxfac(i) = lxfac_temp(i) 1807 nlxfac(i) = nlxfac_temp(i) 1808 enddo 1809 1810 return 1811 end 1812c 1813cc AJL/End 1814 1815 subroutine calc_rgamma(ipol,nq,delrho,rgamma) 1816 implicit none 1817cinclude "nwxc_param.fh" 1818#define G_TT 1 1819#define G_AA 1 1820#define G_AB 2 1821#define G_BB 3 1822 integer ipol !< [Input] The number of spin channels 1823 integer nq !< [Input] The number of grid points 1824 double precision delrho(nq,3,ipol) !< [Input] The density gradient 1825 double precision rgamma(nq,3) !< [Output] The density gradient norm 1826c 1827 integer iq 1828c 1829 if (ipol.eq.1) then 1830 do iq = 1, nq 1831 rgamma(iq,G_TT) = delrho(iq,1,1)*delrho(iq,1,1) 1832 & + delrho(iq,2,1)*delrho(iq,2,1) 1833 & + delrho(iq,3,1)*delrho(iq,3,1) 1834 enddo 1835 else 1836 do iq = 1, nq 1837 rgamma(iq,G_AA) = delrho(iq,1,1)*delrho(iq,1,1) 1838 & + delrho(iq,2,1)*delrho(iq,2,1) 1839 & + delrho(iq,3,1)*delrho(iq,3,1) 1840 rgamma(iq,G_BB) = delrho(iq,1,2)*delrho(iq,1,2) 1841 & + delrho(iq,2,2)*delrho(iq,2,2) 1842 & + delrho(iq,3,2)*delrho(iq,3,2) 1843 rgamma(iq,G_AB) = delrho(iq,1,1)*delrho(iq,1,2) 1844 & + delrho(iq,2,1)*delrho(iq,2,2) 1845 & + delrho(iq,3,1)*delrho(iq,3,2) 1846 enddo 1847 endif 1848c 1849 end 1850C> 1851C> @} 1852c $Id$ 1853