1c 2c == compute the matrix elements for the zora corrections == 3 subroutine zora_getv_sf(rtdb, g_dens, 4 & g_zora_sf, g_zora_scale_sf, 5 & ofinite,zetanuc_arr, 6 & Knucl, nexc) 7c 8C$Id$ 9c Modified from zora_getv_sf() by FA 10-31-10 10c 11 implicit none 12c 13#include "rtdb.fh" 14#include "bas.fh" 15#include "cdft.fh" 16#include "errquit.fh" 17#include "mafdecls.fh" 18#include "global.fh" 19#include "geom.fh" 20#include "msgtypesf.h" 21#include "msgids.fh" 22#include "stdio.fh" 23#include "cgridfile.fh" 24#include "grid_cube.fh" 25#include "modelpotential.fh" 26c 27c == arguments == 28 integer rtdb 29 integer g_dens(2) 30 integer g_zora_sf(2) 31 integer g_zora_scale_sf(2) 32 logical ofinite 33 double precision zetanuc_arr(*)!for Gaussian Nuclear Model 34 logical Knucl 35 integer nexc 36c 37c == local variables == 38 integer i,j,nij 39 double precision rho_n 40 double precision tmat 41 double precision dummy(2) 42 integer ilo,ihi,jlo,jhi 43 integer iqsh, istep, nxyz, ncontrset 44 integer ixyz, lxyz, icharge, lcharge, itags, ltags 45 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 46 & ictr_buf,iqpts 47c 48 double precision rad,ke 49 integer lbas_cset_info, ibas_cset_info, 50 & lbas_cent_info, ibas_cent_info, 51 & ldocset, idocset, 52 & l_rchi_atom,i_rchi_atom, 53 & l_rq,i_rq,lniz, iniz, 54 & lchi_ao, ichi_ao, 55 & ldelchi_ao, idelchi_ao, 56 & lzora0, izora0, 57 & lzora0scal, izora0scal 58c 59 integer inntsize,ddblsize,ok 60c 61 logical dorepl 62 integer g_densrep(2),g_dens0(2),dorep_glob 63 logical docopy,dozero,dorepd,dorepon 64 logical util_mirrmat 65 external util_mirrmat 66c 67 logical grid_file_rewind 68 external grid_file_rewind 69c 70 logical int_normalize 71 external int_normalize 72c 73c model potential parameters 74 character*2 gelem(ncenters) 75 double precision gexpo(ncenters,50) 76 double precision gcoef(ncenters,50) 77c 78c == allocate memory == 79 if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao,'zora0',lzora0,izora0)) 80 & call errquit('zora_getv: zorasf',0, MA_ERR) 81 if (.not.MA_Push_Get(mt_dbl, nbf_ao*nbf_ao,'zora0scal', 82 & lzora0scal, izora0scal)) 83 & call errquit('zora_getv: zorasf',0, MA_ERR) 84c 85c == preliminaries == 86 do i= 1, nbf_ao*nbf_ao 87 dbl_mb(izora0+i-1)=0.d0 88 dbl_mb(izora0scal+i-1)=0.d0 89 enddo 90c 91c mirroring 92 dorepon=.true. 93 if (.not.rtdb_get(rtdb,'fock:mirrmat',mt_log,1,dorepon)) then 94 dorepon=nbf_ao.lt.3000 95 endif 96c 97 dorepl=.false. 98 if(ga_cluster_nnodes().gt.1.and.dorepon) then 99 docopy=.true. 100 dozero=.false. 101 dorepd=util_mirrmat(ipol,g_dens,g_densrep,docopy,dozero) 102 dorepl=dorepd 103 dorep_glob=0 104 if(dorepl) dorep_glob=1 105 call ga_igop(375,dorep_glob,1, '+') 106 dorepl=dorep_glob.eq.ga_nnodes() 107 if(dorepl) then 108 g_dens0(1)=g_dens(1) 109 g_dens(1)=g_densrep(1) 110 if(ipol.eq.2) then 111 g_dens0(2)=g_dens(2) 112 g_dens(2)=g_densrep(2) 113 endif 114 else 115 if(dorepd) then 116 call util_mirrstop(g_densrep(1)) 117 if(ipol.eq.2) call util_mirrstop(g_densrep(2)) 118 endif 119 if(ga_nodeid().eq.0) then 120 write(6,*) ' no mirroring in zora_getv_so' 121 call util_flush(6) 122 endif 123 endif 124 endif 125c 126c get zora model potential parameters of geometry 127 if (use_modelpotential) 128 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 129c 130c == generate the grid == 131 dummy(1) = 0.d0 132 dummy(2) = 0.d0 133 call grid_quadv0(rtdb, g_dens, g_zora_sf, nexc, rho_n, dummy, 134 & tmat) 135c 136c == initialization == 137 call ga_zero(g_zora_sf(1)) 138 if (ipol.gt.1) call ga_zero(g_zora_sf(2)) 139 call ga_zero(g_zora_scale_sf(1)) 140 if (ipol.gt.1) call ga_zero(g_zora_scale_sf(2)) 141c 142c == ao basis set info used by xc_eval_basis == 143 if(.not.int_normalize(rtdb,AO_bas_han)) 144 & call errquit('zora_getv_sf: int_normalize failed',0, BASIS_ERR) 145 if (.not.bas_numcont(AO_bas_han, ncontrset)) 146 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 147 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 148 & lbas_cent_info, ibas_cent_info)) 149 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 150 & MA_ERR) 151 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 152 & lbas_cset_info, ibas_cset_info)) 153 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 154 & MA_ERR) 155 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 156 & int_mb(ibas_cset_info), ncenters) 157c 158 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 159 & ldocset, idocset)) 160 & call errquit('zora_getv_sf: cannot allocate ccdocset', 161 . ncontrset, MA_ERR) 162 do i=1,ncontrset 163 log_mb(idocset+i-1)=.true. 164 enddo 165c 166 if (.not.MA_push_get(MT_int, ncenters, 'iniz', 167 & lniz, iniz)) call errquit("zora_getv_sf:iniz",0, MA_ERR) 168 do i= 1, ncenters 169 int_mb(iniz+i-1)=1 170 enddo 171c 172 nxyz = 3*ncenters 173 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 174 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 175 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 176 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 177 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 178 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 179 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 180 & Dbl_MB(ixyz), Dbl_MB(icharge))) 181 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 182 183 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 184 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 185 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 186 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 187 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 188 & 'quad pts buffer', lrqbuf, irqbuf)) 189 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 190c 191c == rewind grid file == 192 if (.not. grid_file_rewind()) 193 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 194 & UNKNOWN_ERR) 195c == loop over records in the grid file == 196 do iqsh = 1, n_rec_in_file 197c 198c == define the current range of radial shells and integration center == 199 call grid_file_read(n_per_rec, nqpts, ictr_buf, 200 & rad,dbl_mb(irqbuf),nsubb) 201 202 if(nqpts.gt.buffer_size) 203 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 204c 205c == loop over a subset of the grid == 206 istep=0 207 do ncube=1,nsubb 208c 209c put buf into currently used arrays qxyz and qwght 210 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 211 & dbl_mb(iqwght), nqpts, rad,istep) 212 213 if(nqpts.ne.0) then 214c 215c == compute the basis functions over the grid == 216 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 217 & l_rchi_atom,i_rchi_atom)) 218 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 219c 220 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 221 & l_rq,i_rq)) 222 & call errquit("zora_getv_sf:rq",0, MA_ERR) 223c 224c == delchi == 225 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 226 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 227 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 228c 229c == chi == 230 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 231 & 'chi_ao', lchi_ao, ichi_ao)) 232 & call errquit('zora_getv: chi_ao',0, MA_ERR) 233c 234 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 235 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 236 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 237 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 238 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 239 & int_mb(iniz), log_mb(idocset), 240 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 241c 242c calculate numerical overlap 243c call num_ovlp(dbl_mb(ichi_ao), dbl_mb(iqwght), 244c & nbf_ao, nqpts, dbl_mb(izora0)) 245c 246c == calculate spin-free zora contribution == 247 call calc_zora_sf(ao_bas_han, geom, ipol, g_dens, 248 & dbl_mb(ichi_ao), dbl_mb(idelchi_ao), dbl_mb(iqxyz), 249 & dbl_mb(iqwght), nbf_ao, nqpts, ncenters, 250 & dbl_mb(izora0), dbl_mb(izora0scal), 251 & ofinite,zetanuc_arr, 252 & Knucl, ! = .true. if including ONLY nuclear part in K ZORA 253 & use_modelpotential,gexpo,gcoef) 254c 255c == delete memory == 256 if(.not.MA_pop_stack(lchi_ao)) 257 & call errquit("zora_getv: pop chi_ao", 100, MA_ERR) 258 if(.not.MA_pop_stack(ldelchi_ao)) 259 & call errquit("zora_getv: pop delchi_ao", 100, MA_ERR) 260 if(.not.MA_pop_stack(l_rq)) 261 & call errquit("zora_getv: pop rq", 100, MA_ERR) 262 if(.not.MA_pop_stack(l_rchi_atom)) 263 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 264 endif ! nqpts 265 enddo ! ncube 266 end do ! iqsh 267c 268c == delete memory == 269 if(.not.MA_pop_stack(lrqbuf)) 270 & call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR) 271 if(.not.MA_pop_stack(lqwght)) 272 & call errquit("zora_getv_sf: pop qwght", 100, MA_ERR) 273 if(.not.MA_pop_stack(lqxyz)) 274 & call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR) 275 if(.not.MA_pop_stack(ltags)) 276 & call errquit("zora_getv_sf: pop tags", 100, MA_ERR) 277 if(.not.MA_pop_stack(lcharge)) 278 & call errquit("zora_getv_sf: pop charge", 100, MA_ERR) 279 if(.not.MA_pop_stack(lxyz)) 280 & call errquit("zora_getv_sf: pop xyz", 100, MA_ERR) 281 if(.not.MA_pop_stack(lniz)) 282 & call errquit("zora_getv_sf: pop niz", 100, MA_ERR) 283 if(.not.MA_pop_stack(ldocset)) 284 & call errquit("zora_getv_sf: pop docset", 100, MA_ERR) 285 if(.not.MA_pop_stack(lbas_cset_info)) 286 & call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR) 287 if(.not.MA_pop_stack(lbas_cent_info)) 288 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 289c 290 if(dorepl) then 291 call util_mirrstop(g_densrep(1)) 292 if(ipol.eq.2) call util_mirrstop(g_densrep(2)) 293 g_dens(1)=g_dens0(1) 294 if(ipol.eq.2) g_dens(2)=g_dens0(2) 295 endif 296c 297c == tally up over all the nodes == 298 call ga_dgop(msg_excrho, dbl_mb(izora0), nbf_ao*nbf_ao, '+') 299 call ga_dgop(msg_excrho, dbl_mb(izora0scal), nbf_ao*nbf_ao, '+') 300c 301c == scalar contribution == 302 call ga_zero(g_zora_sf(1)) 303 call ga_distribution(g_zora_sf(1), 304 & ga_nodeid(), ilo, ihi, jlo, jhi) 305 if (ilo.gt.0 .and. ilo.le.ihi) then 306 do j=jlo,jhi 307 call ga_put(g_zora_sf(1),ilo, ihi, jlo, jhi, 308 & dbl_mb(izora0+(jlo-1)*nbf_ao+ilo-1),nbf_ao) 309 enddo 310 endif 311 call ga_symmetrize(g_zora_sf(1)) 312c 313 call ga_zero(g_zora_scale_sf(1)) 314 call ga_distribution(g_zora_scale_sf(1), 315 & ga_nodeid(), ilo, ihi, jlo, jhi) 316 if (ilo.gt.0 .and. ilo.le.ihi) then 317 do j=jlo,jhi 318 call ga_put(g_zora_scale_sf(1),ilo, ihi, jlo, jhi, 319 & dbl_mb(izora0scal+(jlo-1)*nbf_ao+ilo-1),nbf_ao) 320 enddo 321 endif 322 call ga_symmetrize(g_zora_scale_sf(1)) 323c 324 if(ipol.eq.2) then 325 call ga_copy(g_zora_sf(1),g_zora_sf(2)) 326 call ga_copy(g_zora_scale_sf(1),g_zora_scale_sf(2)) 327 endif 328c 329c call ga_print(g_zora_sf(1)) 330c 331 call ga_sync() 332c 333 return 334 end 335c 336 subroutine num_ovlp(chi_ao,wght,nbf,npts,ovlp) 337 338 implicit none 339 340#include "stdio.fh" 341 342 integer nbf, npts 343 double precision chi_ao(npts,nbf),wght(npts) 344 double precision ovlp(nbf,nbf) 345 integer i, j, k 346 347 do i = 1, nbf 348 do j = 1, nbf 349 do k = 1, npts 350 ovlp(i,j)=ovlp(i,j)+chi_ao(k,i)*wght(k)*chi_ao(k,j) 351 enddo 352 enddo 353 enddo 354c 355 return 356 end 357