1c 2c == EFGZ4 scalar and spin-orbit components (x,y,z) 3 subroutine zora_getv_EFGZ4_SO(rtdb, 4 & g_dens, ! in : atomic density 5 & zora_calc_type,! in : =4 EFG-NUM, =3 EFGZ4-SR 6 & zora_Qpq, ! in : type of EFG potential 7 & xyz_EFGcoords, ! in : EFG-nuclear coordinates 8 & g_efgz4_sf, ! out: munu matrix 9 & g_efgz4_so, ! out: munu matrix 10 & nexc) 11c 12C$Id$ 13c Adapted from zora_getv_so 14 15 implicit none 16#include "rtdb.fh" 17#include "bas.fh" 18#include "cdft.fh" 19#include "errquit.fh" 20#include "mafdecls.fh" 21#include "global.fh" 22#include "geom.fh" 23#include "msgtypesf.h" 24#include "msgids.fh" 25#include "stdio.fh" 26#include "cgridfile.fh" 27#include "grid_cube.fh" 28#include "modelpotential.fh" 29c 30c == arguments == 31 integer rtdb 32 integer g_dens(2) 33 integer g_efgz4_sf(2) 34 integer g_efgz4_so(3) 35 integer zora_Qpq,zora_calc_type 36 double precision xyz_EFGcoords(3) 37 integer nexc 38c 39c == local variables == 40 integer i,j,k,ind,nij 41 double precision rho_n 42 double precision tmat 43 double precision dummy(2) 44 integer iqsh, istep, nxyz, ncontrset 45 integer ixyz, lxyz, icharge, lcharge, itags, ltags 46 integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube, 47 & ictr_buf,iqpts 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 & lscal0, iscal0, 57 & lscalx, iscalx, 58 & lscaly, iscaly, 59 & lscalz, iscalz 60 integer lscal(4),iscal(4) 61 integer inntsize,ddblsize,ok 62 logical grid_file_rewind 63 external grid_file_rewind,calc_zora_EFGZ4_SO,calc_EFG, 64 & ga_antisymmetrize 65c 66c model potential parameters 67 character*2 gelem(ncenters) 68 double precision gexpo(ncenters,50) 69 double precision gcoef(ncenters,50) 70c 71c == allocate memory == 72 do i=1,4 73 if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao, 74 & 'scal0',lscal(i),iscal(i))) 75 & call errquit('zora_getv_so: scali',0, MA_ERR) 76 enddo 77c == preliminaries == 78 do i= 1, nbf_ao*nbf_ao 79 do j=1,4 80 dbl_mb(iscal(j)+i-1)=0.d0 81 enddo 82 enddo 83c 84c get zora model potential parameters of geometry 85 if (use_modelpotential) 86 & call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef) 87c 88c == generate the grid == 89 dummy(1) = 0.d0 90 dummy(2) = 0.d0 91 call grid_quadv0(rtdb,g_dens,g_efgz4_sf,nexc,rho_n,dummy,tmat) 92c == ao basis set info used by xc_eval_basis == 93 if (.not.bas_numcont(AO_bas_han, ncontrset)) 94 & call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR) 95 if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info', 96 & lbas_cent_info, ibas_cent_info)) 97 & call errquit('zora_getv_sf: cannot allocate bas_cent_info',0, 98 & MA_ERR) 99 if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info', 100 & lbas_cset_info, ibas_cset_info)) 101 & call errquit('zora_getv_sf: cannot allocate bas_cset_info',0, 102 & MA_ERR) 103 call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info), 104 & int_mb(ibas_cset_info), ncenters) 105 if (.not.MA_Push_Get(mt_log, ncontrset, 'docset', 106 & ldocset, idocset)) 107 & call errquit('zora_getv_sf: cannot allocate ccdocset', 108 . ncontrset, MA_ERR) 109 do i=1,ncontrset 110 log_mb(idocset+i-1)=.true. 111 enddo 112c 113 if(.not.MA_push_get(MT_int, ncenters, 'iniz', 114 & lniz, iniz)) 115 & call errquit("zora_getv_sf:iniz",0, MA_ERR) 116 do i= 1, ncenters 117 int_mb(iniz+i-1)=1 118 enddo 119c 120 nxyz = 3*ncenters 121 if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz)) 122 & call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR) 123 if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge)) 124 & call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR) 125 if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags)) 126 & call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR) 127 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 128 & Dbl_MB(ixyz), Dbl_MB(icharge))) 129 & call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR) 130 131 if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz)) 132 & call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR) 133 if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght)) 134 & call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR) 135 if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4, 136 & 'quad pts buffer', lrqbuf, irqbuf)) 137 & call errquit('zora_getv_sf: quad buffer', 3, MA_ERR) 138 139 if (.not. grid_file_rewind()) 140 $ call errquit('zora_getv_sf: rewinding gridpts?', 0, 141 & UNKNOWN_ERR) 142c 143c == loop over records in the grid file == 144 do iqsh = 1, n_rec_in_file 145c 146c == define the current range of radial shells and integration center == 147 call grid_file_read(n_per_rec, nqpts, ictr_buf, 148 & rad,dbl_mb(irqbuf),nsubb) 149 150 if(nqpts.gt.buffer_size) 151 & call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR) 152c 153c == loop over a subset of the grid == 154 istep=0 155 do ncube=1,nsubb 156c 157c put buf into currently used arrays qxyz and qwght 158 call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz), 159 & dbl_mb(iqwght), nqpts, rad,istep) 160 161 if(nqpts.ne.0) then 162c 163c == compute the basis functions over the grid == 164 if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom', 165 & l_rchi_atom,i_rchi_atom)) 166 & call errquit("zora_getv:rchi_atom",0, MA_ERR) 167c 168 if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq', 169 & l_rq,i_rq)) 170 & call errquit("zora_getv_sf:rq",0, MA_ERR) 171c 172c == delchi == 173 if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao, 174 & 'delchi_ao', ldelchi_ao, idelchi_ao)) 175 & call errquit('zora_getv: delchi_ao',0, MA_ERR) 176c 177c == chi == 178 if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao, 179 & 'chi_ao', lchi_ao, ichi_ao)) 180 & call errquit('zora_getv: chi_ao',0, MA_ERR) 181 call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq), 182 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters) 183 call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao), 184 & dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq), 185 & dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters, 186 & int_mb(iniz), log_mb(idocset), 187 & int_mb(ibas_cent_info), int_mb(ibas_cset_info)) 188 if (zora_calc_type.eq.3) then 189c write(*,*) 'enter calc_zora_EFGZ4_so() ...' 190 call calc_zora_EFGZ4_so(ao_bas_han,geom,ipol,g_dens, 191 & dbl_mb(idelchi_ao), 192 & dbl_mb(iqxyz),dbl_mb(iqwght), 193 & nbf_ao,nqpts,ncenters, 194 & zora_Qpq,xyz_EFGcoords, 195 & use_modelpotential, 196 & gexpo, gcoef, 197 & dbl_mb(iscal(1)), ! out 198 & dbl_mb(iscal(2)), ! out 199 & dbl_mb(iscal(3)), ! out 200 & dbl_mb(iscal(4))) ! out 201 else if (zora_calc_type.eq.4) then 202c write(*,*) 'enter calc_EFG() ...' 203 call calc_EFG(ao_bas_han,geom,ipol,g_dens, 204 & dbl_mb(ichi_ao), 205 & dbl_mb(iqxyz),dbl_mb(iqwght), 206 & nbf_ao,nqpts,ncenters, 207 & zora_Qpq,xyz_EFGcoords, 208 & dbl_mb(iscal(1))) ! out 209 endif 210c == delete memory == 211 if(.not.MA_chop_stack(l_rchi_atom)) 212 & call errquit("zora_getv: pop rchi_atom",100,MA_ERR) 213 endif ! nqpts 214 enddo ! ncube 215 end do ! iqsh 216c 217c == delete memory == 218 if(.not.MA_chop_stack(lbas_cent_info)) 219 & call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR) 220c 221c == tally up over all the nodes == 222 do k=1,4 223 call ga_dgop(msg_excrho, dbl_mb(iscal(k)), nbf_ao*nbf_ao, '+') 224 enddo 225c 226c == pack into a ga == 227c 228c == scalar contribution == 229 do i=1,ipol 230 call ga_zero(g_efgz4_sf(i)) 231 call ga_put(g_efgz4_sf(i),1,nbf_ao,1,nbf_ao,dbl_mb(iscal(1)) 232 & ,nbf_ao) 233 call ga_symmetrize(g_efgz4_sf(i)) 234 enddo 235c == spin-orbit contributions == 236 ind=4 ! = 4,3,2 = z,y,x in g_zora_so(i) 237 do i=1,3 ! = z,y,x 238 call ga_zero(g_efgz4_so(i)) 239 call ga_put(g_efgz4_so(i), 240 & 1,nbf_ao,1,nbf_ao,dbl_mb(iscal(ind)),nbf_ao) 241 call ga_antisymmetrize(g_efgz4_so(i)) 242 ind=ind-1 243 enddo 244 if (.not.MA_chop_stack(lscal(1))) 245 & call errquit('zora_getv_so: scali chop_stack',0, MA_ERR) 246 call ga_sync() 247 return 248 end 249