1c $Id$ 2 3* ********************************* 4* * * 5* * nwpw_hartree_init * 6* * * 7* ********************************* 8 subroutine nwpw_hartree_init(nkatm, 9 > nprj,nbasis,psp_type,lmax, 10 > nprj_max0,l_prj,m_prj,b_prj, 11 < hartree_tag) 12 implicit none 13 integer nkatm 14 integer nprj(*),nbasis(*),psp_type(*),lmax(*) 15 integer nprj_max0 16 integer l_prj(nprj_max0,*),m_prj(nprj_max0,*),b_prj(nprj_max0,*) 17 integer hartree_tag 18 19 20#include "bafdecls.fh" 21#include "errquit.fh" 22#include "nwpw_hartree.fh" 23 24c **** local variables **** 25 logical ok 26 integer ii,ia,nsize 27 integer ic 28 integer l,m,li,mi,lj,mj,bi,bj,iprj,jprj,lm 29 integer li1,mi1,lj1,mj1,bi1,bj1,iprj1,jprj1 30 integer i_p,i_t 31 integer i_tlm,i_plm,i_rlm 32 integer matr_ptr 33 integer nb 34 real*8 tmp_theta,cs_theta,tmp_phi,tmp_gaunt 35 36c **** external functions **** 37 integer psi_data_get_chnk 38 external psi_data_get_chnk 39 real*8 rtheta_lm,drtheta_lm,nwpw_gaunt 40 external rtheta_lm,drtheta_lm,nwpw_gaunt 41 42 call nwpw_timing_start(4) 43 ok = .true. 44 45* ***** set up index arrays for nwpw_hartree_solve2 ***** 46 ok = BA_alloc_get(mt_int,nkatm,"nindx_hartree", 47 > nindx_hartree(2),nindx_hartree(1)) 48 ok = ok.and. 49 > BA_alloc_get(mt_int,nkatm,"shift_hartree", 50 > shift_hartree(2),shift_hartree(1)) 51 if (.not.ok) 52 >call errquit("nwpw_hartree_init:error allocating work arrays",0,0) 53 54 nsize = 0 55 do ia=1,nkatm 56 nb = nbasis(ia) 57 matr_ptr = psi_data_get_chnk(hartree_tag(ia)) 58 int_mb(shift_hartree(1)+ia-1) = nsize 59 do jprj = 1,nprj(ia) 60 lj = l_prj(jprj,ia) 61 mj = m_prj(jprj,ia) 62 bj = b_prj(jprj,ia) 63 do iprj = 1,nprj(ia) 64 li = l_prj(iprj,ia) 65 mi = m_prj(iprj,ia) 66 bi = b_prj(iprj,ia) 67 do jprj1 = 1,nprj(ia) 68 lj1 = l_prj(jprj1,ia) 69 mj1 = m_prj(jprj1,ia) 70 bj1 = b_prj(jprj1,ia) 71 do iprj1 = 1,nprj(ia) 72 li1 = l_prj(iprj1,ia) 73 mi1 = m_prj(iprj1,ia) 74 bi1 = b_prj(iprj1,ia) 75 do l=0,2*lmax(ia) 76 do m=-l,l 77 tmp_gaunt = 78 > nwpw_gaunt(.false.,l,m,li,mi,lj,mj) 79 > *nwpw_gaunt(.false.,l,m,li1,mi1,lj1,mj1) 80 > *dbl_mb(matr_ptr 81 > + (bi-1) 82 > + (bj-1)*nb 83 > + (bi1-1)*nb*nb 84 > + (bj1-1)*nb*nb*nb 85 > + l*nb*nb*nb*nb) 86 if (dabs(tmp_gaunt).gt.1.0d-11) then 87 nsize = nsize + 1 88 end if 89 end do 90 end do 91 92 end do 93 end do 94 95 end do 96 end do 97 int_mb(nindx_hartree(1)+ia-1)= nsize 98 > - int_mb(shift_hartree(1)+ia-1) 99 end do 100 101 ok = BA_alloc_get(mt_int,nsize,"iprj_hartree", 102 > iprj_hartree(2),iprj_hartree(1)) 103 ok = ok.and. 104 > BA_alloc_get(mt_int,nsize,"jprj_hartree", 105 > jprj_hartree(2),jprj_hartree(1)) 106 ok = ok.and. 107 > BA_alloc_get(mt_int,nsize,"iprj1_hartree", 108 > iprj1_hartree(2),iprj1_hartree(1)) 109 ok = ok.and. 110 > BA_alloc_get(mt_int,nsize,"jprj1_hartree", 111 > jprj1_hartree(2),jprj1_hartree(1)) 112 ok = ok.and. 113 > BA_alloc_get(mt_dbl,nsize,"coeff_hartree", 114 > coeff_hartree(2),coeff_hartree(1)) 115 if (.not.ok) 116 >call errquit("nwpw_hartree_init:error allocating work arrays",0,0) 117 118 119 nsize = 0 120 do ia=1,nkatm 121 nb = nbasis(ia) 122 matr_ptr = psi_data_get_chnk(hartree_tag(ia)) 123 do jprj = 1,nprj(ia) 124 lj = l_prj(jprj,ia) 125 mj = m_prj(jprj,ia) 126 bj = b_prj(jprj,ia) 127 do iprj = 1,nprj(ia) 128 li = l_prj(iprj,ia) 129 mi = m_prj(iprj,ia) 130 bi = b_prj(iprj,ia) 131 do jprj1 = 1,nprj(ia) 132 lj1 = l_prj(jprj1,ia) 133 mj1 = m_prj(jprj1,ia) 134 bj1 = b_prj(jprj1,ia) 135 do iprj1 = 1,nprj(ia) 136 li1 = l_prj(iprj1,ia) 137 mi1 = m_prj(iprj1,ia) 138 bi1 = b_prj(iprj1,ia) 139 do l=0,2*lmax(ia) 140 do m=-l,l 141 tmp_gaunt = 142 > nwpw_gaunt(.false.,l,m,li,mi,lj,mj) 143 > *nwpw_gaunt(.false.,l,m,li1,mi1,lj1,mj1) 144 > *dbl_mb(matr_ptr 145 > + (bi-1) 146 > + (bj-1)*nb 147 > + (bi1-1)*nb*nb 148 > + (bj1-1)*nb*nb*nb 149 > + l*nb*nb*nb*nb) 150 if (dabs(tmp_gaunt).gt.1.0d-11) then 151 int_mb(iprj_hartree(1)+nsize) = iprj 152 int_mb(jprj_hartree(1)+nsize) = jprj 153 int_mb(iprj1_hartree(1)+nsize) = iprj1 154 int_mb(jprj1_hartree(1)+nsize) = jprj1 155 nsize = nsize + 1 156 end if 157 end do 158 end do 159 160 end do 161 end do 162 163 end do 164 end do 165 end do 166 167 call nwpw_timing_end(4) 168 return 169 end 170 171* ******************************************** 172* * * 173* * nwpw_hartree_solve * 174* * * 175* ******************************************** 176 subroutine nwpw_hartree_solve(ia, 177 > ispin,ne,nprj,sw1,sw2) 178 implicit none 179 integer ii,ia 180 integer ispin,ne(2),nprj 181 real*8 sw1(ne(1)+ne(2),nprj) 182 real*8 sw2(ne(1)+ne(2),nprj) 183 184#include "bafdecls.fh" 185#include "errquit.fh" 186#include "nwpw_hartree.fh" 187 188 integer shift 189 190 call nwpw_timing_start(4) 191 call nwpw_timing_start(22) 192 193* **** hartree potential non-local matrix elements **** 194 shift = int_mb(shift_hartree(1)+ia-1) 195 call nwpw_hartree_sw1tosw2(ispin,ne,nprj,sw1,sw2, 196 > int_mb(nindx_hartree(1)+ia-1), 197 > int_mb(iprj_hartree(1)+shift), 198 > int_mb(jprj_hartree(1)+shift), 199 > int_mb(iprj1_hartree(1)+shift), 200 > int_mb(jprj1_hartree(1)+shift), 201 > dbl_mb(coeff_hartree(1)+shift)) 202 203 call nwpw_timing_end(4) 204 call nwpw_timing_end(22) 205 206 return 207 end 208 209 210* ***************************************** 211* * * 212* * nwpw_hartree_end * 213* * * 214* ***************************************** 215 subroutine nwpw_xc_end() 216 implicit none 217 218#include "bafdecls.fh" 219#include "errquit.fh" 220#include "nwpw_hartree.fh" 221 222 logical ok 223 224 call nwpw_timing_start(4) 225 ok = .true. 226 ok = ok.and.BA_free_heap(coeff_hartree(2)) 227 ok = ok.and.BA_free_heap(jprj_hartree(2)) 228 ok = ok.and.BA_free_heap(iprj_hartree(2)) 229 ok = ok.and.BA_free_heap(jprj1_hartree(2)) 230 ok = ok.and.BA_free_heap(iprj1_hartree(2)) 231 ok = ok.and.BA_free_heap(shift_hartree(2)) 232 ok = ok.and.BA_free_heap(nindx_hartree(2)) 233 234 if (.not.ok) 235 > call errquit("nwpw_hartree_end: error freeing heap",0,MA_ERR) 236 237 call nwpw_timing_end(4) 238 return 239 end 240 241c ********************************************* 242c * * 243c * nwpw_hartree_sw1tosw2 * 244c * * 245c ********************************************* 246 subroutine nwpw_hartree_sw1tosw2(ispin,ne,nprj,sw1,sw2, 247 > nindx, 248 > iprj_hartree, jprj_hartree, 249 > iprj1_hartree,jprj1_hartree, 250 > coeff_hartree) 251 implicit none 252 integer ispin,ne(2),nprj 253 real*8 sw1(ne(1)+ne(2),nprj) 254 real*8 sw2(ne(1)+ne(2),nprj) 255 integer n1dgrid,nbasis,lmax2 256 integer nindx,iprj_hartree(*),jprj_hartree(*) 257 integer iprj1_hartree(*),jprj1_hartree(*) 258 real*8 coeff_hartree(*) 259 260 integer n,i,iprj,jprj,iprj1,jprj1 261 real*8 coeff,w,scal 262 263* **** external functions **** 264 real*8 lattice_omega 265 external lattice_omega 266 267 call nwpw_timing_start(21) 268 scal = 1.0d0/dsqrt(lattice_omega()) 269 scal = scal/lattice_omega() 270 271* ***init to zero*** 272 do i=1,nindx 273 iprj = iprj_hartree(i) 274 jprj = jprj_hartree(i) 275 iprj1 = iprj1_hartree(i) 276 jprj1 = jprj1_hartree(i) 277 coeff = coeff_hartree(i)*scal 278 w = 0.0d0 279 do n=1,ne(1)+ne(2) 280 w = w + sw1(n,iprj1)*sw1(n,jprj1) 281 end do 282 do n=1,ne(1)+ne(2) 283 sw2(n,iprj) = sw2(n,iprj) + coeff*w*sw1(n,jprj) 284 end do 285 end do 286 call nwpw_timing_end(21) 287 return 288 end 289