1 subroutine zora_getv_HFine_slow( 2 & rtdb, 3 & g_dens, ! in : atomic density 4 & ofinite, ! in : = .true. if Gaussian Nucl. Model requested 5 & atmass, ! in : atomic mass 6 & xyz_NMRcoords, ! in : nuclear coordinates 7 & g_zpso, ! out : ZPSO term 8 & g_fcsd, ! out : FC+SD (v,u) term 9 & nexc) 10c 11C$Id$ 12c Adapted from zora_getv_so 13 14 implicit none 15#include "rtdb.fh" 16#include "bas.fh" 17#include "cdft.fh" 18#include "errquit.fh" 19#include "mafdecls.fh" 20#include "global.fh" 21#include "geom.fh" 22#include "msgtypesf.h" 23#include "msgids.fh" 24#include "stdio.fh" 25#include "cgridfile.fh" 26#include "grid_cube.fh" 27#include "modelpotential.fh" 28c 29c == arguments == 30 integer rtdb 31 integer g_dens(2) 32 integer g_zpso(3),g_fcsd(3,3) 33 integer nexc 34c 35c == local variables == 36 integer i,j,k,t,n,ind,nij 37 double precision rho_n 38 double precision tmat 39 double precision dummy(2) 40 integer iqsh, istep, nxyz, ncontrset 41 integer ixyz, lxyz, icharge, lcharge, itags, ltags 42 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 43 & ictr_buf,iqpts 44 double precision rad,ke 45 integer lbas_cset_info, ibas_cset_info, 46 & lbas_cent_info, ibas_cent_info, 47 & ldocset, idocset, 48 & l_rchi_atom,i_rchi_atom, 49 & l_rq,i_rq,l_iniz, k_iniz, 50 & lchi_ao, ichi_ao, 51 & ldelchi_ao, idelchi_ao 52 integer lzpso(3),izpso(3), 53 & lfcsd(3,3),ifcsd(3,3) 54 integer inntsize,ddblsize,ok 55 double precision xyz_NMRcoords(3),atmass 56 double precision chi_cntr(3,nbf_ao) 57 58 logical grid_file_rewind,ofinite 59 external grid_file_rewind,calc_zora_HFine_slow, 60 & ga_antisymmetrize 61c 62c model potential parameters 63 character*2 gelem(ncenters) 64 double precision gexpo(ncenters,50) 65 double precision gcoef(ncenters,50) 66c 67c == allocate memory == 68 do i=1,3 69 if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao, 70 & 'lnmr',lzpso(i),izpso(i))) 71 & call errquit('zora_getv_HFine: zpso',911,MA_ERR) 72 enddo 73 do i=1,3 74 do j=i,3 75 if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao, 76 & 'lfcsd',lfcsd(i,j),ifcsd(i,j))) 77 & call errquit('zora_getv_HFine: fcsd',911,MA_ERR) 78 enddo 79 enddo ! end-loop-i 80c == preliminaries == 81 call dfill(3*nbf_ao*nbf_ao,0d0,dbl_mb(izpso(j)),1) 82 do k=1,3 83 do j=k,3 84 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(ifcsd(k,j)),1) 85 enddo 86 enddo 87c 88c get zora model potential parameters of geometry 89 if (use_modelpotential) 90 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 91c 92c == generate the grid == 93 dummy(1) = 0.d0 94 dummy(2) = 0.d0 95 call grid_quadv0(rtdb,g_dens,g_zpso,nexc,rho_n,dummy,tmat) 96c == ao basis set info used by xc_eval_basis == 97 if (.not.bas_numcont(AO_bas_han, ncontrset)) 98 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 99 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 100 & lbas_cent_info, ibas_cent_info)) 101 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 102 & MA_ERR) 103 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 104 & lbas_cset_info, ibas_cset_info)) 105 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 106 & MA_ERR) 107 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 108 & int_mb(ibas_cset_info), ncenters) 109 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 110 & ldocset, idocset)) 111 & call errquit('zora_getv_sf: cannot allocate ccdocset', 112 . ncontrset, MA_ERR) 113 do i=1,ncontrset 114 log_mb(idocset+i-1)=.true. 115 enddo 116 if(.not.MA_push_get(MT_int, ncenters, 'iniz', 117 & l_iniz, k_iniz)) 118 & call errquit("zora_getv_sf:iniz",0, MA_ERR) 119 do i= 1, ncenters 120 int_mb(k_iniz+i-1)=1 121 enddo 122 nxyz = 3*ncenters 123 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 124 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 125 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 126 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 127 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 128 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 129 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 130 & Dbl_MB(ixyz), Dbl_MB(icharge))) 131 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 132 133 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 134 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 135 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 136 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 137 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 138 & 'quad pts buffer', lrqbuf, irqbuf)) 139 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 140 141 if (.not. grid_file_rewind()) 142 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 143 & UNKNOWN_ERR) 144c 145c == loop over records in the grid file == 146 do iqsh = 1, n_rec_in_file 147c 148c == define the current range of radial shells and integration center == 149 call grid_file_read(n_per_rec, nqpts, ictr_buf, 150 & rad,dbl_mb(irqbuf),nsubb) 151 152 if(nqpts.gt.buffer_size) 153 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 154c 155c == loop over a subset of the grid == 156 istep=0 157 do ncube=1,nsubb 158c 159c put buf into currently used arrays qxyz and qwght 160 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 161 & dbl_mb(iqwght), nqpts, rad,istep) 162 163 if(nqpts.ne.0) then 164c 165c == compute the basis functions over the grid == 166 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 167 & l_rchi_atom,i_rchi_atom)) 168 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 169c 170 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 171 & l_rq,i_rq)) 172 & call errquit("zora_getv_sf:rq",0, MA_ERR) 173c 174c == delchi == 175 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 176 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 177 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 178c 179c == chi == 180 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 181 & 'chi_ao', lchi_ao, ichi_ao)) 182 & call errquit('zora_getv: chi_ao',0, MA_ERR) 183 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 184 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 185 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 186 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 187 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 188 & int_mb(k_iniz), log_mb(idocset), 189 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 190 call calc_zora_HFine_slow( 191 & ao_bas_han,geom,ipol,g_dens, 192 & dbl_mb(ichi_ao), 193 & dbl_mb(idelchi_ao), 194 & dbl_mb(iqxyz),dbl_mb(iqwght), 195 & nbf_ao,nqpts,ncenters, 196 & ofinite, 197 & atmass, 198 & xyz_NMRcoords, 199 & use_modelpotential, 200 & gexpo, 201 & gcoef, 202 & dbl_mb(izpso(1)), ! out 203 & dbl_mb(izpso(2)), ! out 204 & dbl_mb(izpso(3)), ! out 205 & dbl_mb(ifcsd(1,1)),! out 206 & dbl_mb(ifcsd(1,2)),! out 207 & dbl_mb(ifcsd(1,3)),! out 208c & dbl_mb(ifcsd(2,1)),! out 209 & dbl_mb(ifcsd(2,2)),! out 210 & dbl_mb(ifcsd(2,3)),! out 211c & dbl_mb(ifcsd(3,1)),! out 212c & dbl_mb(ifcsd(3,2)),! out 213 & dbl_mb(ifcsd(3,3)))! out 214 215c == delete memory == 216 if(.not.MA_chop_stack(l_rchi_atom)) 217 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 218 endif ! nqpts 219 enddo ! ncube 220 end do ! iqsh 221c 222c == delete memory == 223 if(.not.MA_chop_stack(lbas_cent_info)) 224 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 225c 226c == tally up over all the nodes == 227 do k=1,3 228 call ga_dgop(msg_excrho+k, 229 D dbl_mb(izpso(k)), nbf_ao*nbf_ao, '+') 230 do t=k,3 231c do t=1,3 232 call ga_dgop(msg_excrho+t*k, 233 D dbl_mb(ifcsd(k,t)),nbf_ao*nbf_ao, '+') 234 enddo ! end-loop-t 235 enddo ! end-loop-k 236c 237c == pack into a ga all HFine AOs == 238 if(ga_nodeid().eq.0) then 239 do i=1,3 240 call ga_zero(g_zpso(i)) 241 call ga_put(g_zpso(i), 242 & 1,nbf_ao,1,nbf_ao,dbl_mb(izpso(i)),nbf_ao) 243 do j=i,3 244 call ga_zero(g_fcsd(i,j)) 245 call ga_put(g_fcsd(i,j), 246 & 1,nbf_ao,1,nbf_ao,dbl_mb(ifcsd(i,j)),nbf_ao) 247 if(i.ne.j) then 248 call ga_copy(g_fcsd(i,j),g_fcsd(j,i)) 249 endif 250 enddo ! end-loop-j 251 enddo ! end-loop-i 252 endif 253 call ga_sync() 254c ---- Free memory 255 if (.not.ma_chop_stack(lzpso(1))) call errquit 256 & ('zora_getv_HFine: ma_chop_stack of lzpso failed', 257 & 911,MA_ERR) 258 return 259 end 260c ++++++++++++++++++++++++++++++++ 261c ++++++++++++++++++++++++++++++++ 262 subroutine zora_getv_HFine_fast( 263 & rtdb, 264 & g_dens, ! in : atomic density 265 & ofinite, ! in : = .true. if Gaussian Nucl. Model requested 266 & zetanuc_arr, ! in : sqrt(zetanuc(i)) i=1,natoms 267 & zetanuc_slc, ! in : zetanuc(i) 268 & Knucl, 269 & xyz_NMRcoords, ! in : nuclear coordinates 270 & g_zpso, ! out : ZPSO term 271 & g_fcsd, ! out : FC+SD (v,u) term 272 & nexc) 273c 274C$Id$ 275c Adapted from zora_getv_so 276 277 implicit none 278#include "rtdb.fh" 279#include "bas.fh" 280#include "cdft.fh" 281#include "errquit.fh" 282#include "mafdecls.fh" 283#include "global.fh" 284#include "geom.fh" 285#include "msgtypesf.h" 286#include "msgids.fh" 287#include "stdio.fh" 288#include "cgridfile.fh" 289#include "grid_cube.fh" 290#include "modelpotential.fh" 291c 292c == arguments == 293 integer rtdb 294 integer g_dens(2) 295 integer g_zpso(6),g_fcsd(3,3) 296 integer nexc 297c 298c == local variables == 299 integer i,j,k,t,m,n,ind,nij,ac 300 double precision rho_n 301 double precision tmat 302 double precision dummy(2) 303 integer iqsh, istep, nxyz, ncontrset 304 integer ixyz, lxyz, icharge, lcharge, itags, ltags 305 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 306 & ictr_buf,iqpts 307 double precision rad,ke 308 double precision zetanuc_arr(*) ! For Gaussian Nuclear Model 309 double precision zetanuc_slc ! For Gaussian Nuclear Model 310 integer lbas_cset_info, ibas_cset_info, 311 & lbas_cent_info, ibas_cent_info, 312 & ldocset, idocset, 313 & l_rchi_atom,i_rchi_atom, 314 & l_rq,i_rq, 315 & lchi_ao, ichi_ao, 316 & ldelchi_ao, idelchi_ao 317 integer lzpso(3,3),izpso(3,3), 318 & lfcsd(3,3),ifcsd(3,3) 319 integer inntsize,ddblsize,ok,atn 320 double precision xyz_NMRcoords(3),atmass 321 double precision chi_cntr(3,nbf_ao) 322 logical grid_file_rewind,ofinite,Knucl 323 integer ind_tmn(2,3) 324 data ind_tmn / 2, 3, ! tmn=123 325 & 3, 1, ! tmn=231 326 & 1, 2 / ! tmn=312 327 logical dft_mirrdens_start,dorepl 328 integer g_dens0(2),g_densrep(2),ii 329 external dft_mirrdens_start 330 external grid_file_rewind 331c 332c model potential parameters 333 character*2 gelem(ncenters) 334 double precision gexpo(ncenters,50) 335 double precision gcoef(ncenters,50) 336c mbf 337 integer grid_nbfm 338 external grid_nbfm 339 character*32 pname 340 double precision acc_ao_gauss 341 integer mbf_ao,k_expo,l_expo,l_ifin,k_ifin,l_iniz,k_iniz, 342 k k_ibf_ao,l_ibf_ao 343 pname = 'zora_getv_hfine: ' 344c 345c == allocate memory == 346 do t=1,3 347 m=ind_tmn(1,t) 348 n=ind_tmn(2,t) 349 if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao, 350 & 'lnmr',lzpso(m,n),izpso(m,n))) 351 & call errquit('zora_getv_NMR: zpso',911,MA_ERR) 352 if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao, 353 & 'lnmr',lzpso(n,m),izpso(n,m))) 354 & call errquit('zora_getv_NMR: zpso',911,MA_ERR) 355 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(izpso(m,n)),1) 356 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(izpso(n,m)),1) 357 enddo ! end-loop-t 358 do i=1,3 359 do j=i,3 360 if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao, 361 & 'lfcsd',lfcsd(i,j),ifcsd(i,j))) 362 & call errquit('zora_getv_HFine: fcsd',911,MA_ERR) 363 call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(ifcsd(i,j)),1) 364 enddo 365 enddo ! end-loop-i 366c 367c get zora model potential parameters of geometry 368 if (use_modelpotential) 369 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 370c 371c == generate the grid == 372 dummy(1) = 0.d0 373 dummy(2) = 0.d0 374 dorepl = dft_mirrdens_start(g_dens,g_densrep,g_dens0, 375 i ipol) 376 call grid_quadv0(rtdb,g_dens,g_zpso,nexc,rho_n,dummy,tmat) 377c == ao basis set info used by xc_eval_basis == 378 if (.not.bas_numcont(AO_bas_han, ncontrset)) 379 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 380 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 381 & lbas_cent_info, ibas_cent_info)) 382 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 383 & MA_ERR) 384 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 385 & lbas_cset_info, ibas_cset_info)) 386 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 387 & MA_ERR) 388 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 389 & int_mb(ibas_cset_info), ncenters) 390 nxyz = 3*ncenters 391 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 392 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 393 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 394 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 395 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 396 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 397 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 398 & Dbl_MB(ixyz), Dbl_MB(icharge))) 399 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 400 401 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 402 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 403 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 404 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 405 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 406 & 'quad pts buffer', lrqbuf, irqbuf)) 407 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 408 409 if (.not. grid_file_rewind()) 410 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 411 & UNKNOWN_ERR) 412 if (.not.MA_Push_Get(MT_Dbl,nbf_ao_mxprim,'expo',l_expo,k_expo)) 413 & call errquit(pname//'cannot allocate expo',0, MA_ERR) 414 if (.not.MA_Push_Get(mt_int, nbf_ao, 'ibf_ao', l_ibf_ao, 415 & k_ibf_ao)) 416 & call errquit(pname//'cannot allocate ibf_ao',2, 417 & MA_ERR) 418 if (.not.MA_Push_get(MT_int,ncenters,'atom list',l_iniz,k_iniz)) 419 & call errquit(pname//'cannot allocate iniz',0, MA_ERR) 420 if (.not.MA_Push_get(MT_int,ncenters,'atom list',l_ifin,k_ifin)) 421 & call errquit(pname//'cannot allocate fin',0, MA_ERR) 422 do i= 1, ncenters 423 int_mb(k_iniz+i-1)=1 424 enddo 425 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 426 & ldocset, idocset)) 427 & call errquit('zora_getv_sf: cannot allocate ccdocset', 428 . ncontrset, MA_ERR) 429 do i=1,ncontrset 430 log_mb(idocset+i-1)=.true. 431 enddo 432c 433c == loop over records in the grid file == 434 do iqsh = 1, n_rec_in_file 435c 436c == define the current range of radial shells and integration center == 437 call grid_file_read(n_per_rec, nqpts, ictr_buf, 438 & rad,dbl_mb(irqbuf),nsubb) 439 440 if(nqpts.gt.buffer_size) 441 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 442c 443c == loop over a subset of the grid == 444 istep=0 445 do ncube=1,nsubb 446c 447c put buf into currently used arrays qxyz and qwght 448 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 449 & dbl_mb(iqwght), nqpts, rad,istep) 450 451 if(nqpts.ne.0) then 452 mbf_ao=nbf_ao 453 454crestrict to mbf_ao subset 455c acc_ao_gauss = dble(max(iaoacc,25)) 456 acc_ao_gauss = dble(iaoacc) 457 call icopy(mbf_ao, 0,0, int_mb(k_ibf_ao), 1) 458 mbf_ao=grid_nbfm(ao_bas_han, ncenters, 459 & ictr_buf,rad, 460 Q dbl_mb(ixyz),dbl_mb(iqxyz),nqpts, 461 & int_mb(k_ibf_ao), log_mb(idocset), 462 I int_mb(k_iniz), int_mb(k_ifin), dbl_mb(k_expo), 463 & minexp,ldiff,acc_ao_gauss,iatype_pt_chg) 464 465c 466c == compute the basis functions over the grid == 467 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 468 & l_rchi_atom,i_rchi_atom)) 469 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 470c 471 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 472 & l_rq,i_rq)) 473 & call errquit("zora_getv_sf:rq",0, MA_ERR) 474c 475c == delchi == 476 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 477 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 478 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 479c 480c == chi == 481 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 482 & 'chi_ao', lchi_ao, ichi_ao)) 483 & call errquit('zora_getv: chi_ao',0, MA_ERR) 484 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 485 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 486 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 487 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 488 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 489 & int_mb(k_iniz), log_mb(idocset), 490 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 491 call calc_zora_HFine_fast( 492 & ao_bas_han,geom,ipol,g_dens, 493 & dbl_mb(ichi_ao), 494 & dbl_mb(idelchi_ao), 495 & dbl_mb(iqxyz),dbl_mb(iqwght), 496 & nbf_ao,mbf_ao,int_mb(k_ibf_ao), 497 N nqpts,ncenters, 498 & ofinite, 499 & zetanuc_arr, 500 & zetanuc_slc, 501 & Knucl, 502 & xyz_NMRcoords, 503 & use_modelpotential, 504 & gexpo, 505 & gcoef, 506 & dbl_mb(izpso(1,2)),! out 507 & dbl_mb(izpso(2,3)),! out 508 & dbl_mb(izpso(3,1)),! out 509 & dbl_mb(ifcsd(1,1)),! out 510 & dbl_mb(ifcsd(1,2)),! out 511 & dbl_mb(ifcsd(1,3)),! out 512c & dbl_mb(ifcsd(2,1)),! out 513 & dbl_mb(ifcsd(2,2)),! out 514 & dbl_mb(ifcsd(2,3)),! out 515c & dbl_mb(ifcsd(3,1)),! out 516c & dbl_mb(ifcsd(3,2)),! out 517 & dbl_mb(ifcsd(3,3)))! out 518c == delete memory == 519 if(.not.MA_chop_stack(l_rchi_atom)) 520 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 521 endif ! nqpts 522 enddo ! ncube 523 end do ! iqsh 524c 525c == delete memory == 526 if(.not.MA_chop_stack(lbas_cent_info)) 527 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 528 call ga_sync() 529 if(dorepl) then 530 do ii=1,ipol 531 call util_mirrstop(g_densrep(ii)) 532 g_dens(ii)=g_dens0(ii) 533 enddo 534 endif 535c 536c == tally up over all the nodes == 537 do k=1,3 538 m=ind_tmn(1,k) 539 n=ind_tmn(2,k) 540 if(k.eq.1) then 541 call ga_mask_sync(.true.,.false.) 542 else 543 call ga_mask_sync(.false.,.false.) 544 endif 545 call ga_dgop(msg_excrho,dbl_mb(izpso(m,n)),nbf_ao*nbf_ao,'+') 546 call dcopy( 547 D nbf_ao*nbf_ao,dbl_mb(izpso(n,m)),1,dbl_mb(izpso(m,n)),1) 548c call ga_dgop(msg_excrho,dbl_mb(izpso(n,m)),nbf_ao*nbf_ao,'+') 549 do t=k,3 550 if(t.eq.3) then 551 call ga_mask_sync(.false.,.true.) 552 else 553 call ga_mask_sync(.false.,.false.) 554 endif 555 call ga_dgop(msg_excrho,dbl_mb(ifcsd(k,t)),nbf_ao*nbf_ao,'+') 556c ifcsd upper triangle gone 557c if(t.ne.k) call dcopy( 558c D nbf_ao*nbf_ao,dbl_mb(ifcsd(k,t)),1,dbl_mb(ifcsd(t,k)),1) 559 enddo ! end-loop-t 560 enddo ! end-loop-k 561c 562c == pack into a ga all HFine AOs == 563 ac=1 564 do t=1,3 565 m=ind_tmn(1,t) 566 n=ind_tmn(2,t) 567 call ga_zero(g_zpso(ac)) 568 if(ga_nodeid().eq.0) 569 c call ga_put(g_zpso(ac), 570 & 1,nbf_ao,1,nbf_ao,dbl_mb(izpso(m,n)),nbf_ao) 571 ac=ac+1 572 enddo ! end-loop-t 573 do t=1,3 574 m=ind_tmn(2,t) 575 n=ind_tmn(1,t) 576 call ga_zero(g_zpso(ac)) 577 if(ga_nodeid().eq.0) 578 c call ga_put(g_zpso(ac), 579 & 1,nbf_ao,1,nbf_ao,dbl_mb(izpso(m,n)),nbf_ao) 580 ac=ac+1 581 enddo ! end-loop-t 582 call ga_sync() 583 do i=1,3 584 do j=i,3 585 call ga_zero(g_fcsd(i,j)) 586 if(ga_nodeid().eq.0) 587 c call ga_put(g_fcsd(i,j), 588 & 1,nbf_ao,1,nbf_ao,dbl_mb(ifcsd(i,j)),nbf_ao) 589 if(i.ne.j) call ga_copy(g_fcsd(i,j),g_fcsd(j,i)) 590 enddo ! end-loop-j 591 enddo ! end-loop-i 592 call ga_sync() 593c ----- free memory ----------- START 594 m=ind_tmn(1,1) 595 n=ind_tmn(2,1) 596 if (.not.ma_chop_stack(lzpso(m,n))) call errquit 597 & ('zora_getv_HFine: ma_chop_stack of lzpso failed', 598 & 911,MA_ERR) 599 600c ----- free memory ----------- END 601 return 602 end 603 604 subroutine zora_getv_NMRHFine_F1ji( 605 & rtdb, 606 & g_dens, ! in : atomic density 607 & g_hfine, ! out: 608 & natoms, ! in: nr. atoms 609 & ofinite, ! in: = .true. if Gaussian Nucl. Model of charges requested 610 & zetanuc_arr, ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model 611 & Knucl, ! in: = .true. if K_ZORA(V=Nuclear pot. only) 612 & nexc) 613c 614C$Id$ 615c Adapted from zora_getv_sf 616 617 implicit none 618#include "rtdb.fh" 619#include "bas.fh" 620#include "cdft.fh" 621#include "errquit.fh" 622#include "mafdecls.fh" 623#include "global.fh" 624#include "geom.fh" 625#include "msgtypesf.h" 626#include "msgids.fh" 627#include "stdio.fh" 628#include "cgridfile.fh" 629#include "grid_cube.fh" 630#include "modelpotential.fh" 631c 632c == arguments == 633 integer rtdb 634 integer g_dens(2) 635 integer g_hfine(3) 636 integer nexc 637c == local variables == 638 integer i,j,k,t,a,m,n,ind,nij,ac 639 double precision rho_n 640 double precision tmat 641 double precision dummy(2) 642 integer iqsh, istep, nxyz, ncontrset 643 integer ixyz, lxyz, icharge, lcharge, itags, ltags 644 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 645 & ictr_buf,iqpts 646 double precision rad,ke 647 integer lbas_cset_info, ibas_cset_info, 648 & lbas_cent_info, ibas_cent_info, 649 & ldocset, idocset, 650 & l_rchi_atom,i_rchi_atom, 651 & l_rq,i_rq,lniz, iniz, 652 & lchi_ao, ichi_ao, 653 & ldelchi_ao, idelchi_ao 654 double precision xyz_NMRcoords(3) 655 integer lhfine(3) ,ihfine(3) 656 integer inntsize,ddblsize,ok 657 logical grid_file_rewind 658 integer natoms 659 logical ofinite,Knucl 660 double precision zetanuc_arr(natoms) 661 external grid_file_rewind,calc_NMRHFine_F1ij, 662 & ga_antisymmetrize,grid_quadv0 663c 664c model potential parameters 665 character*2 gelem(ncenters) 666 double precision gexpo(ncenters,50) 667 double precision gcoef(ncenters,50) 668c 669c == allocate memory == 670 do i=1,3 671 if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao, 672 & 'lhfine',lhfine(i),ihfine(i))) 673 & call errquit('zora_getv_NMR: hfine',911,MA_ERR) 674 enddo ! end-loop-i 675c == preliminaries == 676 do j=1,3 677 do i= 1, nbf_ao*nbf_ao 678 dbl_mb(ihfine(j)+i-1) =0.d0 679 enddo 680 enddo 681c 682c get zora model potential parameters of geometry 683 if (use_modelpotential) 684 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 685c 686c == generate the grid == 687 dummy(1) = 0.d0 688 dummy(2) = 0.d0 689 call grid_quadv0(rtdb,g_dens,g_hfine(1),nexc,rho_n,dummy,tmat) 690c == ao basis set info used by xc_eval_basis == 691 if (.not.bas_numcont(AO_bas_han, ncontrset)) 692 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 693 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 694 & lbas_cent_info, ibas_cent_info)) 695 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 696 & MA_ERR) 697 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 698 & lbas_cset_info, ibas_cset_info)) 699 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 700 & MA_ERR) 701 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 702 & int_mb(ibas_cset_info), ncenters) 703 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 704 & ldocset, idocset)) 705 & call errquit('zora_getv_sf: cannot allocate ccdocset', 706 . ncontrset, MA_ERR) 707 do i=1,ncontrset 708 log_mb(idocset+i-1)=.true. 709 enddo 710 if(.not.MA_push_get(MT_int, ncenters, 'iniz', 711 & lniz, iniz)) 712 & call errquit("zora_getv_sf:iniz",0, MA_ERR) 713 do i= 1, ncenters 714 int_mb(iniz+i-1)=1 715 enddo 716 nxyz = 3*ncenters 717 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 718 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 719 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 720 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 721 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 722 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 723 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 724 & Dbl_MB(ixyz), Dbl_MB(icharge))) 725 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 726 727 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 728 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 729 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 730 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 731 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 732 & 'quad pts buffer', lrqbuf, irqbuf)) 733 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 734 735 if (.not. grid_file_rewind()) 736 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 737 & UNKNOWN_ERR) 738c 739c == loop over records in the grid file == 740 do iqsh = 1, n_rec_in_file 741c == define the current range of radial shells and integration center == 742 call grid_file_read(n_per_rec, nqpts, ictr_buf, 743 & rad,dbl_mb(irqbuf),nsubb) 744 if(nqpts.gt.buffer_size) 745 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 746c == loop over a subset of the grid == 747 istep=0 748 do ncube=1,nsubb 749c put buf into currently used arrays qxyz and qwght 750 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 751 & dbl_mb(iqwght), nqpts, rad,istep) 752 if(nqpts.ne.0) then 753c == compute the basis functions over the grid == 754 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 755 & l_rchi_atom,i_rchi_atom)) 756 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 757 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 758 & l_rq,i_rq)) 759 & call errquit("zora_getv_sf:rq",0, MA_ERR) 760c == delchi == 761 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 762 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 763 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 764c == chi == 765 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 766 & 'chi_ao', lchi_ao, ichi_ao)) 767 & call errquit('zora_getv: chi_ao',0, MA_ERR) 768 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 769 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 770 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 771 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 772 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 773 & int_mb(iniz), log_mb(idocset), 774 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 775 call calc_NMRHFine_F1ij(ao_bas_han,geom,ipol,g_dens, 776 & dbl_mb(idelchi_ao), 777 & dbl_mb(iqxyz),dbl_mb(iqwght), 778 & nbf_ao,nqpts,ncenters, 779 & ofinite, ! in: = .true. if Gaussian Nucl. Model of charges requested 780 & zetanuc_arr, ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model 781 & Knucl, 782 & use_modelpotential, 783 & gexpo, 784 & gcoef, 785 & dbl_mb(ihfine(1)), ! out 786 & dbl_mb(ihfine(2)), ! out 787 & dbl_mb(ihfine(3))) ! out 788c == delete memory == 789 if(.not.MA_pop_stack(lchi_ao)) 790 & call errquit("zora_getv: pop chi_ao", 100, MA_ERR) 791 if(.not.MA_pop_stack(ldelchi_ao)) 792 & call errquit("zora_getv: pop delchi_ao", 100, MA_ERR) 793 if(.not.MA_pop_stack(l_rq)) 794 & call errquit("zora_getv: pop rq", 100, MA_ERR) 795 if(.not.MA_pop_stack(l_rchi_atom)) 796 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 797 endif ! nqpts 798 enddo ! ncube 799 end do ! iqsh 800c == delete memory == 801 if(.not.MA_pop_stack(lrqbuf)) 802 & call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR) 803 if(.not.MA_pop_stack(lqwght)) 804 & call errquit("zora_getv_sf: pop qwght", 100, MA_ERR) 805 if(.not.MA_pop_stack(lqxyz)) 806 & call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR) 807 if(.not.MA_pop_stack(ltags)) 808 & call errquit("zora_getv_sf: pop tags", 100, MA_ERR) 809 if(.not.MA_pop_stack(lcharge)) 810 & call errquit("zora_getv_sf: pop charge", 100, MA_ERR) 811 if(.not.MA_pop_stack(lxyz)) 812 & call errquit("zora_getv_sf: pop xyz", 100, MA_ERR) 813 if(.not.MA_pop_stack(lniz)) 814 & call errquit("zora_getv_sf: pop niz", 100, MA_ERR) 815 if(.not.MA_pop_stack(ldocset)) 816 & call errquit("zora_getv_sf: pop docset", 100, MA_ERR) 817 if(.not.MA_pop_stack(lbas_cset_info)) 818 & call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR) 819 if(.not.MA_pop_stack(lbas_cent_info)) 820 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 821c 822c == tally up over all the nodes == 823 do k=1,3 824 call ga_dgop(msg_excrho,dbl_mb(ihfine(k)), nbf_ao*nbf_ao, '+') 825 enddo ! end-loop-k 826c 827c == pack into a ga all NMR AOs == 828 do i=1,3 829 call ga_zero(g_hfine(i)) 830 call ga_put(g_hfine(i), 831 & 1,nbf_ao,1,nbf_ao,dbl_mb(ihfine(i)),nbf_ao) 832 call ga_antisymmetrize(g_hfine(i)) 833 enddo ! end-loop-i 834 call ga_sync() 835c ----- free memory ----------- START 836 do i=1,3 837 if (.not.ma_free_heap(lhfine(i))) call errquit 838 & ('zora_getv_NMR: ma_free_heap of lhfine failed',911,MA_ERR) 839 enddo 840c ----- free memory ----------- END 841 return 842 end 843 logical function dft_mirrdens_start(g_dens,g_densrep,g_dens0, 844 i ipol) 845 implicit none 846#include "global.fh" 847 integer g_dens(2) ! [in/out] 848 integer g_dens0(2),g_densrep(2) ! [in/out] 849 integer ipol ! [in] 850c 851 integer dorep_glob,ii 852 logical dorepd,dorepl,docopy,dozero 853 logical util_mirrmat 854 external util_mirrmat 855 dorepl=.false. 856 call ga_sync() 857 if(ga_cluster_nnodes().gt.1) then 858 docopy=.true. 859 dozero=.false. 860 dorepd=util_mirrmat(ipol,g_dens,g_densrep,docopy,dozero) 861 dorep_glob=0 862 if(dorepd) dorep_glob=1 863 call ga_igop(375,dorep_glob,1, '+') 864 dorepl=dorep_glob.eq.ga_nnodes() 865 if(dorepl) then 866 do ii=1,ipol 867 g_dens0(ii)=g_dens(ii) 868 g_dens(ii)=g_densrep(ii) 869 enddo 870 else 871 if(dorepd) then 872 call util_mirrstop(g_densrep(1)) 873 if(ipol.eq.2) call util_mirrstop(g_densrep(2)) 874 endif 875 if(ga_nodeid().eq.0) then 876 write(6,*) ' no DM mirroring in zora_getv' 877 call util_flush(6) 878 endif 879 endif 880 endif 881 dft_mirrdens_start=dorepl 882 return 883 end 884