1C> 2C> \brief Calculate Fukui indeces 3C> 4C> Calculate Fukui indeces for open shell wavefunctions [1,2]. 5C> 6C> [1] E. Chamorro, P. Pérez, 7C> <i>"Condensed-to-atoms electronic Fukui functions within the 8C> framework of spin-polarized density-functional theory"</i>, 9C> J. Chem. Phys. <b>123</b>, 114107 (2005), DOI: 10C> <a href="https://doi.org/10.1063/1.2033689"> 11C> 10.1063/1.2033689</a>. 12C> 13C> [2] R.G. Parr, W. Yang, 14C> <i>"Density functional approach to the frontier-electron theory 15C> of chemical reactivity"</i>, 16C> J. Am. Chem. Soc. <b>106</b>, 4049-4050 (1984), DOI: 17C> <a href="https://doi.org/10.1021/ja00326a036"> 18C> 10.1021/ja00326a036</a>. 19C> 20 Subroutine fukui(g_movecs, k_eval, tol2e, rtdb, 21 & nExc, iVxc_opt, g_xcinv, IOLGC, g_wght, 22 & g_xyz, g_nq, wght_GA, rho_n, irdens_atom, 23 & icetobfr, natoms) 24c 25c $Id$ 26c 27 implicit none 28#include "errquit.fh" 29 integer g_dens_HOMO(2), g_dens_LUMO(2), g_orb, g_dens_ss, ik, 30 & g_movecs(2), isp, k_eval(2), g_s, 31 & me, l_temp_vec, i_temp_vec, 32 & g_dens(2), noc_aux_1, noc_aux_2, noc_test 33 integer irdens_atom, icetobfr 34c 35 integer nExc 36 integer iVxc_opt 37 integer g_xcinv, g_vxc(4), g_wght, g_xyz,g_nq 38 integer natoms 39 logical IOLGC, wght_GA 40 integer rtdb 41c 42 double precision eig_lumo(2), eig_homo(2), jfac(4), kfac(4), 43 & mu_n_mas, mu_n_men, mu_n_cer, tol2e, 44 & mu_s_mas, mu_s_men, mu_s_cer, 45 & int_HaHa, int_HbHb, int_LaLa, int_LbLb, 46 & int_HaLb, int_LaHb, Exc(2), ecoul, rho_n, 47 & int_vxc_H(2), int_vxc_L(2), Exc_zero, 48 & Exc_pert, diff_Exc, ion_pot, ele_afi, 49 & high_mult, low_mult, e_orbital, e_coul, 50 & e_xc 51c 52c #include "stdio.fh" 53c #include "mafdecls.fh" 54c #include "global.fh" 55c #include "util.fh" 56c 57#include "mafdecls.fh" 58#include "rtdb.fh" 59#include "bas.fh" 60#include "global.fh" 61#include "tcgmsg.fh" 62#include "cdft.fh" 63#include "util.fh" 64#include "sym.fh" 65c 66#include "stdio.fh" 67c 68 integer ga_create_atom_blocked 69 external ga_create_atom_blocked 70 logical oprint_fukui 71 me = ga_nodeid() 72 oprint_fukui = util_print('Fukui information', print_high) 73 g_orb = ga_create_atom_blocked(geom, AO_bas_han, 'ga_orb') 74 if (.not.MA_Push_Get(MT_Dbl, nbf_ao, 'temp vec', 75 & l_temp_vec, i_temp_vec)) 76 & call errquit('fukui: cannot allocate temp vec',0, MA_ERR) 77 do isp = 1, 2 78 g_dens_HOMO(isp) = ga_create_atom_blocked(geom, 79 & AO_bas_han, 'ga_dens_orb_homo') 80 g_dens_LUMO(isp) = ga_create_atom_blocked(geom, 81 & AO_bas_han, 'ga_dens_orb_lumo') 82 call ga_zero(g_dens_HOMO(isp)) 83 call ga_zero(g_dens_LUMO(isp)) 84 end do 85c 86 g_dens_ss = ga_create_atom_blocked(geom, AO_bas_han, 87 & 'ga_dens_orb_ss') 88c 89 do isp = 1, ipol 90 if (noc(isp).eq.0) then 91 noc_aux_1 = 1 92 noc_aux_2 = 1 93 else 94 noc_aux_1 = noc(isp) 95 noc_aux_2 = noc(isp) + 1 96 end if 97C do ik = noc(isp), (noc(isp)+1) 98 do ik = noc_aux_1, noc_aux_2 99 call ga_zero(g_orb) 100 call ga_get(g_movecs(isp), 1, nbf_ao, ik, ik, 101 & Dbl_MB(i_temp_vec), nbf_ao) 102 call ga_put(g_orb, 1, nbf_ao, ik, ik, 103 & Dbl_MB(i_temp_vec), nbf_ao) 104 if (ik.eq.noc(isp)) then 105 eig_homo(isp) = dbl_mb(k_eval(isp) + 106 & noc(isp) - 1) 107 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, nbf_ao, 108 & 1.0d00, g_orb, g_orb, 0.d00, 109 & g_dens_HOMO(isp)) 110 else 111 eig_lumo(isp) = dbl_mb(k_eval(isp) + noc(isp)) 112 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, nbf_ao, 113 & 1.0d00, g_orb, g_orb, 0.d00, 114 & g_dens_LUMO(isp)) 115 end if 116 end do 117 end do 118c 119 if (.not.ma_pop_stack(l_temp_vec)) 120 & call errquit('fukui: cannot pop stack',0, MA_ERR) 121c 122 if (ipol.eq.1) then 123 call ga_copy(g_dens_HOMO(1), g_dens_HOMO(2)) 124 call ga_copy(g_dens_LUMO(1), g_dens_LUMO(2)) 125 eig_homo(2) = eig_homo(1) 126 eig_lumo(2) = eig_lumo(1) 127 end if 128c 129 if (me.eq.0.and.oprint_fukui) 130 & call dft_header(' Reactivity Parameters ') 131c 132 mu_n_mas = 0.5d00*(eig_lumo(1) + eig_lumo(2)) 133 mu_n_men = 0.5d00*(eig_homo(1) + eig_homo(2)) 134 mu_n_cer = 0.5*(mu_n_mas + mu_n_men) 135 mu_s_mas = 0.5d00*(eig_lumo(1) - eig_homo(2)) 136 mu_s_men = 0.5d00*(eig_homo(1) - eig_lumo(2)) 137 mu_s_cer = 0.5*(mu_s_mas + mu_s_men) 138c 139 g_s = ga_create_atom_blocked(geom, AO_bas_han, 'AO ovl') 140 call ga_zero(g_s) 141 call int_1e_ga(ao_bas_han,ao_bas_han,g_s,'overlap',.false.) 142c 143 call ga_zero(g_dens_ss) 144 call ga_add(0.5d00, g_dens_LUMO(2), 0.5d00, 145 & g_dens_LUMO(1), g_dens_ss) 146 if (me.eq.0.and.oprint_fukui) 147 & call dft_header(' Condensed Fukui function [fnn(+)]') 148 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 149c 150 call ga_zero(g_dens_ss) 151 call ga_add(0.5d00, g_dens_HOMO(2), 0.5d00, 152 & g_dens_HOMO(1), g_dens_ss) 153 if (me.eq.0.and.oprint_fukui) 154 & call dft_header(' Condensed Fukui function [fnn(-)]') 155 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 156c 157 call ga_zero(g_dens_ss) 158 call ga_add(-0.5d00, g_dens_LUMO(2), 0.5d00, 159 & g_dens_LUMO(1), g_dens_ss) 160 if (me.eq.0.and.oprint_fukui) 161 & call dft_header(' Condensed Fukui function [fsn(+)]') 162 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 163c 164 call ga_zero(g_dens_ss) 165 call ga_add(0.5d00, g_dens_HOMO(2), 0.5d00, 166 & g_dens_LUMO(1), g_dens_ss) 167 if (me.eq.0.and.oprint_fukui) 168 & call dft_header(' Condensed Fukui function [fss(+)]') 169 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 170c 171 call ga_zero(g_dens_ss) 172 call ga_add(-0.5d00, g_dens_HOMO(2), 0.5d00, 173 & g_dens_LUMO(1), g_dens_ss) 174 if (me.eq.0.and.oprint_fukui) 175 & call dft_header(' Condensed Fukui function [fns(+)]') 176 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 177c 178 call ga_zero(g_dens_ss) 179 call ga_add(0.5d00, g_dens_HOMO(1), 0.5d00, 180 & g_dens_LUMO(2), g_dens_ss) 181 if (me.eq.0.and.oprint_fukui) 182 & call dft_header(' Condensed Fukui function [fss(-)]') 183 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 184c 185 call ga_zero(g_dens_ss) 186 call ga_add(0.5d00, g_dens_HOMO(1), -0.5d00, 187 & g_dens_LUMO(2), g_dens_ss) 188 if (me.eq.0.and.oprint_fukui) 189 & call dft_header(' Condensed Fukui function [fns(-)]') 190 call mull_pop_fuk(geom, AO_bas_han, g_dens_ss, g_s) 191c 192 mu_n_men = 0.5d00*(eig_homo(1) + eig_homo(2)) 193 mu_n_cer = 0.5*(mu_n_mas + mu_n_men) 194 mu_s_mas = 0.5d00*(eig_lumo(1) - eig_homo(2)) 195 mu_s_men = 0.5d00*(eig_homo(1) - eig_lumo(2)) 196 mu_s_cer = 0.5*(mu_s_mas + mu_s_men) 197 if (me.eq.0.and.oprint_fukui)then 198 write(LuOut,*) ' ------------------------------------' 199 write(LuOut,*) ' mu_n(+) mu_n(-) mu_n(0)' 200 write(LuOut,'(3f11.4)') mu_n_mas, mu_n_men, mu_n_cer 201 write(LuOut,*) ' ------------------------------------' 202 write(LuOut,*) ' mu_s(+) mu_s(-) mu_s(0)' 203 write(LuOut,'(3f11.4)') mu_s_mas, mu_s_men,mu_s_cer 204 write(LuOut,*) ' ------------------------------------' 205 endif 206c 207cc Evaluating Coulomb integrals for HOMO, LUMO and differences 208c 209 kfac(1) = 0.d00 210 jfac(1) = 1.0d0 211 jfac(2) = 1.0d0 212 kfac(2) = 0.0d0 213 call ga_zero(g_orb) 214 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 215 & tol2e, oskel, g_dens_HOMO(1), g_orb, .false.) 216 int_HaHa = ga_ddot(g_dens_HOMO(1), g_orb) 217c 218 int_HaLb = ga_ddot(g_dens_LUMO(2), g_orb) 219c 220 call ga_zero(g_orb) 221 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 222 & tol2e, oskel, g_dens_HOMO(2), g_orb, .false.) 223 int_HbHb = ga_ddot(g_dens_HOMO(2), g_orb) 224c 225 int_LaHb = ga_ddot(g_dens_LUMO(1), g_orb) 226c 227 call ga_zero(g_orb) 228 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 229 & tol2e, oskel, g_dens_LUMO(1), g_orb, .false.) 230 int_LaLa = ga_ddot(g_dens_LUMO(1), g_orb) 231c 232 call ga_zero(g_orb) 233 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 234 & tol2e, oskel, g_dens_LUMO(2), g_orb, .false.) 235 int_LbLb = ga_ddot(g_dens_LUMO(2), g_orb) 236c 237c 238cc Evaluating exchange-correlation integrals 239c 240 g_dens(1) = ga_create_atom_blocked(geom, AO_bas_han, 241 & 'ga_dens(1)') 242 g_vxc(1) = ga_create_atom_blocked(geom, AO_bas_han, 243 & 'g_vxc(1)') 244 call ga_zero(g_dens(1)) 245 call ga_zero(g_vxc(1)) 246 if (ipol.eq.2) then 247 g_dens(2) = ga_create_atom_blocked(geom, AO_bas_han, 248 & 'ga_dens(2)') 249 g_vxc(2) = ga_create_atom_blocked(geom, AO_bas_han, 250 & 'g_vxc(2)') 251 call ga_zero(g_dens(2)) 252 call ga_zero(g_vxc(2)) 253 endif 254 do isp=1,ipol 255 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, 256 & noc(isp), 2d0/dble(ipol), g_movecs(isp), 257 & g_movecs(isp), 0.0d00, g_dens(isp)) 258 enddo 259 Exc(1) = 0.0d00 260 Exc(2) = 0.0d00 261 Ecoul = 0.0d00 262 call xc_pot(rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens, 263 & g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n, 264 & dbl_mb(irdens_atom), 265 & int_mb(icetobfr), natoms) 266 Exc_zero = Exc(1) 267 call ga_zero(g_orb) 268 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 269 & tol2e, oskel, g_dens_HOMO(1), g_orb, .false.) 270 do isp = 1, ipol 271 int_vxc_H(isp) = ga_ddot(g_dens_HOMO(isp), g_vxc(isp)) 272 int_vxc_L(isp) = ga_ddot(g_dens_LUMO(isp), g_vxc(isp)) 273 end do 274 if(ipol .eq.1) then 275 int_vxc_H(2) = int_vxc_H(1) 276 int_vxc_L(2) = int_vxc_L(1) 277 endif 278c 279cc Approximating the ionization potential 280c 281 noc_aux_1 = noc(1) 282 noc_aux_2 = noc(2) 283C Think about it. How can you excite only one electron and still be RKS? 284C Code produces rubbish with RKS - compare to UKS/ROKS and see 285C This code assumes once closed shell, then always closed shell 286 if(ipol .eq.1) then 287 noc_aux_2 = noc_aux_1 288 call errquit('fukui: UKS or ROKS required. No RKS', 0, 0) 289 end if 290 291 noc_test = noc_aux_1 - noc_aux_2 292 293 if (noc_test.eq.0) then 294 noc(2) = noc_aux_2 - 1 295 e_orbital = -eig_homo(2) 296 e_coul = 0.5d00*int_HbHb 297 e_xc = int_vxc_H(2) 298 else 299 noc(1) = noc_aux_1 - 1 300 e_orbital = -eig_homo(1) 301 e_coul = 0.5d00*int_HaHa 302 e_xc = int_vxc_H(1) 303 end if 304 ion_pot = e_orbital + e_coul + e_xc 305 if (.not.rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 306 & call errquit('fukui: rtdb_put of noc failed', 0, RTDB_ERR) 307 do isp=1,ipol 308 call ga_zero(g_dens(isp)) 309 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, 310 & noc(isp), 2d0/dble(ipol), g_movecs(isp), 311 & g_movecs(isp), 0.0d00, g_dens(isp)) 312 enddo 313 314 Exc(1) = 0.0d00 315 Exc(2) = 0.0d00 316 Ecoul = 0.0d00 317 do isp = 1, ipol 318 call ga_zero(g_vxc(isp)) 319 end do 320 call xc_pot(rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens, 321 & g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n, 322 & dbl_mb(irdens_atom), 323 & int_mb(icetobfr), natoms) 324 Exc_pert = Exc(1) 325 diff_Exc = Exc_pert - Exc_zero 326 ion_pot = ion_pot + diff_Exc 327 if (me.eq.0.and.oprint_fukui)then 328 write(LuOut,'(" Alpha = ",i3," ; Beta = ",i3)') 329 & noc(1), noc(2) 330 write(LuOut,*) ' Contributions in atomic units:' 331 write(LuOut,'(" Orbital Energy = ",f10.4)') e_orbital 332 write(LuOut,'(" Coulomb Integral = ",f10.4)') e_coul 333 write(LuOut,'(" XC Integral = ",f10.4)') e_xc 334 write(LuOut,'(" XC Diff. Energy = ",f10.4)') diff_Exc 335 write(LuOut,'(" Ionization potential = ",f7.4," a.u.")') 336 & ion_pot 337 write(LuOut,'(" = ",f7.2," eV")') 338 & ion_pot*27.211 339 end if 340c 341cc Approximating the electron affinity 342c 343 if (noc_test.eq.0) then 344 noc(1) = noc_aux_1 + 1 345 noc(2) = noc_aux_2 346 e_orbital = -eig_lumo(1) 347 e_coul = - 0.5d00*int_LaLa 348 e_xc = int_vxc_L(1) 349 else 350 noc(1) = noc_aux_1 351 noc(2) = noc_aux_2 + 1 352 e_orbital = -eig_lumo(2) 353 e_coul = - 0.5d00*int_LbLb 354 e_xc = int_vxc_L(2) 355 end if 356 ele_afi = e_orbital + e_coul + e_xc 357 if (.not.rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 358 & call errquit('fukui: rtdb_put of noc failed', 0, RTDB_ERR) 359 do isp=1,ipol 360 call ga_zero(g_dens(isp)) 361 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, 362 & noc(isp), 2d0/dble(ipol), g_movecs(isp), 363 & g_movecs(isp), 0.0d00, g_dens(isp)) 364 enddo 365 366 Exc(1) = 0.0d00 367 Exc(2) = 0.0d00 368 Ecoul = 0.0d00 369 do isp = 1, ipol 370 call ga_zero(g_vxc(isp)) 371 end do 372 call xc_pot(rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens, 373 & g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n, 374 & dbl_mb(irdens_atom), 375 & int_mb(icetobfr), natoms) 376 Exc_pert = Exc(1) 377 diff_Exc = Exc_pert - Exc_zero 378 ele_afi = ele_afi - diff_Exc 379 if (me.eq.0.and.oprint_fukui)then 380 write(LuOut,*) ' ------------------------------------' 381 write(LuOut,'(" Alpha = ",i3," ; Beta = ",i3)') 382 & noc(1), noc(2) 383 write(LuOut,*) ' Contributions in atomic units:' 384 write(LuOut,'(" Orbital Energy = ",f10.4)') e_orbital 385 write(LuOut,'(" Coulomb Integral = ",f10.4)') e_coul 386 write(LuOut,'(" XC Integral = ",f10.4)') e_xc 387 write(LuOut,'(" XC Diff. Energy = ",f10.4)') diff_Exc 388 write(LuOut,'(" Electron Affinity = ",f7.4," a.u.")') 389 & ele_afi 390 write(LuOut,'(" = ",f7.2," eV")') 391 & ele_afi*27.211 392 end if 393 if (me.eq.0.and.oprint_fukui)then 394 write(LuOut,*) ' ------------------------------------' 395 write(LuOut,'(" Electronegativity (I+A)/2 = ",f7.2," eV")') 396 & 0.5d00*(ion_pot + ele_afi)*27.211 397 write(LuOut,'(" Hardness (I-A) = ",f7.2," eV")') 398 & (ion_pot - ele_afi)*27.211 399 end if 400c 401cc Energy difference for high multiplicity 402c 403 if (noc_aux_2.gt.0) then 404 noc(1) = noc_aux_1 + 1 405 noc(2) = noc_aux_2 - 1 406 e_orbital = eig_lumo(1) - eig_homo(2) 407 e_coul = 0.5d00*(int_LaLa + int_HbHb) - int_LaHb 408 e_xc = int_vxc_H(2) - int_vxc_L(1) 409 high_mult = e_orbital + e_coul + e_xc 410 if (.not.rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 411 & call errquit('fukui: rtdb_put of noc failed', 0, RTDB_ERR) 412 do isp=1,ipol 413 call ga_zero(g_dens(isp)) 414 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, 415 & noc(isp), 2d0/dble(ipol), g_movecs(isp), 416 & g_movecs(isp), 0.0d00, g_dens(isp)) 417 enddo 418 Exc(1) = 0.0d00 419 Exc(2) = 0.0d00 420 Ecoul = 0.0d00 421 do isp = 1, ipol 422 call ga_zero(g_vxc(isp)) 423 end do 424 call xc_pot(rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens, 425 & g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n, 426 & dbl_mb(irdens_atom), 427 & int_mb(icetobfr), natoms) 428 Exc_pert = Exc(1) 429 diff_Exc = Exc_pert - Exc_zero 430 high_mult = high_mult + diff_Exc 431 if (me.eq.0.and.oprint_fukui)then 432 write(LuOut,*) ' ------------------------------------' 433 write(LuOut,'(" Alpha = ",i3," ; Beta = ",i3)') 434 & noc(1), noc(2) 435 write(LuOut,*) ' Contributions in atomic units:' 436 write(LuOut,'(" Orbital Energy = ",f10.4)') e_orbital 437 write(LuOut,'(" Coulomb Integrals = ",f10.4)') e_coul 438 write(LuOut,'(" XC Integrals = ",f10.4)') e_xc 439 write(LuOut,'(" XC Diff. Energy = ",f10.4)') diff_Exc 440 write(LuOut,'(" High Multiplicity = ",f7.4," a.u.")') 441 & high_mult 442 write(LuOut,'(" = ",f7.2," eV")') 443 & high_mult*27.211 444 end if 445 end if 446c 447cc Energy difference for low multiplicity 448c 449 if (noc_test.ge.2) then 450 noc(1) = noc_aux_1 - 1 451 noc(2) = noc_aux_2 + 1 452 low_mult = eig_lumo(2) - eig_homo(1) + 453 & 0.5d00*(int_LbLb + int_HaHa) - int_HaLb + 454 & int_vxc_H(1) - int_vxc_L(2) 455 e_orbital = eig_lumo(2) - eig_homo(1) 456 e_coul = 0.5d00*(int_LbLb + int_HaHa) - int_HaLb 457 e_xc = int_vxc_H(1) - int_vxc_L(2) 458 low_mult = e_orbital + e_coul + e_xc 459 if (.not.rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 460 & call errquit('fukui: rtdb_put of noc failed', 0, RTDB_ERR) 461 do isp=1,ipol 462 call ga_zero(g_dens(isp)) 463 call ga_dgemm('N', 'T', nbf_ao, nbf_ao, 464 & noc(isp), 2d0/dble(ipol), g_movecs(isp), 465 & g_movecs(isp), 0.0d00, g_dens(isp)) 466 enddo 467 468 Exc(1) = 0.0d00 469 Exc(2) = 0.0d00 470 Ecoul = 0.0d00 471 do isp = 1, ipol 472 call ga_zero(g_vxc(isp)) 473 end do 474 call xc_pot(rtdb, Exc, Ecoul,nExc, iVxc_opt, g_xcinv, g_dens, 475 & g_vxc, IOLGC, g_wght, g_xyz, g_nq,wght_GA, rho_n, 476 & dbl_mb(irdens_atom), 477 & int_mb(icetobfr), natoms) 478 Exc_pert = Exc(1) 479 diff_Exc = Exc_pert - Exc_zero 480 low_mult = low_mult + diff_Exc 481 if (me.eq.0.and.oprint_fukui)then 482 write(LuOut,*) ' ------------------------------------' 483 write(LuOut,'(" Alpha = ",i3," ; Beta = ",i3)') 484 & noc(1), noc(2) 485 write(LuOut,*) ' Contributions:' 486 write(LuOut,'(" Orbital Energy = ",f10.4)') e_orbital 487 write(LuOut,'(" Coulomb Integrals = ",f10.4)') e_coul 488 write(LuOut,'(" XC Integrals = ",f10.4)') e_xc 489 write(LuOut,'(" XC Diff. Energy = ",f10.4)') diff_Exc 490 write(LuOut,'(" Low Multiplicity = ",f7.4," a.u.")') 491 & low_mult 492 write(LuOut,'(" = ",f7.2," eV")') 493 & low_mult*27.211 494 write(LuOut,*) ' ------------------------------------' 495 end if 496 end if 497cc 498c 499 noc(1) = noc_aux_1 500 noc(2) = noc_aux_2 501 if (.not.rtdb_put(rtdb, 'dft:noc', mt_int, 2, noc)) 502 & call errquit('fukui: rtdb_put of noc failed', 0, RTDB_ERR) 503 if (.not. ga_destroy(g_dens(1))) call errquit 504 & ('fukui: could not destroy g_dens(1)',0, GA_ERR) 505 if (.not. ga_destroy(g_vxc(1))) call errquit 506 & ('fukui: could not destroy g_vxc(1)',0, GA_ERR) 507 if (ipol.eq.2) then 508 if (.not. ga_destroy(g_dens(2))) call errquit 509 & ('fukui: could not destroy g_dens(2)',0, GA_ERR) 510 if (.not. ga_destroy(g_vxc(2))) call errquit 511 & ('fukui: could not destroy g_vxc(2)',0, GA_ERR) 512 end if 513c 514 if (.not. ga_destroy(g_orb)) call errquit 515 & ('fukui: could not destroy g_orb', 0, GA_ERR) 516 if (.not. ga_destroy(g_dens_HOMO(1))) call errquit 517 & ('fukui: could not destroy g_dens_HOMO', 0, GA_ERR) 518 if (.not. ga_destroy(g_dens_LUMO(1))) call errquit 519 & ('fukui: could not destroy g_dens_LUMO', 0, GA_ERR) 520 if (.not. ga_destroy(g_dens_HOMO(2))) call errquit 521 & ('fukui: could not destroy g_dens_HOMO_bet', 0, GA_ERR) 522 if (.not. ga_destroy(g_dens_LUMO(2))) call errquit 523 & ('fukui: could not destroy g_dens_LUMO_bet', 0, GA_ERR) 524 if (.not. ga_destroy(g_dens_ss)) call errquit 525 & ('fukui: could not destroy g_dens_ss', 0, GA_ERR) 526 if (.not. ga_destroy(g_s)) call errquit 527 & ('fukui: could not destroy g_s', 0, GA_ERR) 528 return 529 end 530 531 Subroutine mull_pop_fuk( geom, basis, iga_dens,iga_s) 532 533C$Id$ 534 Implicit none 535#include "errquit.fh" 536 integer geom,basis 537 integer iga_s ! overlap GA handle 538 integer iga_dens ! dens. mat GA handle 539 integer iga_ps ! product GA handle 540 541 integer natoms,nshells 542 integer lPSmat,iPSmat,lqatom,iqatom,lqshell,iqshell 543 integer iatom,ilo,ihi,nbf,max_at_bf2 544c 545 integer ga_create_atom_blocked 546 external ga_create_atom_blocked 547 logical status 548 549#include "bas.fh" 550#include "geom.fh" 551#include "global.fh" 552#include "cscfps.fh" 553#include "mafdecls.fh" 554 555 if (oscfps) call pstat_on(ps_mull) 556 557 558c***************************************************************************** 559 560c 561 if(.not.geom_ncent(geom, natoms)) 562 & call errquit(' exiting in mull_pop',0, GEOM_ERR) 563 if( .not. bas_numcont(basis,nshells) ) 564 & call errquit(' exiting in mull_pop',1, BASIS_ERR) 565 if ( .not. bas_numbf(basis,nbf) ) 566 & call errquit(' exiting in mull_op',1, BASIS_ERR) 567 max_at_bf2 = 0 568 do iatom = 1, natoms 569 if (.not. bas_ce2bfr(basis, iatom, ilo, ihi)) 570 $ call errquit('mul_pop: bas_ce2bfr failed', iatom, 571 & BASIS_ERR) 572 max_at_bf2 = max(max_at_bf2, ihi-ilo+1) 573 enddo 574 max_at_bf2 = max_at_bf2*max_at_bf2 575 576 if(.not.MA_Push_Get(mt_dbl,max_at_bf2,'PS',lPSmat,iPSmat)) 577 & call errquit(' exiting in mull_pop: insuff stack',21, MA_ERR) 578 if(.not.MA_Push_Get(mt_dbl,natoms,'q atom',lqatom,iqatom)) 579 & call errquit(' exiting in mull_pop: insuff stack',22, MA_ERR) 580 if(.not.MA_Push_Get(mt_dbl,nshells,'q shell',lqshell,iqshell)) 581 & call errquit(' exiting in mull_pop: insuff stack',3, MA_ERR) 582 583 iga_PS=ga_create_atom_blocked(geom, basis, 'PS product') 584 585 call ga_dgemm('N','N',nbf,nbf,nbf,1.d0, 586 & iga_dens,iga_s,0.d0,iga_PS) 587 call mull_calc_fuk(basis,natoms, nshells,max_at_bf2,iga_PS, 588 & dbl_mb(iqatom),dbl_mb(iqshell),dbl_mb(iPSmat)) 589 590 call ga_sync 591 592 if(.not.MA_Pop_Stack(lqshell)) 593 & call errquit(' exiting in mull_pop',33, MA_ERR) 594 if(.not.MA_Pop_Stack(lqatom)) 595 & call errquit(' exiting in mull_pop',34, MA_ERR) 596 if(.not.MA_Pop_Stack(lPSmat)) 597 & call errquit(' exiting in mull_pop',35, MA_ERR) 598 status= ga_destroy(iga_PS) 599c 600 if (oscfps) call pstat_off(ps_mull) 601c 602 return 603 end 604c 605c 606c 607 Subroutine mull_calc_fuk(basis, natoms, nshells,max_at_bf2,iga_PS, 608 & qatom,qshell,PSmat) 609 610 Implicit none 611#include "errquit.fh" 612#include "geom.fh" 613#include "bas.fh" 614#include "global.fh" 615#include "tcgmsg.fh" 616#include "stdio.fh" 617#include "msgids.fh" 618#include "mafdecls.fh" 619#include "inp.fh" 620 integer basis 621 integer natoms,nshells 622 integer iga_PS ! product GA handle 623 integer ifirst,ilast,nbf_at,max_at_bf2 624 integer ish1,ish2,ish,nn,iat,mu 625 integer me,nproc, geom 626 double precision psmu, coord(3), qnuc 627 double precision qatom(natoms),qshell(nshells),PSmat(max_at_bf2) 628 character*2 symbol 629 character*16 tag 630 character*32 element 631 integer atn 632c 633 me=ga_nodeid() 634 nproc=ga_nnodes() 635 636 call dfill(natoms,0.D0,qatom,1) 637 call dfill(nshells,0.D0,qshell,1) 638 639 if (.not. bas_geom(basis, geom)) call errquit 640 $ ('mull_pop: bas_geom failed',basis, BASIS_ERR) 641 642 do iat=me+1,natoms,nproc 643 if (.not.bas_ce2cnr(basis,iat,ish1,ish2)) 644 & call errquit(' exiting in mull_pop',4, BASIS_ERR) 645 call get_atom_block(iga_PS, basis, 646 $ iat, iat, PSmat, nbf_at, nbf_at) 647 mu=0 648 do ish=ish1,ish2 649 if (.not. bas_cn2bfr(basis,ish,ifirst,ilast)) 650 & call errquit(' exiting in mull_pop.',5, BASIS_ERR) 651 do nn=ifirst,ilast 652 mu=mu+1 653 psmu=PSmat((mu-1)*nbf_at+mu) 654 qshell(ish)=qshell(ish)+psmu 655 enddo 656 qatom(iat)=qatom(iat)+qshell(ish) 657 enddo 658 enddo 659 call ga_sync 660 call ga_dgop(Msg_Mull1,qatom,natoms,'+') 661 call ga_dgop(Msg_Mull2,qshell,nshells,'+') 662 if(me.eq.0) then 663 write(LuOut,1) 664 1 format(/' Atom Condensed Fukui ') 665 write(luout,11) 666 11 format( ' ----------- ----------------') 667 do iat=1,natoms 668 if (.not.bas_ce2cnr(basis,iat,ish1,ish2)) 669 & call errquit(' exiting in mull_pop',4, BASIS_ERR) 670c 671 if (.not. geom_cent_get(geom, iat, tag, coord, qnuc)) 672 $ call errquit('mull_pop: geom_cent_tag failed',0, 673 & GEOM_ERR) 674c 675 if (.not. geom_tag_to_element(tag, symbol, element, atn)) 676 $ symbol = 'X' 677 if (ish2.ge.ish1) then 678 write(LuOut,2) iat,symbol,nint(qnuc),qatom(iat) 679 2 format(1x,i4,1x,a2,i4,1x,f10.4) 680 endif 681 enddo 682 endif 683c 684 call ga_sync 685 686 return 687 end 688 689 Subroutine xc_pot(rtdb, Exc, ecoul,nExc, iVxc_opt, g_xcinv, 690 & g_dens, g_vxc, IOLGC, g_wght, g_xyz,g_nq, 691 & wght_GA, rho_n, rdens_atom, 692 & cetobfr, natoms) 693c 694C$Id$ 695c 696 implicit none 697#include "errquit.fh" 698#include "stdio.fh" 699c 700 integer nExc 701 integer iVxc_opt 702 integer g_xcinv, g_dens(2), g_vxc(4), g_wght, g_xyz,g_nq 703 integer natoms 704 logical IOLGC, wght_GA 705 integer rtdb 706c 707#include "mafdecls.fh" 708#include "rtdb.fh" 709#include "bas.fh" 710#include "global.fh" 711#include "tcgmsg.fh" 712#include "cdft.fh" 713#include "util.fh" 714#include "sym.fh" 715c 716 integer cetobfr(2,natoms) 717 double precision rho_n, rdens_atom(ipol*natoms*natoms) 718 double precision jfac(4),kfac(4) 719 integer g_jk(4), g_d(4) 720c 721 integer ga_create_atom_blocked 722 external ga_create_atom_blocked 723c 724c--> XC Energy 725c 726 double precision Exc(2) 727 double precision ecoul ! [output] 728c 729c This driver routine solves for the XC energy and potential (Vxc) via 730c numerical quadrature methods. The results are obtained either by direct 731c numerical integration or by means of a LSQ fit of the Vxc to a set of 732c Gaussian functions. This fitted function can be used to evaluate Vxc 733c via a summation of a series of 3-center overlap integrals (3OIs). The 734c algorithms are formulated in terms of matrix products. See subsequent 735c subroutines for further explanation. 736c 737c XC Energy and Potential Index Key, Vxc(pq,i) 738c 739c Value of | Definition of index "i" 740c ipol nExc | 1 2 3 4 741c -------------------------------------------------- 742c 1 1 | Vxc 743c 2 1 | Vxc^up Vxc^dw 744c 1 2 | Vxc 745c 2 2 | Vxc^up Vxc^dw 746c 747c nTcols = ipol 748c 749 integer me,nproc,i,nTrows,nTcols 750 integer lTmat,iTmat,g_oep 751 double precision zero,one,onem 752 logical oprint_intermediate_xc, oprint_time, grid_on_file 753 parameter(zero=0.d0,one=1.d0,onem=-1.d0) 754 double precision tol2e,tot 755c****************************************************************************** 756c 757c Compute the matrix elements for the XC potential and energy. 758c 759 oprint_intermediate_xc = util_print('intermediate XC matrix', 760 $ print_debug) 761 oprint_time = util_print('dft timings', print_high) 762 Exc(1)=0.d0 763 Exc(2)=0.d0 764 me=ga_nodeid() 765 nproc=ga_nnodes() 766c 767 if (oprint_intermediate_xc)then 768 write(*,*)' rtdb, Exc, nExc, iVxc_opt, g_xcinv: ', 769 & rtdb, Exc, nExc, iVxc_opt, g_xcinv 770 write(*,*)' g_dens(1),g_vxc(1),IOLGC,g_wght,g_xyz,wght_GA: ', 771 & g_dens(1),g_vxc(1),IOLGC,g_wght,g_xyz,wght_GA 772 write(*,*)' Fock XC matrix entering xc_getv: ' 773 call ga_print(g_vxc(1)) 774 if(ipol.eq.2)call ga_print(g_vxc(2)) 775 endif 776c 777 if (abs(xfac(1)).gt.1e-8 .or. (.not. CDFIT))then 778c 779c Compute the exact exchange potential (as in Hartree-Fock calculations). 780c 781 tol2e=10.d0**(-itol2e) 782 call ga_sync 783 if (oprint_time)call dft_tstamp(' Before call to fock_2e. ') 784 if (ipol.eq.1) then 785 kfac(1) = -0.5d0*xfac(1) 786 jfac(1)=0.0d0 787 if (.not. CDFIT) then 788 jfac(2) = 1.0d0 789 kfac(2) = 0d0 790 g_vxc(2) = ga_create_atom_blocked(geom,ao_bas_han,'jk') 791 g_dens(2)=g_dens(1) 792 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 793 & tol2e, oskel, g_dens(1), g_vxc(1), .false.) 794 Exc(1) = Exc(1)+0.5d0*ga_ddot(g_dens(1),g_vxc(1)) 795 call ga_zero(g_vxc(2)) 796 ecoul = 0.5d0*ga_ddot(g_dens(1),g_vxc(2)) 797 call ga_dadd(1d0,g_vxc(1),1d0,g_vxc(2),g_vxc(1)) 798 if (.not. ga_destroy(g_vxc(2))) call errquit 799 $ ('xc_getv: ga corrupt?',0, GA_ERR) 800 else 801 call fock_2e(geom, AO_bas_han, 1, jfac, kfac, 802 & tol2e, oskel, g_dens(1), g_vxc(1), .false.) 803 Exc(1) = Exc(1)+0.5d0*ga_ddot(g_dens(1),g_vxc(1)) 804 endif 805 else 806 if (CDFIT) then 807 jfac(1)=0.d0 808 jfac(2)=0.d0 809 kfac(1)=-1.0d0*xfac(1) 810 kfac(2)=-1.0d0*xfac(1) 811 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 812 & tol2e, oskel, g_dens, g_vxc, .false.) 813 Exc(1) = Exc(1)+0.5d0*(ga_ddot(g_dens(1),g_vxc(1)) + 814 & ga_ddot(g_dens(2),g_vxc(2))) 815 else 816 jfac(1) = 1.0d0 817 jfac(2) = 0.0d0 818 jfac(3) = 1.0d0 819 jfac(4) = 0.0d0 820 kfac(1) = 0.0d0 821 kfac(2) = 1.0d0 822 kfac(3) = 0.0d0 823 kfac(4) = 1.0d0 824 g_jk(1) = g_vxc(1) ! This assignment is assumed 825 g_jk(2) = g_vxc(2) 826 g_jk(3) = ga_create_atom_blocked(geom, ao_bas_han, 'jk') 827 g_jk(4) = ga_create_atom_blocked(geom, ao_bas_han, 'jk') 828 call ga_zero(g_jk(3)) 829 call ga_zero(g_jk(4)) 830 g_d(1) = g_dens(1) 831 g_d(2) = g_dens(1) 832 g_d(3) = g_dens(2) 833 g_d(4) = g_dens(2) 834 call fock_2e(geom, AO_bas_han, 4, jfac, kfac, 835 & tol2e, oskel, g_d(1), g_jk(1),.false. ) 836 call ga_zero(g_jk(1)) 837 call ga_zero(g_jk(3)) 838 ecoul = 0.5d0*( ! Alpha coulomb energy 839 $ ga_ddot(g_dens(1),g_jk(1)) + 840 $ ga_ddot(g_dens(1),g_jk(3))) 841 ecoul = ecoul + 0.5d0*( ! Beta coulomb energy 842 $ ga_ddot(g_dens(2),g_jk(1)) + 843 $ ga_ddot(g_dens(2),g_jk(3))) 844 exc(1) = exc(1) - xfac(1)*0.5d0*( ! All exchange energy 845 $ ga_ddot(g_dens(1),g_jk(2)) + 846 $ ga_ddot(g_dens(2),g_jk(4))) 847 call ga_dadd(1.0d0, g_jk(1), 1.0d0, g_jk(3), g_jk(1)) 848 call ga_copy(g_jk(1), g_jk(3)) 849 call ga_dadd(1.0d0, g_jk(1), -xfac(1), g_jk(2), 850 $ g_jk(1)) 851 call ga_dadd(1.0d0, g_jk(3), -xfac(1), g_jk(4), 852 $ g_jk(2)) 853 if (.not. ga_destroy(g_jk(3))) call errquit 854 $ ('xc_getv: ga corrupt?',0, GA_ERR) 855 if (.not. ga_destroy(g_jk(4))) call errquit 856 $ ('xc_getv: ga corrupt?',1, GA_ERR) 857 endif 858 endif 859 if (oprint_time)call dft_tstamp(' After call to fock_2e. ') 860 call ga_sync 861c 862c Symmetrize Vxc? 863c 864c if (oskel)then 865c call sym_symmetrize(geom, AO_bas_han, .false., g_vxc(1)) 866c if (ipol.gt.1)then 867c call sym_symmetrize(geom, AO_bas_han, .false., 868c & g_vxc(2)) 869c endif 870c endif 871c 872c Compute the exact exchange energy. 873c 874 endif 875c 876 tot=-xfac(1) 877 do i=1,numfunc 878 tot=tot+xfac(i)+cfac(i) 879 enddo 880c 881 if (.not. rtdb_get(rtdb, 'dft:grid_on_file', mt_log, 1, 882 & grid_on_file))then 883 grid_on_file = .false. 884 endif 885 if (abs(tot).gt.1e-8) then 886 if(xcfit) then 887 if (.not. bas_numbf(XC_bas_han,nbf_xc) )then 888 call errquit('Exiting in getvxc.',1, BASIS_ERR) 889 endif 890 nTrows = nbf_xc 891 nTcols = ipol 892c 893c Allocate scratch space for the "T" matrix. 894c 895 if (.not.ma_push_get(MT_Dbl,nTrows*nTcols,'Tmat',lTmat, 896 & iTmat))call errquit('xc_getv: cannot allocate Tmat',0, 897 & MA_ERR) 898 call dfill(nTrows*nTcols,0.D0,dbl_mb(iTmat),1) 899 endif 900 call grid_quadv0(rtdb, g_dens, g_vxc, nexc,rho_n, Exc, 901 . dbl_mb(itmat)) 902 if(xcfit) then 903c 904c symmetrize the "T" vector 905c 906 if (oskel)then 907 call sym_vec_symmetrize(geom, xc_bas_han, Dbl_MB(iTmat)) 908 if (ipol.gt.1)then 909 call sym_vec_symmetrize(geom, xc_bas_han, 910 & Dbl_MB(iTmat+nbf_xc)) 911 endif 912 endif 913c 914 call xc_fitv(rtdb,Dbl_MB(iTmat), nTrows, nTcols, 915 & g_vxc, g_xcinv, g_oep, IOLGC) 916 if (oprint_time)call dft_tstamp(' After call to xc_fitv. ') 917 if (.not.ma_pop_stack(lTmat)) 918 & call errquit('xc_getv: cannot pop stack',0, MA_ERR) 919c 920 endif 921 endif 922c 923 if (oprint_intermediate_xc)then 924 write(*,*)' Fock XC matrix leaving xc_pot: ' 925 call ga_print(g_vxc(1)) 926 if(ipol.eq.2)call ga_print(g_vxc(2)) 927 endif 928c 929c 930 return 931 end 932 933