1 Subroutine xc_sic(rtdb, nExc, iVxc_opt, 2 & g_dens, g_vxc, g_wght, g_xyz,g_nq, 3 & wght_GA, rdens_atom, 4 & cetobfr, natoms, 5 & g_movecs, totsic, i_degen, n_levels) 6c 7C$Id$ 8c 9 implicit none 10#include "errquit.fh" 11c 12 integer iVxc_opt, nExc, g_dens(*), g_vxc(4), g_wght, 13 & g_xyz,g_nq 14c integer noc(2) 15 integer ik, ij, isp, n_states, naux_ipol, size_mat, 16c & g_sic_dens(2),g_orb,g_vxc_orb(4),g_movecs(2), 17 & g_orb,g_vxc_orb(4),g_movecs(2), 18 & g_matm, g_vectu, g_coul_orb, g_vxc_sig(2), 19 & l_vec_aux, i_vec_aux, j, 20 & l_vec_u, i_vec_u, l_consti, i_consti, 21 & i_degen(2), i_gauxsig, l_gauxsig, 22 & n_levels(2), n_orbitals, i_level, i_orb, j_level, j_orb, 23 & m_orbitals, elem_i, elem_j 24 double precision Exc_orb(2), tot_sic_xc, tot_sic_coul, ecoul_orb, 25 & u_array, res_int,const_ci, totsic 26 integer natoms 27 logical wght_GA 28 integer rtdb 29c 30#include "mafdecls.fh" 31#include "rtdb.fh" 32#include "bas.fh" 33#include "global.fh" 34#include "tcgmsg.fh" 35#include "cdft.fh" 36#include "util.fh" 37#include "sym.fh" 38#include "stdio.fh" 39c 40 integer cetobfr(2,natoms) 41 double precision rho_n, rdens_atom(ipol*natoms*natoms) 42 double precision jfac(4),kfac(4) 43c 44 integer ga_create_atom_blocked 45 external ga_create_atom_blocked 46c 47 integer me,nproc,ii 48c integer me,nproc,i,nTrows,nTcols 49c integer lTmat,iTmat 50 double precision zero,one,onem 51 logical oprint_intermediate_xc, oprint_time, oprint_sic 52 parameter(zero=0.d0,one=1.d0,onem=-1.d0) 53 double precision tol2e,edodumm 54c****************************************************************************** 55c 56c Compute the matrix elements for the XC potential of the OEP and the 57c energy of the SIC approximation 58c 59 oprint_intermediate_xc = util_print('intermediate XC matrix', 60 $ print_debug) 61 oprint_time = util_print('dft timings', print_high) 62 oprint_sic = util_print('SIC information', print_high) 63c 64 me=ga_nodeid() 65 nproc=ga_nnodes() 66c 67c tol2e def 68c 69 tol2e = 10.d0**(-itol2e) 70c 71 if (iVxc_opt.eq.0) then 72c 73c Option A: Compute via direct numerical quadrature. 74c 75 if (ipol.eq.2) then 76 if (noc(2).eq.0) call ga_zero(g_vxc(2)) 77 end if 78c 79 if (me.eq.0.and.oprint_sic) write(LuOut,*) 80 & ' Starting SIC by orbital..' 81 g_vxc_orb(1) = ga_create_atom_blocked(geom, ao_bas_han, 82 & 'Vxc_orb_alpha') 83 g_vxc_orb(2) = ga_create_atom_blocked(geom, ao_bas_han, 84 & 'Vxc_orb_beta') 85 g_orb = ga_create_atom_blocked(geom, ao_bas_han, 'Orbsic') 86 g_sic_dens(1) = ga_create_atom_blocked(geom, ao_bas_han, 87 & 'Densic_alpha') 88 g_sic_dens(2) = ga_create_atom_blocked(geom, ao_bas_han, 89 & 'Densic_beta') 90 g_coul_orb= ga_create_atom_blocked(geom, ao_bas_han, 91 & 'V_coul_orb') 92 g_vxc_sig(1) = ga_create_atom_blocked(geom, ao_bas_han, 93 & 'Vxc_sig_alpha') 94 if (ipol.eq.2) then 95 g_vxc_sig(2) = ga_create_atom_blocked(geom, ao_bas_han, 96 & 'Vxc_sig_beta') 97 end if 98c 99 tot_sic_xc = 0.d0 100 tot_sic_coul = 0.d0 101c 102 if (ipol.eq.2) then 103 naux_ipol = 2 104 sic_orb_occ = 1 105 else 106 naux_ipol = 1 107 sic_orb_occ = 2 108 end if 109c 110 do 300 isp = 1,ipol 111c 112 n_states = noc(isp) 113 size_mat = 0 114 if (noc(isp).ne.0) then 115 do j = 2, n_levels(isp) 116 n_orbitals = Int_MB(i_degen(isp) + j - 1) 117 do ik = 1, n_orbitals 118 size_mat = size_mat + 1 119 end do 120 end do 121 end if 122 call ga_zero(g_vxc_sig(isp)) 123 if (size_mat.gt.0) then 124 if (.not. ga_create(mt_dbl, size_mat,size_mat,'Mat_m', 125 & size_mat, size_mat, g_matm)) 126 & call errquit('dft_main0d: error creating g_matm',0, 127 & GA_ERR) 128 if (.not. ga_create(mt_dbl, size_mat, 1, 'Vet_u', 129 & size_mat, 1, g_vectu)) 130 & call errquit('dft_main0d: error creating Vet_u',0, 131 & GA_ERR) 132 if (.not.MA_Push_Get(MT_Dbl, size_mat, 'const_i', 133 & l_consti, i_consti)) 134 & call errquit('xc_sic: cannot allocate vec aux',0, MA_ERR) 135 end if 136 if (.not.MA_Push_Get(MT_Dbl, nbf_ao, 'vec aux', 137 & l_vec_aux, i_vec_aux)) 138 & call errquit('xc_sic: cannot allocate vec aux',0, MA_ERR) 139 if (size_mat.gt.0) then 140 if (.not.MA_Push_Get(MT_Dbl, size_mat, 'vec u', 141 & l_vec_u, i_vec_u)) 142 & call errquit('xc_sic: cannot allocate vec u',0, MA_ERR) 143 if (.not.ma_push_get(mt_int, size_mat, 'g_aux_sig', 144 & l_gauxsig, i_gauxsig)) 145 & call errquit('xc_sic:push_get failed', 0, MA_ERR) 146 end if 147c 148c Start loop of total states by spin 149c 150cc 151 ik = noc(isp) + 1 152 elem_i = size_mat + 1 153 do 100 i_level = 1,n_levels(isp) 154 n_orbitals = Int_MB(i_degen(isp) + i_level - 1) 155 do 90 i_orb = 1,n_orbitals 156 ik = ik - 1 157 if (me.eq.0.and.oprint_sic) write(LuOut,*) ' Orbital ',ik 158 call ga_zero(g_sic_dens(1)) 159 call ga_zero(g_sic_dens(2)) 160 call ga_zero(g_orb) 161 call ga_zero(g_vxc_orb(1)) 162 call ga_zero(g_vxc_orb(2)) 163 call ga_zero(g_coul_orb) 164 call ga_get(g_movecs(isp), 1, nbf_ao, ik, ik, 165 & Dbl_MB(i_vec_aux), nbf_ao) 166 call ga_put(g_orb, 1, nbf_ao, ik, ik, 167 & Dbl_MB(i_vec_aux), nbf_ao) 168 call ga_dgemm('n', 't', nbf_ao, nbf_ao, nbf_ao, 1.0d00, 169 & g_orb, g_orb, 0.d00, g_sic_dens(1)) 170c 171c g_sic_dens(1) is the orbital density of the state ik 172c 173 if (abs(xfac(1)).gt.1e-8 )then 174 if (me.eq.0.and.oprint_sic) write(LuOut,*) 175 & ' NO SIC in Coulomb term' 176 else 177 if (me.eq.0.and.oprint_sic) write(LuOut,*) 178 & ' SIC approximation in Coulomb term' 179 kfac(1) = 0.d00 180 jfac(1) = 1.0d0 181 jfac(2) = 1.0d0 182 kfac(2) = 0d0 183 tol2e=10.d0**(-itol2e) 184 call fock_2e(geom, AO_bas_han, 2, jfac, kfac, 185 & tol2e, oskel, g_sic_dens(1), g_vxc_orb(1),.false.) 186 ecoul_orb = 0.0d0 187 ecoul_orb = 0.5d0*ga_ddot(g_sic_dens(1), g_vxc_orb(1)) 188c 189c ecoul_orb is the orbital coulomb energy 190c 191 tot_sic_coul = tot_sic_coul - 192 & dble(sic_orb_occ)*ecoul_orb 193 call ga_copy(g_vxc_orb(1), g_coul_orb) 194 end if 195c 196 ipol = 2 197 call ga_zero(g_vxc_orb(1)) 198 Exc_orb(1) = 0.0d0 199 Exc_orb(2) = 0.0d0 200 if (oprint_time)call dft_tstamp 201 & ('Before call to grid_quadv0.') 202 sic_orb_spin = isp 203 sic_orb_index = 0 204 rho_n = 0.0d00 205 if(naux_ipol.eq.1) then 206 do ii=1,natoms*natoms 207 rdens_atom(natoms*natoms+ii)=rdens_atom(ii) 208 enddo 209 endif 210 call grid_quadv0(rtdb, g_sic_dens, g_vxc_orb, 211 , nexc,rho_n, Exc_orb,edodumm) 212c 213c Exc_orb(1) is the exchange-correlation energy by orbital 214c 215 tot_sic_xc = tot_sic_xc - dble(sic_orb_occ)*Exc_orb(1) 216c 217 if (i_level.ne.1) then 218 elem_i = elem_i - 1 219 call ga_zero(g_orb) 220 call ga_add(1.0d0, g_coul_orb, 1.0d0, g_vxc_orb(1), 221 & g_orb) 222c 223c g_orb contains the coulomb + exchange-correlation potential evaluated with 224c the ik-orbital 225c 226 u_array = -ga_ddot(g_sic_dens(1), g_orb) 227 Dbl_MB(i_vec_u + elem_i - 1) = u_array 228 end if 229 230 call ga_zero(g_vxc_orb(1)) 231 call ga_zero(g_vxc_orb(2)) 232 sic_orb_spin = isp 233 sic_orb_index = 1 234 rho_n = 0.0d00 235 call grid_quadv0(rtdb, g_dens, g_vxc_orb, 236 & nexc,rho_n, Exc_orb,edodumm) 237 sic_orb_index = 0 238 ij = noc(isp) - Int_MB(i_degen(isp)) + 1 239 elem_j = size_mat + 1 240 do 80 j_level = 2, n_levels(isp) 241 m_orbitals = Int_MB(i_degen(isp) + j_level - 1) 242 do 70 j_orb = 1, m_orbitals 243 if (i_level.ne.1) then 244 ij = ij -1 245 elem_j = elem_j - 1 246 call ga_zero(g_orb) 247 call ga_zero(g_sic_dens(2)) 248 call ga_get(g_movecs(isp), 1, nbf_ao, ij, ij, 249 & Dbl_MB(i_vec_aux) ,nbf_ao) 250 call ga_put(g_orb, 1, nbf_ao, ij, ij, 251 & Dbl_MB(i_vec_aux) ,nbf_ao) 252 call ga_dgemm('n', 't', nbf_ao, nbf_ao, nbf_ao, 253 & 1.0d00, g_orb, g_orb, 0.d00, 254 & g_sic_dens(2)) 255 res_int = -ga_ddot(g_sic_dens(2), g_vxc_orb(2)) 256 if (ij.eq.ik) res_int = 1.0d00 + res_int 257 call ga_fill_patch(g_matm, elem_i, elem_i, 258 & elem_j, elem_j, res_int) 259 end if 260 70 continue 261 80 continue 262c 263 if (i_level.ne.1) then 264 if (.not. ga_create(mt_dbl, nbf_ao, nbf_ao,'i_gauxsig', 265 & nbf_ao, nbf_ao, 266 & int_mb(i_gauxsig + elem_i - 1))) 267 & call errquit( 268 ' ' xc_sic: error creating i_gauxsig',0, GA_ERR) 269 call ga_copy(g_vxc_orb(2), 270 & int_mb(i_gauxsig + elem_i - 1)) 271 end if 272 273 call ga_add(-1.0d0, g_vxc_orb(1), 1.0d00, 274 & g_vxc_sig(isp), g_vxc_sig(isp)) 275 if (naux_ipol.eq.1) then 276 ipol=1 277 end if 278 if (oprint_time)call dft_tstamp 279 & (' After call to grid_quadv0.') 280 90 continue 281 100 continue 282c 283 ik = noc(isp) - Int_MB(i_degen(isp)) + 1 284 elem_i = size_mat + 1 285 do 200 i_level = 2, n_levels(isp) 286 n_orbitals = Int_MB(i_degen(isp) + i_level - 1) 287 do 190 i_orb = 1,n_orbitals 288 ik = ik - 1 289 elem_i = elem_i - 1 290 call ga_zero(g_orb) 291 call ga_zero(g_sic_dens(2)) 292 call ga_get(g_movecs(isp), 1, nbf_ao, ik, ik, 293 & Dbl_MB(i_vec_aux) ,nbf_ao) 294 call ga_put(g_orb, 1, nbf_ao, ik, ik, 295 & Dbl_MB(i_vec_aux), nbf_ao) 296 call ga_dgemm('n', 't', nbf_ao, nbf_ao, nbf_ao, 297 & 1.0d00, g_orb, g_orb, 0.d00, 298 & g_sic_dens(2)) 299 res_int = ga_ddot(g_sic_dens(2), g_vxc_sig(isp)) 300 Dbl_MB(i_consti + elem_i - 1) = res_int-Dbl_MB(i_vec_u + 301 & elem_i - 1) 302 190 continue 303 200 continue 304c 305 if (size_mat.gt.0) then 306 call ga_zero(g_vectu) 307 call ga_put(g_vectu, 1, size_mat, 1, 1, Dbl_MB(i_consti), 308 & size_mat) 309 call ga_lu_solve_seq('n', g_matm, g_vectu) 310 call ga_get(g_vectu, 1, size_mat, 1, 1, Dbl_MB(i_consti), 311 & size_mat) 312 call ga_zero(g_orb) 313 elem_i = size_mat + 1 314 do 400 i_level = 2, n_levels(isp) 315 n_orbitals = Int_MB(i_degen(isp) + i_level - 1) 316 do 390 i_orb = 1,n_orbitals 317 elem_i = elem_i - 1 318 const_ci= Dbl_MB(i_consti + elem_i - 1) 319 call ga_add(1.0d0, g_orb, const_ci, 320 & int_mb(i_gauxsig + elem_i - 1), g_orb) 321 390 continue 322 400 continue 323 call ga_add(1.0d0, g_orb, 1.0d00, g_vxc_sig(isp), 324 & g_vxc_sig(isp)) 325 do ik = 1, size_mat 326 if (.not. ga_destroy(int_mb(i_gauxsig + ik - 1))) 327 & call errquit ('xc_sic: could not destroy i_gauxsig', 328 & 0, GA_ERR) 329 end do 330 if (size_mat.gt.0) then 331 if (.not.ma_pop_stack(l_gauxsig)) 332 & call errquit('xc_sic: cannot pop stack',0, MA_ERR) 333 if (.not.ma_pop_stack(l_vec_u)) 334 & call errquit('xc_sic: cannot pop stack',0, MA_ERR) 335 end if 336 if (.not.ma_pop_stack(l_vec_aux)) 337 & call errquit('xc_sic: cannot pop stack',0, MA_ERR) 338c 339c g_vxc_sig is the exchange-correlation potential associate to the OEP 340c 341 end if 342 call ga_add(1.0d0, g_vxc_sig(isp), 1.0d00, g_vxc(isp), 343 & g_vxc(isp)) 344c 345 if (size_mat.gt.0) then 346 if (.not.ma_pop_stack(l_consti)) 347 & call errquit('xc_sic: cannot pop stack',0, MA_ERR) 348 if (.not. ga_destroy(g_matm)) call errquit 349 & ('xc_getv: could not destroy g_matm', 0, GA_ERR) 350 if (.not. ga_destroy(g_vectu)) call errquit 351 & ('xc_getv: could not destroy g_vectu', 0, GA_ERR) 352 end if 353 300 continue 354 totsic = tot_sic_xc + tot_sic_coul 355 sic_orb_occ = 0 356 if (me.eq.0.and.oprint_sic) write(LuOut,*) 357 & ' tot_sic_coul, tot_sic_xc, tot_sic: ',tot_sic_coul, 358 & tot_sic_xc, totsic 359c 360 if (.not. ga_destroy(g_vxc_orb(1))) call errquit 361 & ('xc_sic: could not destroy g_vxc_orb(1)', 0, GA_ERR) 362 if (.not. ga_destroy(g_vxc_orb(2))) call errquit 363 & ('xc_sic: could not destroy g_vxc_orb(2)', 0, GA_ERR) 364 if (.not. ga_destroy(g_orb)) call errquit 365 & ('xc_sic: could not destroy g_orb', 0, GA_ERR) 366 if (.not. ga_destroy(g_sic_dens(1))) call errquit 367 & ('xc_sic: could not destroy g_sic_dens(1)', 0, GA_ERR) 368 if (.not. ga_destroy(g_sic_dens(2))) call errquit 369 & ('xc_sic: could not destroy g_sic_dens(2)', 0, GA_ERR) 370 if (.not. ga_destroy(g_coul_orb)) call errquit 371 & ('xc_sic: could not destroy g_coul_orb', 0, GA_ERR) 372 if (.not. ga_destroy(g_vxc_sig(1))) call errquit 373 & ('xc_sic: could not destroy g_vxc_sig', 0, GA_ERR) 374 if (ipol.eq.2) then 375 if (.not. ga_destroy(g_vxc_sig(2))) call errquit 376 & ('xc_sic: could not destroy g_vxc_sig', 0, GA_ERR) 377 end if 378c 379 elseif (iVxc_opt.eq.1 )then 380 call errquit 381 & ('xc_sic: SIC + XC fitting not implemented', 0, INPUT_ERR) 382 endif 383c 384 if (oprint_intermediate_xc)then 385 write(*,*)' Fock XC matrix leaving xc_sic: ' 386 call ga_print(g_vxc(1)) 387 if(ipol.eq.2)call ga_print(g_vxc(2)) 388 endif 389c 390 return 391 end 392 393