1 subroutine nwchem_libxc_compute(nq,ex,ec,qwght,rho,delrho,ttau, 2 $ laprho,Amat,Amat2,Amat3,Cmat, 3 $ Cmat2,Cmat3, 4 $ Mmat,Lmat,func, 5 $ ksgrad,kske,kslap,ldew, 6 $ do_2nd, do_3rd) 7 use, intrinsic :: iso_c_binding 8 9#ifdef USE_LIBXC 10 use xc_f03_lib_m 11#endif 12 implicit none 13 14#include "libxc.fh" 15#include "mafdecls.fh" 16#include "errquit.fh" 17#include "cdft.fh" 18#include "dft2drv.fh" 19#include "dft3drv.fh" 20 21 logical :: ksgrad, kske, kslap, ldew, do_2nd, do_3rd 22 23 integer :: q,nq 24 integer :: ifunc 25 integer :: size1, size2 26 integer :: lexc,kexc 27 integer :: lrho,krho,lvrho,kvrho 28 integer :: lsigma,ksigma,lvsigma,kvsigma 29 integer :: ltau,ktau,lvtau,kvtau 30 integer :: llapl,klapl,lvlapl,kvlapl 31 32 integer :: kv2rho2,lv2rho2 33 integer :: kv2sig2,lv2sig2 34 integer :: kv2rhosig,lv2rhosig 35 36 integer :: kv3rho3,lv3rho3 37 integer :: kv3sig3,lv3sig3 38 integer :: kv3rho2sig,lv3rho2sig 39 integer :: kv3rhosig2,lv3rhosig2 40 41 integer(c_int) :: polarized 42 43 double precision :: ex,ec,func(nq),qwght(nq) 44 double precision :: rho(nq,(ipol*(ipol+1))/2) 45 double precision :: delrho(nq,3,ipol) 46 double precision :: laprho(nq,ipol) 47 double precision :: ttau(nq,ipol) 48 49 double precision :: Amat(nq,ipol) 50 double precision :: Cmat(nq,3) 51 double precision :: Mmat(nq,ipol) 52 double precision :: Lmat(nq,ipol) 53 54 double precision :: Amat2(nq,NCOL_AMAT2) 55 double precision :: Cmat2(nq,NCOL_CMAT2) 56 57 double precision :: Amat3(nq,NCOL_AMAT3) 58 double precision :: Cmat3(nq,NCOL_CMAT3) 59 60 double precision :: fac 61 62 double precision, external :: ddot 63 logical gga,mgga,dolap 64 logical,external :: nwchem_libxc_family 65 66 integer(c_size_t) :: nqs 67#ifdef USE_LIBXC 68 type(xc_f03_func_t) :: xcfunc 69 70 gga = .false. 71 mgga = .false. 72 dolap = .false. 73 74 nqs = nq 75 size1 = nq*ipol 76 size2 = nq*(ipol*(ipol+1))/2 77 78 if (.not.(ma_push_get(mt_dbl,size1,'rhoval',lrho,krho))) 79 $ call errquit(" could not allocate",0,ma_err) 80 if (.not.(ma_push_get(mt_dbl,size1,'vrho',lvrho,kvrho))) 81 $ call errquit(" could not allocate",0,ma_err) 82 if (.not.(ma_push_get(mt_dbl,nq,'exc',lexc,kexc))) 83 $ call errquit(" could not allocate",0,ma_err) 84 85 if (ksgrad) then 86 if (.not.(ma_push_get(mt_dbl,size2,'sigma',lsigma,ksigma))) 87 $ call errquit(" could not allocate",0,ma_err) 88 if (.not.(ma_push_get(mt_dbl,size2,'vsigma',lvsigma,kvsigma))) 89 $ call errquit(" could not allocate",0,ma_err) 90 endif 91 92 if (kske.or.kslap) then 93 if (.not.(ma_push_get(mt_dbl,size2,'tau',ltau,ktau))) 94 $ call errquit(" could not allocate",0,ma_err) 95 if (.not.(ma_push_get(mt_dbl,size2,'vtau',lvtau,kvtau))) 96 $ call errquit(" could not allocate",0,ma_err) 97 if (.not.(ma_push_get(mt_dbl,size2,'lapl',llapl,klapl))) 98 $ call errquit(" could not allocate",0,ma_err) 99 if (.not.(ma_push_get(mt_dbl,size2,'vlapl',lvlapl,kvlapl))) 100 $ call errquit(" could not allocate",0,ma_err) 101 endif 102 103 if (do_2nd .or. do_3rd) then 104 if (.not.(ma_push_get(mt_dbl,3*nq,'v2rho2',lv2rho2,kv2rho2))) 105 $ call errquit(" could not allocate v2rho2",0,ma_err) 106 if (ksgrad) then 107 if (.not.(ma_push_get(mt_dbl,6*nq,'v2rhosig',lv2rhosig, 108 $ kv2rhosig))) 109 $ call errquit(" could not allocate v2rhosig",0,ma_err) 110 if (.not.(ma_push_get(mt_dbl,6*nq,'v2sig2',lv2sig2,kv2sig2))) 111 $ call errquit(" could not allocate v2sig2",0,ma_err) 112 endif 113 endif 114 115 if (do_3rd) then 116 if (.not.(ma_push_get(mt_dbl,4*nq,'v3rho3',lv3rho3,kv3rho3))) 117 $ call errquit(" could not allocate v3rho3",0,ma_err) 118 if (ksgrad) then 119 if (.not.(ma_push_get(mt_dbl,9*nq,'v3rho2sig',lv3rho2sig, 120 $ kv3rho2sig))) 121 $ call errquit(" could not allocate v3rho2sig",0,ma_err) 122 if (.not.(ma_push_get(mt_dbl,12*nq,'v3rhosig2',lv3rhosig2, 123 $ kv3rhosig2))) 124 $ call errquit(" could not allocate v3rhosig2",0,ma_err) 125 if (.not.(ma_push_get(mt_dbl,10*nq,'v3sig3',lv3sig3,kv3sig3))) 126 $ call errquit(" could not allocate v3sig3",0,ma_err) 127 endif 128 endif 129 130 if (kske) call dcopy(nq,ttau,1,dbl_mb(ktau),ipol) 131 if (kslap) call dcopy(nq,laprho,1,dbl_mb(klapl),ipol) 132 133 if (ipol.eq.1) then 134 polarized = xc_unpolarized 135 136 call dcopy(nq,rho,1,dbl_mb(krho),1) 137 if (ksgrad) then 138 call nwchem_libxc_util1(nq,dbl_mb(ksigma),delrho) 139 endif 140 else 141 polarized = xc_polarized 142 call dcopy(nq,rho(1,2),1,dbl_mb(krho),2) 143 call dcopy(nq,rho(1,3),1,dbl_mb(krho+1),2) 144 if (ksgrad) then 145 call nwchem_libxc_util2(nq,dbl_mb(ksigma),delrho) 146 endif 147 if (kske) call dcopy(nq,ttau(1,2),1,dbl_mb(ktau+1),ipol) 148 if (kslap) call dcopy(nq,laprho(1,2),1,dbl_mb(klapl+1),ipol) 149 endif 150 151 do ifunc=1,libxc_nfuncs 152 153 fac = libxc_facts(ifunc) 154 155 call xc_f03_func_init(xcfunc,libxc_funcs(ifunc),polarized) 156 call xc_f03_func_set_dens_threshold(xcfunc, tol_rho) 157 call xc_f03_func_set_sigma_threshold(xcfunc, tol_rho**2) 158 call xc_f03_func_set_zeta_threshold(xcfunc, 1d-10) 159 160 select case(libxc_family(ifunc)) 161 case (XC_FAMILY_LDA, XC_FAMILY_HYB_LDA) 162 163 if ((.not.do_2nd) .and. (.not.do_3rd)) then 164 call xc_f03_lda_exc_vxc(xcfunc,nqs, 165 $ dbl_mb(krho),dbl_mb(kexc),dbl_mb(kvrho)) 166 elseif (.not.do_3rd) then 167 call xc_f03_lda_exc_vxc_fxc(xcfunc,nqs,dbl_mb(krho), 168 $ dbl_mb(kexc),dbl_mb(kvrho),dbl_mb(kv2rho2)) 169 else 170 call xc_f03_lda_exc_vxc_fxc_kxc(xcfunc,nqs,dbl_mb(krho), 171 $ dbl_mb(kexc),dbl_mb(kvrho),dbl_mb(kv2rho2), 172 $ dbl_mb(kv3rho3)) 173 endif 174 175 176 case (XC_FAMILY_GGA, XC_FAMILY_HYB_GGA) 177 178 gga = .true. 179 180 if (iand(libxc_flags(ifunc),xc_flags_have_exc).eq. 181 $ xc_flags_have_exc) then 182 183 if ((.not.do_2nd).and.(.not.do_3rd)) then 184 call xc_f03_gga_exc_vxc(xcfunc,nqs, 185 $ dbl_mb(krho),dbl_mb(ksigma), 186 $ dbl_mb(kexc), 187 $ dbl_mb(kvrho),dbl_mb(kvsigma)) 188 elseif (.not.do_3rd) then 189 call xc_f03_gga_exc_vxc_fxc(xcfunc,nqs, 190 $ dbl_mb(krho),dbl_mb(ksigma),dbl_mb(kexc), 191 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kv2rho2), 192 $ dbl_mb(kv2rhosig),dbl_mb(kv2sig2)) 193 else 194 call xc_f03_gga_exc_vxc_fxc_kxc(xcfunc,nqs, 195 $ dbl_mb(krho),dbl_mb(ksigma),dbl_mb(kexc), 196 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kv2rho2), 197 $ dbl_mb(kv2rhosig),dbl_mb(kv2sig2), 198 $ dbl_mb(kv3rho3),dbl_mb(kv3rho2sig), 199 $ dbl_mb(kv3rhosig2),dbl_mb(kv3sig3)) 200 endif 201 else 202 if ((.not.do_2nd).and.(.not.do_3rd)) then 203 call xc_f03_gga_vxc(xcfunc,nqs, 204 $ dbl_mb(krho),dbl_mb(ksigma), 205 $ dbl_mb(kvrho),dbl_mb(kvsigma)) 206 elseif (.not.do_3rd) then 207 call xc_f03_gga_vxc_fxc(xcfunc,nqs, 208 $ dbl_mb(krho),dbl_mb(ksigma), 209 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kv2rho2), 210 $ dbl_mb(kv2rhosig),dbl_mb(kv2sig2)) 211 else 212 call xc_f03_gga_vxc_fxc_kxc(xcfunc,nqs, 213 $ dbl_mb(krho),dbl_mb(ksigma), 214 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kv2rho2), 215 $ dbl_mb(kv2rhosig),dbl_mb(kv2sig2), 216 $ dbl_mb(kv3rho3),dbl_mb(kv3rho2sig), 217 $ dbl_mb(kv3rhosig2),dbl_mb(kv3sig3)) 218 endif 219 call dfill(nq,0d0,dbl_mb(kexc),1) 220 endif 221 222 case (XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA) 223 gga = .true. 224 mgga = .true. 225 dolap = .true. 226 227 228 if (iand(libxc_flags(ifunc),xc_flags_have_exc).eq. 229 $ xc_flags_have_exc) then 230 231 call xc_f03_mgga_exc_vxc(xcfunc,nqs, 232 $ dbl_mb(krho),dbl_mb(ksigma),dbl_mb(klapl),dbl_mb(ktau), 233 $ dbl_mb(kexc), 234 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kvlapl),dbl_mb(kvtau)) 235 else 236 call xc_f03_mgga_vxc(xcfunc,nqs, 237 $ dbl_mb(krho),dbl_mb(ksigma),dbl_mb(klapl),dbl_mb(ktau), 238 $ dbl_mb(kvrho),dbl_mb(kvsigma),dbl_mb(kvlapl),dbl_mb(kvtau)) 239 call dfill(nq,0d0,dbl_mb(kexc),1) 240 endif 241 242 end select 243 244 call xc_f03_func_end(xcfunc) 245 call nwchem_libxc_util3(ldew,nq,fac, 246 e dbl_mb(kexc),func,qwght,rho) 247 248 if (libxc_kind(ifunc).eq.xc_exchange .or. 249 $ libxc_kind(ifunc).eq.xc_exchange_correlation) then 250 ex = ex + fac*ddot(nq,dbl_mb(kexc),1,rho,1) 251 elseif(libxc_kind(ifunc).eq.xc_correlation) then 252 ec = ec + fac*ddot(nq,dbl_mb(kexc),1,rho,1) 253 endif 254 255 call daxpy(nq,fac,dbl_mb(kvrho),ipol,Amat,1) 256 if (mgga) call daxpy(nq,0.5d0*fac,dbl_mb(kvtau),ipol,Mmat,1) 257 if (dolap) call daxpy(nq,fac,dbl_mb(kvlapl),ipol,Lmat,1) 258 259 if ((polarized.eq.xc_unpolarized) .and. (gga)) then 260 call daxpy(nq,fac,dbl_mb(kvsigma),1,Cmat(1,D1_GAA),1) 261 call daxpy(nq,2d0*fac,dbl_mb(kvsigma),1,Cmat(1,D1_GAB),1) 262 elseif (polarized.eq.xc_polarized) then 263 call daxpy(nq,fac,dbl_mb(kvrho+1),2,Amat(1,2),1) 264 if (mgga) call daxpy(nq,0.5d0*fac,dbl_mb(kvtau+1),2, 265 $ Mmat(1,2),1) 266 if (dolap) call daxpy(nq,fac,dbl_mb(kvlapl+1),2,Lmat(1,2),1) 267 if (gga) then 268 call daxpy(nq,fac,dbl_mb(kvsigma),3,Cmat(1,D1_GAA),1) 269 call daxpy(nq,fac,dbl_mb(kvsigma+1),3,Cmat(1,D1_GAB),1) 270 call daxpy(nq,fac,dbl_mb(kvsigma+2),3,Cmat(1,D1_GBB),1) 271 endif 272 endif 273 274 if (do_2nd .or. do_3rd) then 275 if (polarized.eq.xc_unpolarized) then 276 call daxpy(nq,fac,dbl_mb(kv2rho2),1,Amat2(1,D2_RA_RA),1) 277 call daxpy(nq,fac,dbl_mb(kv2rho2),1,Amat2(1,D2_RA_RB),1) 278 if (gga) then 279 call daxpy(nq,fac,dbl_mb(kv2rhosig),1, 280 $ Cmat2(1,D2_RA_GAA),1) 281 call daxpy(nq,2d0*fac,dbl_mb(kv2rhosig),1, 282 $ Cmat2(1,D2_RA_GAB),1) 283 call daxpy(nq,fac,dbl_mb(kv2rhosig),1, 284 $ Cmat2(1,D2_RA_GBB),1) 285 call daxpy(nq,fac,dbl_mb(kv2sig2),1, 286 $ Cmat2(1,D2_GAA_GAA),1) 287 call daxpy(nq,2d0*fac,dbl_mb(kv2sig2),1, 288 $ Cmat2(1,D2_GAA_GAB),1) 289 call daxpy(nq,fac,dbl_mb(kv2sig2),1, 290 $ Cmat2(1,D2_GAA_GBB),1) 291 call daxpy(nq,4d0*fac,dbl_mb(kv2sig2),1, 292 $ Cmat2(1,D2_GAB_GAB),1) 293 endif 294 else 295 call daxpy(nq,fac,dbl_mb(kv2rho2),3,Amat2(1,D2_RA_RA),1) 296 call daxpy(nq,fac,dbl_mb(kv2rho2+1),3,Amat2(1,D2_RA_RB),1) 297 call daxpy(nq,fac,dbl_mb(kv2rho2+2),3,Amat2(1,D2_RB_RB),1) 298 if (gga) then 299 call daxpy(nq,fac,dbl_mb(kv2rhosig),6, 300 $ Cmat2(1,D2_RA_GAA),1) 301 call daxpy(nq,fac,dbl_mb(kv2rhosig+1),6, 302 $ Cmat2(1,D2_RA_GAB),1) 303 call daxpy(nq,fac,dbl_mb(kv2rhosig+2),6, 304 $ Cmat2(1,D2_RA_GBB),1) 305 call daxpy(nq,fac,dbl_mb(kv2rhosig+3),6, 306 $ Cmat2(1,D2_RB_GAA),1) 307 call daxpy(nq,fac,dbl_mb(kv2rhosig+4),6, 308 $ Cmat2(1,D2_RB_GAB),1) 309 call daxpy(nq,fac,dbl_mb(kv2rhosig+5),6, 310 $ Cmat2(1,D2_RB_GBB),1) 311 312 call daxpy(nq,fac,dbl_mb(kv2sig2),6, 313 $ Cmat2(1,D2_GAA_GAA),1) 314 call daxpy(nq,fac,dbl_mb(kv2sig2+1),6, 315 $ Cmat2(1,D2_GAA_GAB),1) 316 call daxpy(nq,fac,dbl_mb(kv2sig2+2),6, 317 $ Cmat2(1,D2_GAA_GBB),1) 318 call daxpy(nq,fac,dbl_mb(kv2sig2+3),6, 319 $ Cmat2(1,D2_GAB_GAB),1) 320 call daxpy(nq,fac,dbl_mb(kv2sig2+4),6, 321 $ Cmat2(1,D2_GAB_GBB),1) 322 call daxpy(nq,fac,dbl_mb(kv2sig2+5),6, 323 $ Cmat2(1,D2_GBB_GBB),1) 324 endif 325 endif 326 endif 327 328 if (do_3rd) then 329 if (polarized.eq.xc_unpolarized) then 330 call daxpy(nq,fac,dbl_mb(kv3rho3),1,Amat3(1,D3_RA_RA_RA),1) 331 call daxpy(nq,fac,dbl_mb(kv3rho3),1,Amat3(1,D3_RA_RA_RB),1) 332 call daxpy(nq,fac,dbl_mb(kv3rho3),1,Amat3(1,D3_RA_RB_RB),1) 333 if (gga) then 334 call daxpy(nq,fac,dbl_mb(kv3rho2sig),1, 335 $ Cmat3(1,D3_RA_RA_GAA),1) 336 call daxpy(nq,2d0*fac,dbl_mb(kv3rho2sig),1, 337 $ Cmat3(1,D3_RA_RA_GAB),1) 338 call daxpy(nq,fac,dbl_mb(kv3rho2sig),1, 339 $ Cmat3(1,D3_RA_RA_GBB),1) 340 341 call daxpy(nq,fac,dbl_mb(kv3rho2sig),1, 342 $ Cmat3(1,D3_RA_RB_GAA),1) 343 call daxpy(nq,2d0*fac,dbl_mb(kv3rho2sig),1, 344 $ Cmat3(1,D3_RA_RB_GAB),1) 345 call daxpy(nq,fac,dbl_mb(kv3rho2sig),1, 346 $ Cmat3(1,D3_RA_RB_GBB),1) 347 348 call daxpy(nq,fac,dbl_mb(kv3rhosig2),1, 349 $ Cmat3(1,D3_RA_GAA_GAA),1) 350 call daxpy(nq,2d0*fac,dbl_mb(kv3rhosig2),1, 351 $ Cmat3(1,D3_RA_GAA_GAB),1) 352 call daxpy(nq,fac,dbl_mb(kv3rhosig2),1, 353 $ Cmat3(1,D3_RA_GAA_GBB),1) 354 call daxpy(nq,4d0*fac,dbl_mb(kv3rhosig2),1, 355 $ Cmat3(1,D3_RA_GAB_GAB),1) 356 call daxpy(nq,2d0*fac,dbl_mb(kv3rhosig2),1, 357 $ Cmat3(1,D3_RA_GAB_GBB),1) 358 call daxpy(nq,fac,dbl_mb(kv3rhosig2),1, 359 $ Cmat3(1,D3_RA_GBB_GBB),1) 360 361 call daxpy(nq,fac,dbl_mb(kv3sig3),1, 362 $ Cmat3(1,D3_GAA_GAA_GAA),1) 363 call daxpy(nq,2d0*fac,dbl_mb(kv3sig3),1, 364 $ Cmat3(1,D3_GAA_GAA_GAB),1) 365 call daxpy(nq,fac,dbl_mb(kv3sig3),1, 366 $ Cmat3(1,D3_GAA_GAA_GBB),1) 367 call daxpy(nq,4d0*fac,dbl_mb(kv3sig3),1, 368 $ Cmat3(1,D3_GAA_GAB_GAB),1) 369 call daxpy(nq,2d0*fac,dbl_mb(kv3sig3),1, 370 $ Cmat3(1,D3_GAA_GAB_GBB),1) 371 call daxpy(nq,fac,dbl_mb(kv3sig3),1, 372 $ Cmat3(1,D3_GAA_GBB_GBB),1) 373 call daxpy(nq,8d0*fac,dbl_mb(kv3sig3),1, 374 $ Cmat3(1,D3_GAB_GAB_GAB),1) 375 endif 376 elseif (polarized.eq.xc_polarized) then 377 call daxpy(nq,fac,dbl_mb(kv3rho3),4,Amat3(1,D3_RA_RA_RA),1) 378 call daxpy(nq,fac,dbl_mb(kv3rho3+1),4,Amat3(1,D3_RA_RA_RB),1) 379 call daxpy(nq,fac,dbl_mb(kv3rho3+2),4,Amat3(1,D3_RA_RB_RB),1) 380 call daxpy(nq,fac,dbl_mb(kv3rho3+3),4,Amat3(1,D3_RB_RB_RB),1) 381 if (gga) then 382 call daxpy(nq,fac,dbl_mb(kv3rho2sig),9, 383 $ Cmat3(1,D3_RA_RA_GAA),1) 384 call daxpy(nq,fac,dbl_mb(kv3rho2sig+1),9, 385 $ Cmat3(1,D3_RA_RA_GAB),1) 386 call daxpy(nq,fac,dbl_mb(kv3rho2sig+2),9, 387 $ Cmat3(1,D3_RA_RA_GBB),1) 388 call daxpy(nq,fac,dbl_mb(kv3rho2sig+3),9, 389 $ Cmat3(1,D3_RA_RB_GAA),1) 390 call daxpy(nq,fac,dbl_mb(kv3rho2sig+4),9, 391 $ Cmat3(1,D3_RA_RB_GAB),1) 392 call daxpy(nq,fac,dbl_mb(kv3rho2sig+5),9, 393 $ Cmat3(1,D3_RA_RB_GBB),1) 394 call daxpy(nq,fac,dbl_mb(kv3rho2sig+6),9, 395 $ Cmat3(1,D3_RB_RB_GAA),1) 396 call daxpy(nq,fac,dbl_mb(kv3rho2sig+7),9, 397 $ Cmat3(1,D3_RB_RB_GAB),1) 398 call daxpy(nq,fac,dbl_mb(kv3rho2sig+8),9, 399 $ Cmat3(1,D3_RB_RB_GBB),1) 400 401 call daxpy(nq,fac,dbl_mb(kv3rhosig2),12, 402 $ Cmat3(1,D3_RA_GAA_GAA),1) 403 call daxpy(nq,fac,dbl_mb(kv3rhosig2+1),12, 404 $ Cmat3(1,D3_RA_GAA_GAB),1) 405 call daxpy(nq,fac,dbl_mb(kv3rhosig2+2),12, 406 $ Cmat3(1,D3_RA_GAA_GBB),1) 407 call daxpy(nq,fac,dbl_mb(kv3rhosig2+3),12, 408 $ Cmat3(1,D3_RA_GAB_GAB),1) 409 call daxpy(nq,fac,dbl_mb(kv3rhosig2+4),12, 410 $ Cmat3(1,D3_RA_GAB_GBB),1) 411 call daxpy(nq,fac,dbl_mb(kv3rhosig2+5),12, 412 $ Cmat3(1,D3_RA_GBB_GBB),1) 413 call daxpy(nq,fac,dbl_mb(kv3rhosig2+6),12, 414 $ Cmat3(1,D3_RB_GAA_GAA),1) 415 call daxpy(nq,fac,dbl_mb(kv3rhosig2+7),12, 416 $ Cmat3(1,D3_RB_GAA_GAB),1) 417 call daxpy(nq,fac,dbl_mb(kv3rhosig2+8),12, 418 $ Cmat3(1,D3_RB_GAA_GBB),1) 419 call daxpy(nq,fac,dbl_mb(kv3rhosig2+9),12, 420 $ Cmat3(1,D3_RB_GAB_GAB),1) 421 call daxpy(nq,fac,dbl_mb(kv3rhosig2+10),12, 422 $ Cmat3(1,D3_RB_GAB_GBB),1) 423 call daxpy(nq,fac,dbl_mb(kv3rhosig2+11),12, 424 $ Cmat3(1,D3_RB_GBB_GBB),1) 425 426 call daxpy(nq,fac,dbl_mb(kv3sig3),10, 427 $ Cmat3(1,D3_GAA_GAA_GAA),1) 428 call daxpy(nq,fac,dbl_mb(kv3sig3+1),10, 429 $ Cmat3(1,D3_GAA_GAA_GAB),1) 430 call daxpy(nq,fac,dbl_mb(kv3sig3+2),10, 431 $ Cmat3(1,D3_GAA_GAA_GBB),1) 432 call daxpy(nq,fac,dbl_mb(kv3sig3+3),10, 433 $ Cmat3(1,D3_GAA_GAB_GAB),1) 434 call daxpy(nq,fac,dbl_mb(kv3sig3+4),10, 435 $ Cmat3(1,D3_GAA_GAB_GBB),1) 436 call daxpy(nq,fac,dbl_mb(kv3sig3+5),10, 437 $ Cmat3(1,D3_GAA_GBB_GBB),1) 438 call daxpy(nq,fac,dbl_mb(kv3sig3+6),10, 439 $ Cmat3(1,D3_GAB_GAB_GAB),1) 440 call daxpy(nq,fac,dbl_mb(kv3sig3+7),10, 441 $ Cmat3(1,D3_GAB_GAB_GBB),1) 442 call daxpy(nq,fac,dbl_mb(kv3sig3+8),10, 443 $ Cmat3(1,D3_GAB_GBB_GBB),1) 444 call daxpy(nq,fac,dbl_mb(kv3sig3+9),10, 445 $ Cmat3(1,D3_GBB_GBB_GBB),1) 446 endif 447 endif 448 endif 449 450 enddo 451 452 if (.not.ma_chop_stack(lrho)) 453 $ call errquit(" could not chop stack",0,ma_err) 454#endif 455 end subroutine 456 subroutine nwchem_libxc_util1(nq,sigma,delrho) 457 implicit none 458 integer nq 459 double precision sigma(*) 460 double precision delrho(nq,3,*) 461c 462 integer q 463c 464 do q=1,nq 465 sigma(q) = delrho(q,1,1)**2 + 466 $ delrho(q,2,1)**2 + 467 $ delrho(q,3,1)**2 468 enddo 469 end subroutine 470 subroutine nwchem_libxc_util2(nq,sigma,delrho) 471 implicit none 472 integer nq 473 double precision sigma(3,*) 474 double precision delrho(nq,3,*) 475c 476 integer q 477c 478 do q=1,nq 479 sigma(1,q) = delrho(q,1,1)**2 + 480 $ delrho(q,2,1)**2 + 481 $ delrho(q,3,1)**2 482 sigma(2,q) = delrho(q,1,1)*delrho(q,1,2) + 483 $ delrho(q,2,1)*delrho(q,2,2) + 484 $ delrho(q,3,1)*delrho(q,3,2) 485 sigma(3,q) = delrho(q,1,2)**2 + 486 $ delrho(q,2,2)**2 + 487 $ delrho(q,3,2)**2 488 enddo 489 end subroutine 490 subroutine nwchem_libxc_util3(ldew,nq,fac, 491 e exc,func,qwght,rho) 492 implicit none 493 logical ldew 494 integer nq 495 double precision fac 496 double precision exc(nq),qwght(nq), 497 r rho(nq,*),func(nq) 498c 499 integer q 500c 501 if(ldew) then 502 do q=1,nq 503 func(q) = func(q) + fac*exc(q)*rho(q,1) 504 enddo 505 endif 506 do q=1,nq 507 exc(q) = exc(q)*qwght(q) 508 enddo 509 end subroutine 510