1c $Id$ 2c 3C> \ingroup nwdft 4C> @{ 5C> 6C> \file xc_getv.F 7C> Calculate exchange-correlation energy 8C> 9C> \brief Calculate the exchange-correlation energy and Fock matrix 10C> contributions 11C> 12C> This driver routine solves for the XC energy and potential (Vxc) via 13C> numerical quadrature methods. The results are obtained either by 14C> direct numerical integration or by means of a LSQ fit of the Vxc to 15C> a set of Gaussian functions. This fitted function can be used to 16C> evaluate Vxc via a summation of a series of 3-center overlap 17C> integrals (3OIs). The algorithms are formulated in terms of matrix 18C> products. See subsequent subroutines for further explanation. 19C> 20 Subroutine xc_getv(rtdb, Exc, ecoul,nExc, iVxc_opt, g_xcinv, 21 & g_dens, g_vxc, IOLGC, g_wght, g_xyz,g_nq, 22 & wght_GA, rho_n, rdens_atom, 23 & cetobfr, natoms) 24c 25 implicit none 26#include "errquit.fh" 27c 28 integer nExc !< [Input] The number of energy terms 29 !< - nExc=1: Exc(1) = exchange + correlation 30 !< - nExc=2: Exc(1) = exchange, 31 !< Exc(2) = correlation 32 integer iVxc_opt !< [Input] If 1 then do density fitting for 33 !< exchange 34 integer g_xcinv !< [Work] GA for the inversion of the fitting 35 !< matrix 36 integer g_dens(2) !< [Input] The density matrices, if ipol=1 37 !< g_dens(1)=\f$D^\alpha+D^\beta\f$, else 38 !< g_dens(1)=\f$D^\alpha\f$ and 39 !< g_dens(2)=\f$D^\beta\f$. 40 integer g_vxc(4) !< [Output] DFT Fock matrix contributions, if 41 !< ipol=1 g_vxc(1)=\f$F^\alpha+F^\beta\f$, else 42 !< g_vxc(1)=\f$F^\alpha\f$ and 43 !< g_vxc(2)=\f$F^\beta\f$. 44 integer g_wght !< [Work] The grid point weights if wght_GA 45 integer g_xyz !< [Work] The grid point coordinates if wght_GA 46 integer g_nq !< [Unused] 47 integer natoms !< [Input] The number of atoms 48 logical IOLGC !< [Input] .TRUE. do not use disk for exchange 49 !< fitting, .FALSE. store data on disk 50 logical wght_GA !< [Input] .TRUE. store grid points in GA, 51 !< .FALSE. store grid points on file 52 integer rtdb !< [Input] The RTDB handle 53c 54#include "mafdecls.fh" 55#include "rtdb.fh" 56#include "bas.fh" 57#include "global.fh" 58#include "tcgmsg.fh" 59#include "cdft.fh" 60#include "oep.fh" 61#include "dftpara.fh" 62#include "util.fh" 63#include "sym.fh" 64#include "stdio.fh" 65#include "case.fh" 66#include "dftps.fh" 67c 68 integer cetobfr(2,natoms) !< [Unused] 69 double precision rho_n !< [Output] The number of electrons 70 !< obtained by integrating the density 71 double precision rdens_atom(ipol*natoms*natoms) !< [Unused] 72 double precision jfac(4),kfac(4) 73 integer g_jk(4), g_d(4) 74 logical havehfxc 75c 76 integer ga_create_atom_blocked 77cc AJL/Begin/FDE 78c logical xc_gotxc 79c external ga_create_atom_blocked,xc_gotxc 80 external ga_create_atom_blocked 81cc AJL/End 82c 83c--> XC Energy 84c 85 double precision Exc(2) !< [Output] The energy terms 86 !< - nExc=1: Exc(1) = exchange + 87 !< correlation 88 !< - nExc=2: Exc(1) = exchange, 89 !< Exc(2) = correlation 90 double precision ecoul !< [Output] The Coulomb energy 91c 92c This driver routine solves for the XC energy and potential (Vxc) via 93c numerical quadrature methods. The results are obtained either by direct 94c numerical integration or by means of a LSQ fit of the Vxc to a set of 95c Gaussian functions. This fitted function can be used to evaluate Vxc 96c via a summation of a series of 3-center overlap integrals (3OIs). The 97c algorithms are formulated in terms of matrix products. See subsequent 98c subroutines for further explanation. 99c 100c XC Energy and Potential Index Key, Vxc(pq,i) 101c 102c Value of | Definition of index "i" 103c ipol nExc | 1 2 3 4 104c -------------------------------------------------- 105c 1 1 | Vxc 106c 2 1 | Vxc^up Vxc^dw 107c 1 2 | Vxc 108c 2 2 | Vxc^up Vxc^dw 109c 110c nTcols = ipol 111c 112cc AJL/Begin/FDE 113c integer me,nTrows,nTcols 114 integer me 115c integer lTmat,iTmat,g_truevxc(2) 116c AJL: These parameters are unused? 117c double precision zero,one,onem 118 logical oprint_intermediate_xc, oprint_time 119cc AJL/Unused 120c , oprint_oep 121c parameter(zero=0.d0,one=1.d0,onem=-1.d0) 122cc AJL/End 123 double precision tol2e 124 integer g_tmp(2) 125c 126c timings 127c 128 double precision time1_2e,time2_2e 129cc AJL/Begin/FDE 130c double precision time1_xc,time2_xc 131cc AJL/End 132c 133c****************************************************************************** 134c 135c Compute the matrix elements for the XC potential and energy. 136c 137 oprint_intermediate_xc = util_print('intermediate XC matrix', 138 $ print_debug) 139 oprint_time = util_print('dft timings', print_high) 140cc AJL/Unused 141c oprint_oep = util_print('oep', print_high) 142 Exc(1)=0.d0 143 Exc(2)=0.d0 144cc AJL/Begin/FDE 145c iTmat=0 146cc AJL/End 147c 148 me=ga_nodeid() 149 havehfxc=abs(xfac(1)).gt.1d-8 150c 151 if (oprint_intermediate_xc)then 152c write(luout,*)' rtdb, Exc, nExc, iVxc_opt, g_xcinv: ', 153c & rtdb, Exc, nExc, iVxc_opt, g_xcinv 154c write(luout,*)'g_dens(1),g_vxc(1),IOLGC,g_wght,g_xyz,wght_GA:', 155c & g_dens(1),g_vxc(1),IOLGC,g_wght,g_xyz,wght_GA 156 write(luout,*)' Fock XC matrix entering xc_getv: ' 157 call ga_print(g_vxc(1)) 158 if(ipol.eq.2)call ga_print(g_vxc(2)) 159c call ga_print(g_dens(1)) 160c if(ipol.eq.2)call ga_print(g_dens(2)) 161 endif 162c 163 if(oprint_time) 164 & time1_2e=util_cpusec() ! start 2e build time 165 if (havehfxc .or. (.not. CDFIT))then 166c 167c Compute the exact exchange potential (as in Hartree-Fock calculations). 168c 169 tol2e=10.d0**(-itol2e) 170 call ga_sync 171 if (odftps) call pstat_on(ps_f2e) 172 if (oprint_time)call dft_tstamp(' Before call to fock_2e. ') 173 if (ipol.eq.1) then 174 if (.not. CDFIT) then 175 if (.not.cam_exch) then ! for regular calculations 176c 177c set up prefactors 178 kfac(1) = -0.5d0*xfac(1) 179 jfac(1) = 0.0d0 180 jfac(2) = 1.0d0 181 kfac(2) = 0.0d0 182c 183c get some work space 184 g_vxc(2) = ga_create_atom_blocked(geom,ao_bas_han,'jk') 185c 186c calculate the exchange and coulomb parts 187 call ga_zero(g_vxc(2)) 188 g_dens(2)=g_dens(1) 189 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 190 & tol2e, oskel, g_dens, g_vxc, .false.) 191 Exc(1) = Exc(1)+0.5d0*ga_ddot(g_dens(1),g_vxc(1)) 192 ecoul = 0.5d0*ga_ddot(g_dens(1),g_vxc(2)) 193 call ga_dadd(1d0,g_vxc(1),1d0,g_vxc(2),g_vxc(1)) 194 if (.not. ga_destroy(g_vxc(2))) call errquit 195 $ ('xc_getv: ga corrupt?',0, GA_ERR) 196 else ! CAM calculations 197c 198c get some work space 199 g_tmp(1)=ga_create_atom_blocked(geom,ao_bas_han,'work') 200 call ga_zero(g_tmp(1)) 201c 202 g_tmp(2)=ga_create_atom_blocked(geom,ao_bas_han,'work') 203 call ga_zero(g_tmp(2)) 204c 205c set up prefactors for exchange 206 kfac(1) = -0.5d0*xfac(1) 207 jfac(1) = 0.0d0 208 kfac(2) = 0.0d0 209 jfac(2) = 0.0d0 210 g_dens(2)=g_dens(1) 211 call case_setflags(.true.) 212 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 213 & tol2e, oskel, g_dens, g_tmp, .false.) 214 Exc(1) = Exc(1)+0.5d0*ga_ddot(g_dens(1),g_tmp(1)) 215 call ga_dadd(1d0,g_vxc(1),1d0,g_tmp(1),g_vxc(1)) 216 call case_setflags(.false.) 217c 218c calculate the full Coulomb 219 call ga_zero(g_tmp(1)) 220 call ga_zero(g_tmp(2)) 221c 222c set up prefactors for coulomb 223 kfac(1) = 0.0d0 224 jfac(1) = 1.0d0 225 kfac(2) = 0.0d0 226 jfac(2) = 0.0d0 227 g_dens(2)=g_dens(1) 228 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 229 & tol2e, oskel, g_dens, g_tmp, .false.) 230 ecoul = 0.5d0*ga_ddot(g_dens(1),g_tmp(1)) 231 call ga_dadd(1d0,g_vxc(1),1d0,g_tmp(1),g_vxc(1)) 232c 233c destroy work space 234 if (.not. ga_destroy(g_tmp(1))) call errquit 235 $ ('xc_getv: ga corrupt?',0, GA_ERR) 236 if (.not. ga_destroy(g_tmp(2))) call errquit 237 $ ('xc_getv: ga corrupt?',0, GA_ERR) 238 end if 239 else ! with CDFIT 240c 241c set up prefactors 242 kfac(1) = -0.5d0*xfac(1) 243 jfac(1) = 0.0d0 244c 245c calculate the non-CAM exchange 246 if (cam_exch) call case_setflags(.true.) 247 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 248 & tol2e, oskel, g_dens(1), g_vxc(1), .false.) 249 if (cam_exch) call case_setflags(.false.) ! turn off attenuation 250 Exc(1) = Exc(1)+0.5d0*ga_ddot(g_dens(1),g_vxc(1)) 251 endif 252 else ! spin-polarized calculations 253 if (CDFIT) then 254 jfac(1)=0.d0 255 jfac(2)=0.d0 256 kfac(1)=-1.0d0*xfac(1) 257 kfac(2)=-1.0d0*xfac(1) 258 if (cam_exch) call case_setflags(.true.) 259 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 260 & tol2e, oskel, g_dens, g_vxc, .false.) 261 if (cam_exch) call case_setflags(.false.) ! turn off attenuation 262 Exc(1) = Exc(1)+0.5d0*(ga_ddot(g_dens(1),g_vxc(1)) + 263 & ga_ddot(g_dens(2),g_vxc(2))) 264 else 265 if (.not.cam_exch) then ! for regular calculations 266 jfac(1) = 1.0d0 267 jfac(2) = 0.0d0 268 jfac(3) = 1.0d0 269 jfac(4) = 0.0d0 270 kfac(1) = 0.0d0 271 kfac(2) = 1.0d0 272 kfac(3) = 0.0d0 273 kfac(4) = 1.0d0 274 g_jk(1) = g_vxc(1) ! This assignment is assumed 275 g_jk(2) = g_vxc(2) 276 g_jk(3) = ga_create_atom_blocked(geom, ao_bas_han, 'jk') 277 g_jk(4) = ga_create_atom_blocked(geom, ao_bas_han, 'jk') 278 call ga_zero(g_jk(3)) 279 call ga_zero(g_jk(4)) 280 g_d(1) = g_dens(1) 281 g_d(2) = g_dens(1) 282 g_d(3) = g_dens(2) 283 g_d(4) = g_dens(2) 284 call fock_2e(geom, AO_bas_han, 4, jfac, kfac, 285 & tol2e, oskel, g_d(1), g_jk(1), .false.) 286 ecoul = 0.5d0*( ! Alpha coulomb energy 287 $ ga_ddot(g_dens(1),g_jk(1)) + 288 $ ga_ddot(g_dens(1),g_jk(3))) 289 ecoul = ecoul + 0.5d0*( ! Beta coulomb energy 290 $ ga_ddot(g_dens(2),g_jk(1)) + 291 $ ga_ddot(g_dens(2),g_jk(3))) 292 exc(1) = exc(1) - xfac(1)*0.5d0*( ! All exchange energy 293 $ ga_ddot(g_dens(1),g_jk(2)) + 294 $ ga_ddot(g_dens(2),g_jk(4))) 295 call ga_dadd(1.0d0, g_jk(1), 1.0d0, g_jk(3), g_jk(1)) 296 call ga_copy(g_jk(1), g_jk(3)) 297 call ga_dadd(1.0d0, g_jk(1), -xfac(1), g_jk(2), 298 $ g_jk(1)) 299 call ga_dadd(1.0d0, g_jk(3), -xfac(1), g_jk(4), 300 $ g_jk(2)) 301 if (.not. ga_destroy(g_jk(3))) call errquit 302 $ ('xc_getv: ga corrupt?',0, GA_ERR) 303 if (.not. ga_destroy(g_jk(4))) call errquit 304 $ ('xc_getv: ga corrupt?',1, GA_ERR) 305 else 306c 307c Allocate some scratch space 308 g_tmp(1)=ga_create_atom_blocked(geom, ao_bas_han,'tmp1') 309 g_tmp(2)=ga_create_atom_blocked(geom, ao_bas_han,'tmp2') 310c 311c Calculate Coulomb 312 jfac(1) = 1.0d0 313 jfac(2) = 1.0d0 314 kfac(1) = 0.0d0 315 kfac(2) = 0.0d0 316 call ga_zero(g_tmp(1)) 317 call ga_zero(g_tmp(2)) 318 call case_setflags(.false.) 319 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 320 & tol2e, oskel, g_dens, g_tmp, .false.) 321c 322c Accumulate contribution 323 call ga_dadd(1.0d0, g_vxc(1), 1.0d0, g_tmp(1), g_vxc(1)) 324 call ga_dadd(1.0d0, g_vxc(2), 1.0d0, g_tmp(2), g_vxc(2)) 325 call ga_dadd(1.0d0, g_vxc(1), 1.0d0, g_vxc(2), g_vxc(1)) 326 call ga_copy(g_vxc(1), g_vxc(2)) 327 ecoul = 0.5d0*( ! Alpha coulomb energy 328 $ ga_ddot(g_dens(1),g_tmp(1)) + 329 $ ga_ddot(g_dens(1),g_tmp(2))) 330 ecoul = ecoul + 0.5d0*( ! Beta coulomb energy 331 $ ga_ddot(g_dens(2),g_tmp(1)) + 332 $ ga_ddot(g_dens(2),g_tmp(2))) 333c 334c Calculate Exchange 335 jfac(1) = 0.0d0 336 jfac(2) = 0.0d0 337 kfac(1) =-1.0d0*xfac(1) 338 kfac(2) =-1.0d0*xfac(1) 339 call ga_zero(g_tmp(1)) 340 call ga_zero(g_tmp(2)) 341 call case_setflags(.true.) ! turn on attenuation 342 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 343 & tol2e, oskel, g_dens, g_tmp, .false.) 344 call case_setflags(.false.) ! turn off attenuation 345c 346c Accumulate contribution 347 call ga_dadd(1.0d0, g_vxc(1), 1.0d0, g_tmp(1), g_vxc(1)) 348 call ga_dadd(1.0d0, g_vxc(2), 1.0d0, g_tmp(2), g_vxc(2)) 349 exc(1) = exc(1) + 0.5d0*( ! Exchange energy 350 $ ga_ddot(g_dens(1),g_tmp(1)) + 351 $ ga_ddot(g_dens(2),g_tmp(2))) 352c 353c Deallocate scratch 354 if (.not. ga_destroy(g_tmp(1))) call errquit 355 $ ('xc_getv: ga corrupt?',0, GA_ERR) 356 if (.not. ga_destroy(g_tmp(2))) call errquit 357 $ ('xc_getv: ga corrupt?',1, GA_ERR) 358c 359 end if 360 endif 361 endif 362 if (odftps) call pstat_off(ps_f2e) 363 if (oprint_time)call dft_tstamp(' After call to fock_2e. ') 364 call ga_sync 365 endif 366 if(oprint_time) 367 & time2_2e=util_cpusec() ! end 2e build time 368c 369c print fock_2e build time 370c 371 if(oprint_time) then 372 if (me.eq.0) then 373 write(luout,"(4x,'Fock_2e Build Time:',F13.1,'s')") 374 & time2_2e-time1_2e 375 endif 376 end if 377c 378c Get the DFT exchange-correlation contribution 379c 380cc AJl/Begin/FDE 381c I have moved this to a new subroutine so the XC evaluation can be 382c called multiple times, as is needed for the FDE evaluation of the 383c non-additive XC energy 384c 385c if(util_print('dft timings', print_high)) 386c & time1_xc=util_cpusec() ! start xc build time 387c if (xc_gotxc()) then 388c if(xcfit) then 389c nTrows = nbf_xc 390c nTcols = ipol 391c if (.not.ma_push_get(MT_Dbl,nTrows*nTcols,'Tmat',lTmat, 392c & iTmat))call errquit('xc_getv: cannot allocate Tmat',0, 393c & MA_ERR) 394c call dfill(nTrows*nTcols,0.D0,dbl_mb(iTmat),1) 395c endif 396c 397c if(havehfxc.or.(.not.cdfit)) then 398c if(.not.ga_duplicate(g_vxc(1),g_truevxc(1),'g vxc 1')) 399c . call errquit('xcgetv: gaduplicate failed',1, GA_ERR) 400c call ga_zero(g_truevxc(1)) 401c if(ipol.eq.2) then 402c if(.not.ga_duplicate(g_vxc(2),g_truevxc(2),'gv21')) 403c . call errquit('xcgetv: gaduplicate failed',1, GA_ERR) 404c call ga_zero(g_truevxc(2)) 405c endif 406c else 407c g_truevxc(1)=g_vxc(1) 408c g_truevxc(2)=g_vxc(2) 409c endif 410cc 411c call grid_quadv0(rtdb, g_dens, g_truevxc, 412c & nexc,rho_n, Exc, dbl_mb(itmat)) 413cc 414c if(havehfxc.or.(.not.cdfit)) then 415c call ga_dadd(1d0,g_vxc(1),1d0,g_truevxc(1),g_vxc(1)) 416c if (.not. ga_destroy(g_truevxc(1))) call errquit( 417c & ' xc_getv: ga_destroy failed ',0, GA_ERR) 418c if(ipol.eq.2) then 419c call ga_dadd(1d0,g_vxc(2),1d0,g_truevxc(2),g_vxc(2)) 420c if (.not. ga_destroy(g_truevxc(2))) call errquit( 421c & ' xc_getv: ga_destroy failed ',0, GA_ERR) 422c endif 423c endif 424c if(util_print('dft timings', print_high)) 425c & time2_xc=util_cpusec() ! end xc build time 426cc 427cc print fock_xc build time 428c if(util_print('dft timings', print_high)) then 429c if (me.eq.0) then 430c write(*,"(4x,'Fock_xc Build Time:',F13.1,'s')") 431c & time2_xc-time1_xc 432c endif 433c end if 434cc 435cc In case we are performing an xc fit calculation 436c if(xcfit) then 437cc 438cc symmetrize the "T" vector 439cc 440c if (oskel)then 441c call sym_vec_symmetrize( 442c . geom,xc_bas_han,Dbl_MB(iTmat)) 443c if (ipol.gt.1)then 444c call sym_vec_symmetrize(geom, xc_bas_han, 445c & Dbl_MB(iTmat+nbf_xc)) 446c endif 447c endif 448c call xc_fitv(rtdb,Dbl_MB(iTmat), nTrows, nTcols, 449c & g_vxc, g_xcinv, IOLGC) 450c if (.not.ma_pop_stack(lTmat)) 451c & call errquit('xc_getv: cannot pop stack',0, MA_ERR) 452cc 453c endif 454c endif 455cc 456 if (odftps) call pstat_on(ps_getvxc) 457 call xc_getvxc(rtdb, Exc, nExc, iVxc_opt, g_xcinv, 458 & g_dens, g_vxc, IOLGC, rho_n, 0, 0) 459 if (odftps) call pstat_off(ps_getvxc) 460c last two inputs -------------------------------->.fde_option., g_dens_fde 461c 462cc AJL/End 463c 464 if (oprint_intermediate_xc)then 465 write(luout,*)' Fock XC matrix leaving xc_getv: ' 466 call util_flush(6) 467 call ga_print(g_vxc(1)) 468 if(ipol.eq.2)call ga_print(g_vxc(2)) 469 call util_flush(6) 470 endif 471c 472 return 473 end 474C> 475C> @} 476