1 subroutine zora_getv_NMRCS_SR12(rtdb, 2 & g_dens, ! in : atomic density 3 & xyz_NMRcoords, ! in 4 & g_dia1, ! out 5 & g_dia2,g_dia3, ! out 6 & g_nmr1,g_nmr2, ! out: munu matrix 7 & nexc) 8c 9C$Id$ 10c Adapted from zora_getv_sf 11 12 implicit none 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_dia1,g_dia2(3),g_dia3(3,3) 31 integer g_nmr1(6),g_nmr2(18) 32 integer nexc 33c 34c == local variables == 35 integer i,j,k,t,a,m,n,ind,nij,ac 36 double precision rho_n 37 double precision tmat 38 double precision dummy(2) 39 integer iqsh, istep, nxyz, ncontrset 40 integer ixyz, lxyz, icharge, lcharge, itags, ltags 41 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 42 & ictr_buf,iqpts 43 double precision rad,ke 44 integer lbas_cset_info, ibas_cset_info, 45 & lbas_cent_info, ibas_cent_info, 46 & ldocset, idocset, 47 & l_rchi_atom,i_rchi_atom, 48 & l_rq,i_rq,lniz, iniz, 49 & lchi_ao, ichi_ao, 50 & ldelchi_ao, idelchi_ao 51 double precision xyz_NMRcoords(3) 52 integer ldia1,idia1, 53 & ldia2(3),idia2(3), 54 & ldia3(3,3),idia3(3,3) 55 integer lnmr1(3,3), inmr1(3,3), 56 & lnmr2(3,3,3),inmr2(3,3,3) 57 integer ind_tmn(2,3) 58 data ind_tmn / 2, 3, ! tmn=123 59 & 3, 1, ! tmn=231 60 & 1, 2 / ! tmn=312 61 integer inntsize,ddblsize,ok 62 logical grid_file_rewind 63 logical dft_mirrdens_start,dorepl 64 integer g_dens0(2),g_densrep(2),ii 65 integer nn 66 external dft_mirrdens_start 67 external grid_file_rewind 68c 69c model potential parameters 70 character*2 gelem(ncenters) 71 double precision gexpo(ncenters,50) 72 double precision gcoef(ncenters,50) 73c 74c == allocate memory == 75 do t=1,3 76 m=ind_tmn(1,t) 77 n=ind_tmn(2,t) 78 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 79 & 'lnmr',lnmr1(m,n),inmr1(m,n))) 80 & call errquit('zora_getv_NMR: nmr1',911,MA_ERR) 81 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 82 & 'lnmr',lnmr1(n,m),inmr1(n,m))) 83 & call errquit('zora_getv_NMR: nmr1',911,MA_ERR) 84 do a=1,3 85 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 86 & 'lnmr',lnmr2(a,m,n),inmr2(a,m,n))) 87 & call errquit('zora_getv_NMR: nmr2',911,MA_ERR) 88 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 89 & 'lnmr',lnmr2(a,n,m),inmr2(a,n,m))) 90 & call errquit('zora_getv_NMR: nmr2',911,MA_ERR) 91 enddo 92 enddo ! end-loop-t 93 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 94 & 'ldia1',ldia1,idia1)) 95 & call errquit('zora_getv_NMR: dia1',911,MA_ERR) 96 do i=1,3 97 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 98 & 'ldia1',ldia2(i),idia2(i))) 99 & call errquit('zora_getv_NMR: dia2',911,MA_ERR) 100 do j=1,3 101 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 102 & 'ldia1',ldia3(i,j),idia3(i,j))) 103 & call errquit('zora_getv_NMR: dia3',911,MA_ERR) 104 enddo 105 enddo ! end-loop-i 106c == preliminaries == 107 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(idia1),1) 108 do k=1,3 109 do j=1,3 110 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(idia3(k,j)),1) 111 enddo 112 enddo 113 do j=1,3 114 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(idia2(j)),1) 115 enddo 116 do t=1,3 117 m=ind_tmn(1,t) 118 n=ind_tmn(2,t) 119 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(inmr1(m,n)),1) 120 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(inmr1(n,m)),1) 121 do a=1,3 122 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(inmr2(a,m,n)),1) 123 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(inmr2(a,n,m)),1) 124 enddo ! end-loop-a 125 enddo ! end-loop-t 126c 127c get zora model potential parameters of geometry 128 if (use_modelpotential) 129 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 130c 131c == generate the grid == 132 dummy(1) = 0.d0 133 dummy(2) = 0.d0 134 dorepl = dft_mirrdens_start(g_dens,g_densrep,g_dens0, 135 i ipol) 136 call grid_quadv0(rtdb,g_dens,g_dia2,nexc,rho_n,dummy,tmat) 137c == ao basis set info used by xc_eval_basis == 138 if (.not.bas_numcont(AO_bas_han, ncontrset)) 139 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 140 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 141 & lbas_cent_info, ibas_cent_info)) 142 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 143 & MA_ERR) 144 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 145 & lbas_cset_info, ibas_cset_info)) 146 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 147 & MA_ERR) 148 149 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 150 & int_mb(ibas_cset_info), ncenters) 151 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 152 & ldocset, idocset)) 153 & call errquit('zora_getv_sf: cannot allocate ccdocset', 154 . ncontrset, MA_ERR) 155 do i=1,ncontrset 156 log_mb(idocset+i-1)=.true. 157 enddo 158 if(.not.MA_push_get(MT_int, ncenters, 'iniz', 159 & lniz, iniz)) 160 & call errquit("zora_getv_sf:iniz",0, MA_ERR) 161 do i= 1, ncenters 162 int_mb(iniz+i-1)=1 163 enddo 164 nxyz = 3*ncenters 165 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 166 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 167 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 168 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 169 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 170 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 171 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 172 & Dbl_MB(ixyz), Dbl_MB(icharge))) 173 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 174 175 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 176 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 177 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 178 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 179 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 180 & 'quad pts buffer', lrqbuf, irqbuf)) 181 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 182 183 if (.not. grid_file_rewind()) 184 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 185 & UNKNOWN_ERR) 186 187c == loop over records in the grid file == 188 do iqsh = 1, n_rec_in_file 189c == define the current range of radial shells and integration center == 190 call grid_file_read(n_per_rec, nqpts, ictr_buf, 191 & rad,dbl_mb(irqbuf),nsubb) 192 if(nqpts.gt.buffer_size) 193 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 194c == loop over a subset of the grid == 195 istep=0 196 do ncube=1,nsubb 197c put buf into currently used arrays qxyz and qwght 198 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 199 & dbl_mb(iqwght), nqpts, rad,istep) 200 201 if(nqpts.ne.0) then 202c == compute the basis functions over the grid == 203 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 204 & l_rchi_atom,i_rchi_atom)) 205 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 206c 207 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 208 & l_rq,i_rq)) 209 & call errquit("zora_getv_sf:rq",0, MA_ERR) 210c 211c == delchi == 212 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 213 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 214 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 215c 216c == chi == 217 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 218 & 'chi_ao', lchi_ao, ichi_ao)) 219 & call errquit('zora_getv: chi_ao',0, MA_ERR) 220 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 221 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 222 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 223 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 224 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 225 & int_mb(iniz), log_mb(idocset), 226 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 227 call calc_zora_NMRCS_SR(ao_bas_han,geom,ipol,g_dens, 228 & dbl_mb(ichi_ao), 229 & dbl_mb(idelchi_ao), 230 & dbl_mb(iqxyz),dbl_mb(iqwght), 231 & nbf_ao,nqpts,ncenters, 232 & xyz_NMRcoords, 233 & use_modelpotential, 234 & gexpo, gcoef, 235 & dbl_mb(idia1), ! out 236 & dbl_mb(idia2(1)), ! out 237 & dbl_mb(idia2(2)), ! out 238 & dbl_mb(idia2(3)), ! out 239 & dbl_mb(idia3(1,1)), ! out 240 & dbl_mb(idia3(1,2)), ! out 241 & dbl_mb(idia3(1,3)), ! out 242 & dbl_mb(idia3(2,1)), ! out 243 & dbl_mb(idia3(2,2)), ! out 244 & dbl_mb(idia3(2,3)), ! out 245 & dbl_mb(idia3(3,1)), ! out 246 & dbl_mb(idia3(3,2)), ! out 247 & dbl_mb(idia3(3,3)), ! out 248 & dbl_mb(inmr1(1,2)), ! out 249 & dbl_mb(inmr1(1,3)), ! out 250 & dbl_mb(inmr1(2,1)), ! out 251 & dbl_mb(inmr1(2,3)), ! out 252 & dbl_mb(inmr1(3,1)), ! out 253 & dbl_mb(inmr1(3,2)), ! out 254 & dbl_mb(inmr2(1,1,2)),! out 255 & dbl_mb(inmr2(1,1,3)),! out 256 & dbl_mb(inmr2(1,2,1)),! out 257 & dbl_mb(inmr2(1,2,3)),! out 258 & dbl_mb(inmr2(1,3,1)),! out 259 & dbl_mb(inmr2(1,3,2)),! out 260 & dbl_mb(inmr2(2,1,2)),! out 261 & dbl_mb(inmr2(2,1,3)),! out 262 & dbl_mb(inmr2(2,2,1)),! out 263 & dbl_mb(inmr2(2,2,3)),! out 264 & dbl_mb(inmr2(2,3,1)),! out 265 & dbl_mb(inmr2(2,3,2)),! out 266 & dbl_mb(inmr2(3,1,2)),! out 267 & dbl_mb(inmr2(3,1,3)),! out 268 & dbl_mb(inmr2(3,2,1)),! out 269 & dbl_mb(inmr2(3,2,3)),! out 270 & dbl_mb(inmr2(3,3,1)),! out 271 & dbl_mb(inmr2(3,3,2)))! out 272c == delete memory == 273 if(.not.MA_pop_stack(lchi_ao)) 274 & call errquit("zora_getv: pop chi_ao", 100, MA_ERR) 275 if(.not.MA_pop_stack(ldelchi_ao)) 276 & call errquit("zora_getv: pop delchi_ao", 100, MA_ERR) 277 if(.not.MA_pop_stack(l_rq)) 278 & call errquit("zora_getv: pop rq", 100, MA_ERR) 279 if(.not.MA_pop_stack(l_rchi_atom)) 280 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 281 endif ! nqpts 282 enddo ! ncube 283 end do ! iqsh 284c == delete memory == 285 if(.not.MA_pop_stack(lrqbuf)) 286 & call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR) 287 if(.not.MA_pop_stack(lqwght)) 288 & call errquit("zora_getv_sf: pop qwght", 100, MA_ERR) 289 if(.not.MA_pop_stack(lqxyz)) 290 & call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR) 291 if(.not.MA_pop_stack(ltags)) 292 & call errquit("zora_getv_sf: pop tags", 100, MA_ERR) 293 if(.not.MA_pop_stack(lcharge)) 294 & call errquit("zora_getv_sf: pop charge", 100, MA_ERR) 295 if(.not.MA_pop_stack(lxyz)) 296 & call errquit("zora_getv_sf: pop xyz", 100, MA_ERR) 297 if(.not.MA_pop_stack(lniz)) 298 & call errquit("zora_getv_sf: pop niz", 100, MA_ERR) 299 if(.not.MA_pop_stack(ldocset)) 300 & call errquit("zora_getv_sf: pop docset", 100, MA_ERR) 301 if(.not.MA_pop_stack(lbas_cset_info)) 302 & call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR) 303 if(.not.MA_pop_stack(lbas_cent_info)) 304 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 305c 306c == tally up over all the nodes == 307 call ga_sync() 308 if(dorepl) then 309 do ii=1,ipol 310 call util_mirrstop(g_densrep(ii)) 311 g_dens(ii)=g_dens0(ii) 312 enddo 313 endif 314 nn=nbf_ao*nbf_ao 315 call ga_dgop(msg_excrho,dbl_mb(idia1),nn, '+') 316 do k=1,3 317 call ga_dgop(msg_excrho,dbl_mb(idia2(k)),nn, '+') 318 do t=1,3 319 call ga_dgop(msg_excrho,dbl_mb(idia3(k,t)),nn,'+') 320 enddo ! end-loop-t 321 enddo ! end-loop-k 322 do t=1,3 323 m=ind_tmn(1,t) 324 n=ind_tmn(2,t) 325 call ga_dgop(msg_excrho,dbl_mb(inmr1(m,n)),nn, '+') 326 call ga_dgop(msg_excrho,dbl_mb(inmr1(n,m)),nn, '+') 327 do a=1,3 328 call ga_dgop(msg_excrho,dbl_mb(inmr2(a,m,n)),nn, '+') 329 call ga_dgop(msg_excrho,dbl_mb(inmr2(a,n,m)),nn, '+') 330 enddo ! end-loop-a 331 enddo ! end-loop-t 332c 333c == pack into a ga all NMR AOs == 334 call ga_zero(g_dia1) 335 call ga_put(g_dia1, 336 & 1,nbf_ao,1,nbf_ao,dbl_mb(idia1),nbf_ao) 337 call ga_symmetrize(g_dia1) 338 do i=1,3 339 call ga_zero(g_dia2(i)) 340 call ga_put(g_dia2(i), 341 & 1,nbf_ao,1,nbf_ao,dbl_mb(idia2(i)),nbf_ao) 342 call ga_symmetrize(g_dia2(i)) 343 do j=1,3 344 call ga_zero(g_dia3(i,j)) 345 call ga_put(g_dia3(i,j), 346 & 1,nbf_ao,1,nbf_ao,dbl_mb(idia3(i,j)),nbf_ao) 347 call ga_symmetrize(g_dia3(i,j)) 348 enddo ! end-loop-j 349 enddo ! end-loop-i 350 ac=1 351 do t=1,3 352 m=ind_tmn(1,t) 353 n=ind_tmn(2,t) 354 call ga_zero(g_nmr1(ac)) 355 if(ga_nodeid().eq.0) 356 G call ga_put(g_nmr1(ac), 357 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr1(m,n)),nbf_ao) 358 ac=ac+1 359 enddo ! end-loop-t 360 do t=1,3 361 m=ind_tmn(2,t) 362 n=ind_tmn(1,t) 363 call ga_zero(g_nmr1(ac)) 364 if(ga_nodeid().eq.0) 365 G call ga_put(g_nmr1(ac), 366 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr1(m,n)),nbf_ao) 367 ac=ac+1 368 enddo ! end-loop-t 369 ac=1 370 do a=1,3 371 do t=1,3 372 m=ind_tmn(1,t) 373 n=ind_tmn(2,t) 374 call ga_zero(g_nmr2(ac)) 375 if(ga_nodeid().eq.0) 376 G call ga_put(g_nmr2(ac), 377 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr2(a,m,n)),nbf_ao) 378 ac=ac+1 379 enddo ! end-loop-t 380 do t=1,3 381 m=ind_tmn(2,t) 382 n=ind_tmn(1,t) 383 call ga_zero(g_nmr2(ac)) 384 if(ga_nodeid().eq.0) 385 G call ga_put(g_nmr2(ac), 386 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr2(a,m,n)),nbf_ao) 387 ac=ac+1 388 enddo ! end-loop-t 389 enddo ! end-loop-a 390 call ga_sync() 391c ----- free memory ----------- START 392 if (.not.ma_free_heap(ldia1)) call errquit 393 & ('zora_getv_NMR: ma_free_heap of ldia1 failed',911,MA_ERR) 394 do i=1,3 395 if (.not.ma_free_heap(ldia2(i))) call errquit 396 & ('zora_getv_NMR: ma_free_heap of ldia2 failed',911,MA_ERR) 397 do j=1,3 398 if (.not.ma_free_heap(ldia3(i,j))) call errquit 399 & ('zora_getv_NMR: ma_free_heap of ldia3 failed',911,MA_ERR) 400 enddo ! end-loop-j 401 enddo 402 do t=1,3 403 m=ind_tmn(1,t) 404 n=ind_tmn(2,t) 405 if (.not.ma_free_heap(lnmr1(m,n))) call errquit 406 & ('zora_getv_NMR: ma_free_heap of lnmr1 failed',911,MA_ERR) 407 if (.not.ma_free_heap(lnmr1(n,m))) call errquit 408 & ('zora_getv_NMR: ma_free_heap of lnmr1 failed',911,MA_ERR) 409 do a=1,3 410 if (.not.ma_free_heap(lnmr2(a,m,n))) call errquit 411 & ('zora_getv_NMR: ma_free_heap of lnmr2 failed',911,MA_ERR) 412 if (.not.ma_free_heap(lnmr2(a,n,m))) call errquit 413 & ('zora_getv_NMR: ma_free_heap of lnmr2 failed',911,MA_ERR) 414 enddo 415 enddo 416c ----- free memory ----------- END 417 return 418 end 419 420 subroutine zora_getv_NMRCS_SR34(rtdb, 421 & g_dens, ! in : atomic density 422 & g_nmr, ! out 423 & g_nmr3,g_nmr4, ! out: munu matrix 424 & nexc) 425c 426C$Id$ 427c Adapted from zora_getv_sf 428 429 implicit none 430#include "rtdb.fh" 431#include "bas.fh" 432#include "cdft.fh" 433#include "errquit.fh" 434#include "mafdecls.fh" 435#include "global.fh" 436#include "geom.fh" 437#include "msgtypesf.h" 438#include "msgids.fh" 439#include "stdio.fh" 440#include "cgridfile.fh" 441#include "grid_cube.fh" 442#include "modelpotential.fh" 443c 444c == arguments == 445 integer rtdb 446 integer g_dens(2) 447 integer g_nmr(3), 448 & g_nmr3(3),g_nmr4(6) 449 integer nexc 450c == local variables == 451 integer i,j,k,t,a,m,n,ind,nij,ac 452 double precision rho_n 453 double precision tmat 454 double precision dummy(2) 455 integer iqsh, istep, nxyz, ncontrset 456 integer ixyz, lxyz, icharge, lcharge, itags, ltags 457 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 458 & ictr_buf,iqpts 459 double precision rad,ke 460 integer lbas_cset_info, ibas_cset_info, 461 & lbas_cent_info, ibas_cent_info, 462 & ldocset, idocset, 463 & l_rchi_atom,i_rchi_atom, 464 & l_rq,i_rq,lniz, iniz, 465 & lchi_ao, ichi_ao, 466 & ldelchi_ao, idelchi_ao 467 double precision xyz_NMRcoords(3) 468 integer lnmr(3) ,inmr(3), 469 & lnmr3(3) ,inmr3(3), 470 & lnmr4(3,3),inmr4(3,3) 471 integer ind_tmn(2,3) 472 data ind_tmn / 2, 3, ! tmn=123 473 & 3, 1, ! tmn=231 474 & 1, 2 / ! tmn=312 475 integer inntsize,ddblsize,ok 476 logical grid_file_rewind 477 logical dft_mirrdens_start,dorepl 478 integer g_dens0(2),g_densrep(2),ii 479 integer nn 480 external grid_file_rewind 481c 482c model potential parameters 483 character*2 gelem(ncenters) 484 double precision gexpo(ncenters,50) 485 double precision gcoef(ncenters,50) 486c 487c == allocate memory == 488 do t=1,3 489 m=ind_tmn(1,t) 490 n=ind_tmn(2,t) 491 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 492 & 'lnmr',lnmr4(m,n),inmr4(m,n))) 493 & call errquit('zora_getv_NMR: nmr4',911,MA_ERR) 494 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 495 & 'lnmr',lnmr4(n,m),inmr4(n,m))) 496 & call errquit('zora_getv_NMR: nmr4',911,MA_ERR) 497 enddo ! end-loop-t 498 do i=1,3 499 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 500 & 'lnmr',lnmr(i),inmr(i))) 501 & call errquit('zora_getv_NMR: nmr',911,MA_ERR) 502 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 503 & 'lnmr',lnmr3(i),inmr3(i))) 504 & call errquit('zora_getv_NMR: nmr3',911,MA_ERR) 505 enddo ! end-loop-i 506c == preliminaries == 507 do j=1,3 508 do i= 1, nbf_ao*nbf_ao 509 dbl_mb(inmr(j)+i-1) =0.d0 510 dbl_mb(inmr3(j)+i-1)=0.d0 511 enddo 512 enddo 513 do t=1,3 514 m=ind_tmn(1,t) 515 n=ind_tmn(2,t) 516 do i= 1, nbf_ao*nbf_ao 517 dbl_mb(inmr4(m,n)+i-1)=0.d0 518 dbl_mb(inmr4(n,m)+i-1)=0.d0 519 enddo ! end-loop-i 520 enddo ! end-loop-t 521c 522c get zora model potential parameters of geometry 523 if (use_modelpotential) 524 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 525c 526c == generate the grid == 527 dummy(1) = 0.d0 528 dummy(2) = 0.d0 529 dorepl = dft_mirrdens_start(g_dens,g_densrep,g_dens0, 530 i ipol) 531 call grid_quadv0(rtdb,g_dens,g_nmr(1),nexc,rho_n,dummy,tmat) 532c == ao basis set info used by xc_eval_basis == 533 if (.not.bas_numcont(AO_bas_han, ncontrset)) 534 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 535 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 536 & lbas_cent_info, ibas_cent_info)) 537 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 538 & MA_ERR) 539 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 540 & lbas_cset_info, ibas_cset_info)) 541 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 542 & MA_ERR) 543 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 544 & int_mb(ibas_cset_info), ncenters) 545 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 546 & ldocset, idocset)) 547 & call errquit('zora_getv_sf: cannot allocate ccdocset', 548 . ncontrset, MA_ERR) 549 do i=1,ncontrset 550 log_mb(idocset+i-1)=.true. 551 enddo 552 if(.not.MA_push_get(MT_int, ncenters, 'iniz', 553 & lniz, iniz)) 554 & call errquit("zora_getv_sf:iniz",0, MA_ERR) 555 do i= 1, ncenters 556 int_mb(iniz+i-1)=1 557 enddo 558 nxyz = 3*ncenters 559 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 560 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 561 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 562 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 563 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 564 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 565 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 566 & Dbl_MB(ixyz), Dbl_MB(icharge))) 567 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 568 569 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 570 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 571 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 572 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 573 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 574 & 'quad pts buffer', lrqbuf, irqbuf)) 575 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 576 577 if (.not. grid_file_rewind()) 578 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 579 & UNKNOWN_ERR) 580c 581c == loop over records in the grid file == 582 do iqsh = 1, n_rec_in_file 583c == define the current range of radial shells and integration center == 584 call grid_file_read(n_per_rec, nqpts, ictr_buf, 585 & rad,dbl_mb(irqbuf),nsubb) 586 if(nqpts.gt.buffer_size) 587 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 588c == loop over a subset of the grid == 589 istep=0 590 do ncube=1,nsubb 591c put buf into currently used arrays qxyz and qwght 592 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 593 & dbl_mb(iqwght), nqpts, rad,istep) 594 if(nqpts.ne.0) then 595c == compute the basis functions over the grid == 596 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 597 & l_rchi_atom,i_rchi_atom)) 598 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 599 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 600 & l_rq,i_rq)) 601 & call errquit("zora_getv_sf:rq",0, MA_ERR) 602c == delchi == 603 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 604 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 605 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 606c == chi == 607 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 608 & 'chi_ao', lchi_ao, ichi_ao)) 609 & call errquit('zora_getv: chi_ao',0, MA_ERR) 610 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 611 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 612 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 613 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 614 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 615 & int_mb(iniz), log_mb(idocset), 616 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 617 call calc_NMRCS_SR_F1ij(ao_bas_han,geom,ipol,g_dens, 618 & dbl_mb(ichi_ao), 619 & dbl_mb(idelchi_ao), 620 & dbl_mb(iqxyz),dbl_mb(iqwght), 621 & nbf_ao,nqpts,ncenters, 622 & use_modelpotential, 623 & gexpo, gcoef, 624 & dbl_mb(inmr(1)), ! out 625 & dbl_mb(inmr(2)), ! out 626 & dbl_mb(inmr(3)), ! out 627 & dbl_mb(inmr3(1)), ! out 628 & dbl_mb(inmr3(2)), ! out 629 & dbl_mb(inmr3(3)), ! out 630 & dbl_mb(inmr4(1,2)), ! out 631 & dbl_mb(inmr4(1,3)), ! out 632 & dbl_mb(inmr4(2,1)), ! out 633 & dbl_mb(inmr4(2,3)), ! out 634 & dbl_mb(inmr4(3,1)), ! out 635 & dbl_mb(inmr4(3,2))) ! out 636c == delete memory == 637 if(.not.MA_pop_stack(lchi_ao)) 638 & call errquit("zora_getv: pop chi_ao", 100, MA_ERR) 639 if(.not.MA_pop_stack(ldelchi_ao)) 640 & call errquit("zora_getv: pop delchi_ao", 100, MA_ERR) 641 if(.not.MA_pop_stack(l_rq)) 642 & call errquit("zora_getv: pop rq", 100, MA_ERR) 643 if(.not.MA_pop_stack(l_rchi_atom)) 644 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 645 endif ! nqpts 646 enddo ! ncube 647 end do ! iqsh 648c == delete memory == 649 if(.not.MA_pop_stack(lrqbuf)) 650 & call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR) 651 if(.not.MA_pop_stack(lqwght)) 652 & call errquit("zora_getv_sf: pop qwght", 100, MA_ERR) 653 if(.not.MA_pop_stack(lqxyz)) 654 & call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR) 655 if(.not.MA_pop_stack(ltags)) 656 & call errquit("zora_getv_sf: pop tags", 100, MA_ERR) 657 if(.not.MA_pop_stack(lcharge)) 658 & call errquit("zora_getv_sf: pop charge", 100, MA_ERR) 659 if(.not.MA_pop_stack(lxyz)) 660 & call errquit("zora_getv_sf: pop xyz", 100, MA_ERR) 661 if(.not.MA_pop_stack(lniz)) 662 & call errquit("zora_getv_sf: pop niz", 100, MA_ERR) 663 if(.not.MA_pop_stack(ldocset)) 664 & call errquit("zora_getv_sf: pop docset", 100, MA_ERR) 665 if(.not.MA_pop_stack(lbas_cset_info)) 666 & call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR) 667 if(.not.MA_pop_stack(lbas_cent_info)) 668 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 669c 670c == tally up over all the nodes == 671 call ga_sync() 672 if(dorepl) then 673 do ii=1,ipol 674 call util_mirrstop(g_densrep(ii)) 675 g_dens(ii)=g_dens0(ii) 676 enddo 677 endif 678 nn=nbf_ao*nbf_ao 679 do k=1,3 680 call ga_dgop(msg_excrho,dbl_mb(inmr(k)), nn, '+') 681 call ga_dgop(msg_excrho,dbl_mb(inmr3(k)), nn, '+') 682 enddo ! end-loop-k 683 do t=1,3 684 m=ind_tmn(1,t) 685 n=ind_tmn(2,t) 686 call ga_dgop(msg_excrho,dbl_mb(inmr4(m,n)), nn, '+') 687 call ga_dgop(msg_excrho,dbl_mb(inmr4(n,m)), nn, '+') 688 enddo ! end-loop-t 689c 690c == pack into a ga all NMR AOs == 691 do i=1,3 692 call ga_zero(g_nmr(i)) 693 if(ga_nodeid().eq.0) 694 G call ga_put(g_nmr(i), 695 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr(i)),nbf_ao) 696 call ga_symmetrize(g_nmr(i)) 697 call ga_zero(g_nmr3(i)) 698 if(ga_nodeid().eq.0) 699 G call ga_put(g_nmr3(i), 700 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr3(i)),nbf_ao) 701 enddo ! end-loop-i 702 ac=1 703 do t=1,3 704 m=ind_tmn(1,t) 705 n=ind_tmn(2,t) 706 call ga_zero(g_nmr4(ac)) 707 if(ga_nodeid().eq.0) 708 G call ga_put(g_nmr4(ac), 709 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr4(m,n)),nbf_ao) 710 ac=ac+1 711 enddo ! end-loop-t 712 do t=1,3 713 m=ind_tmn(2,t) 714 n=ind_tmn(1,t) 715 call ga_zero(g_nmr4(ac)) 716 if(ga_nodeid().eq.0) 717 G call ga_put(g_nmr4(ac), 718 & 1,nbf_ao,1,nbf_ao,dbl_mb(inmr4(m,n)),nbf_ao) 719 ac=ac+1 720 enddo ! end-loop-t 721 call ga_sync() 722 723c ----- free memory ----------- START 724 do i=1,3 725 if (.not.ma_free_heap(lnmr(i))) call errquit 726 & ('zora_getv_NMR: ma_free_heap of lnmr failed',911,MA_ERR) 727 if (.not.ma_free_heap(lnmr3(i))) call errquit 728 & ('zora_getv_NMR: ma_free_heap of lnmr3 failed',911,MA_ERR) 729 enddo 730 do t=1,3 731 m=ind_tmn(1,t) 732 n=ind_tmn(2,t) 733 if (.not.ma_free_heap(lnmr4(m,n))) call errquit 734 & ('zora_getv_NMR: ma_free_heap of lnmr4 failed',911,MA_ERR) 735 if (.not.ma_free_heap(lnmr4(n,m))) call errquit 736 & ('zora_getv_NMR: ma_free_heap of lnmr4 failed',911,MA_ERR) 737 enddo ! end-loop-t 738c ----- free memory ----------- END 739 return 740 end 741 742