1 subroutine dim_grid_quadv0b(ictr_buf, qxyz, qwght, nqpts, rad, 2 $ ncoef_max, natoms, iniz, ifin, rho, 3 $ xyz, expo, Bmat, Fmat, Pmat, rq, 4 $ cetobfr, ibf_ao, bas_cent_info, 5 $ bas_cset_info, docset, rchi_atom, 6 $ g_vdim, nmat, do_scf, do_cpks_l, l3d, 7 $ ipm, imag, dimxyz, muind) 8c 9c$Id: dim_grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $ 10c 11 implicit none 12#include "errquit.fh" 13#include "cgridfile.fh" 14#include "stdio.fh" 15#include "cdft.fh" 16#include "bas.fh" 17#include "mafdecls.fh" 18#include "global.fh" 19#include "util.fh" 20#include "grid_sic.fh" 21#include "dftps.fh" 22#include "geom.fh" 23#include "dimqm.fh" 24c 25c Input Variables 26 integer ictr_buf ! ctr of grid 27 double precision qxyz(3, n_per_rec) ! Quadrature point coordinates 28 double precision qwght(n_per_rec) ! Quadrature point weights 29 integer nqpts ! Number of quadrature points 30 double precision rad 31 integer ncoef_max 32 integer natoms ! Number of atoms 33 integer iniz(natoms) 34 integer ifin(natoms) 35 double precision rho(*) 36 double precision xyz(3, natoms) ! Atom coordinates 37 double precision expo(nbf_ao_mxprim) 38 double precision Bmat(nbf_ao_mxnbf_ce*n_per_rec) 39 double precision Fmat(*) 40 double precision Pmat(*) 41 double precision rq(n_per_rec, natoms) 42 integer cetobfr(2, natoms) 43 integer ibf_ao(nbf_ao) 44 integer bas_cent_info(3,natoms) 45 integer bas_cset_info(6,*) 46 logical docset(*) 47 double precision rchi_atom(natoms) 48 integer g_vdim ! DIM potential for this block of the quadrature 49 integer nmat ! Number of density matrices (number of perturbing directions) 50 logical do_scf 51 logical do_cpks_l 52 logical l3d 53 integer ipm 54 integer imag 55 double precision dimxyz(3, nDIM) ! DIM Coordinates 56 double precision muind(3, nDIM, nmat) ! DIM Coordinates 57c 58c Local Variables 59 integer nbf_ao_mxnbf_ce2, mbf_ao, npol, ncontrset, maxdbas 60 integer lchi_ao, ichi_ao, ldelchi_ao, idelchi_ao 61 integer ldmat,idmat,i, dim_grid_nbfm 62 double precision dabsmax 63 double precision acc_ao_gauss 64 external dabsmax 65 external dim_grid_nbfm 66 integer iscf_rho,iscf_delrho 67 integer iscf_tau,iscf_ttau 68 logical grid_sicinit,grid_sicend 69 external grid_sicinit,grid_sicend 70 double precision acc_xc_gauss 71 integer lemat,iemat,lfmat,ifmat,k_scr,l_scr 72 logical do_2nd 73 logical stat 74 logical dprint 75c 76 npol = (ipol*(ipol+1))/2 77 acc_ao_gauss = dble(iaoacc) 78 nbf_ao_mxnbf_ce2 = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce 79 if(nqpts.eq.0) return 80 dprint = ldebug 81c if(dprint) then 82c write(luout,*) "dim_grid_quadv0b start" 83c end if 84 maxdbas = 0 85 idelchi_ao = 0 86c 87c Evaluate the AO basis set at each of the quad. points. 88c allocate arrays for exponents and contraction coefficients 89c within int_eval_gbsets 90c Evaluate the distances (squared) between each center and the points 91c 92 mbf_ao=nbf_ao 93 if(natoms.gt.1) then 94 call icopy(mbf_ao, 0,0, ibf_ao, 1) 95 96 mbf_ao=dim_grid_nbfm(AO_bas_han, natoms, 97 & ictr_buf,rad,xyz,qxyz,nqpts, 98 & ibf_ao, docset, iniz, ifin, expo, 99 & minexp,ldiff,acc_ao_gauss,iatype_pt_chg) 100 if (mbf_ao.eq.0) return 101 else 102 if (.not.bas_numcont(AO_bas_han, ncontrset)) 103 & call errquit('Exiting in grid_quadv0b',0, BASIS_ERR) 104 iniz(1)=1 105 ifin(1)=nbf_ao 106 do i=1,ncontrset 107 docset(i)=.true. 108 enddo 109 do i=1,nbf_ao 110 ibf_ao(i)=i 111 enddo 112 endif 113 114 if (.not.MA_Push_Get(mt_dbl, nqpts*mbf_ao, 'chi_ao', 115 & lchi_ao, ichi_ao)) 116 & call errquit('grid_quadv0b: cannot allocate chi_ao',0, 117 & MA_ERR) 118 119 call qdist(rchi_atom, rq, qxyz, xyz, nqpts, natoms) 120 121 call xc_eval_basis(ao_bas_han, maxdbas, dbl_mb(ichi_ao), 122 & 0d0, 0d0, 0d0, rq, qxyz, xyz, nqpts, natoms, 123 & iniz, docset, bas_cent_info, bas_cset_info) 124c 125c get reduced Xi(r) and dXi(r) over atoms 126c 127 call util_rmsatbf(nqpts, natoms,iniz,ifin, 128 , dbl_mb(ichi_ao),rchi_atom) 129 130 if(nqpts.eq.0) goto 2010 131 132c 133c Evaluate the DIM potential 134c Set up pointer to the SCF density for the CPKS LHS case. 135c 136 call dim_eval_fnl(rho, nqpts, qwght, 137 & dimxyz, qxyz, 138 $ nmat, ipm, imag, muind) 139c rho(1:nqpts*nmat) = 0.0d0 140c 141c Transform the potential from over the quadrature points to the shape of the Fock Matrix 142 if(do_scf) then 143c 144 call dim_tabcd(0, l3d, Fmat, Pmat, rho, Bmat, 145 $ dbl_mb(ichi_ao), 1, nqpts, mbf_ao, 146 $ nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, 147 $ AO_bas_han, natoms, iniz, ifin, g_vdim, 148 $ ibf_ao, cetobfr) 149c 150 else if (do_cpks_l) then 151c 152 call dim_tabcd(1, l3d, Fmat, Pmat, rho, Bmat, 153 $ dbl_mb(ichi_ao), nmat, nqpts, mbf_ao, 154 $ nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, 155 $ AO_bas_han, natoms, iniz, ifin, g_vdim, 156 $ ibf_ao, cetobfr) 157c 158 endif 159c write(luout,*) "g_dim AFTER dim_tabcd" 160c call ga_print(g_vxc) 161c 162c endif 163 164 2010 continue 165 166 if (sic_orb_index.eq.1) then 167 if(.not.grid_sicend(l_vect1,ldelrho_sig)) 168 , call errquit(' grid_quadv0b: grid_sicend failed',0, 169 & CALC_ERR) 170 endif 171 172 if (.not.ma_pop_stack(lchi_ao)) 173 & call errquit('grid_quadv0b: cannot pop stack', 3, MA_ERR) 174 175c if(dprint) then 176c write(luout,*) "end dim_grid_quadv0b" 177c end if 178 return 179 end subroutine dim_grid_quadv0b 180c 181c function dim_grid_nbfm 182c 183 integer function dim_grid_nbfm(basis_hand, mcenters, 184 & ctr_quad, rad_quad, xyz, qxyz,nq, 185 & ibf, docset, iniz, ifin, zprim, 186 . minexp,ldiff,acc_gauss,iatype_pt_chg) 187c 188C$Id: grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $ 189c 190 implicit none 191#include "errquit.fh" 192c 193 integer basis_hand 194 integer mcenters ! [input] 195 double precision acc_gauss ! [input] 196 double precision xyz(3,*) ! [input] 197 double precision rad_quad ! something about radius of this integration shell [input] 198 integer nq ! [in] # grid pts 199 double precision qxyz(3,nq) ! [input] coord grid pts 200 integer ctr_quad ! grid center [input] 201 integer ibf(*) ! [output] 202 logical docset(*) ! [output] 203 integer ldiff(*) ! [in] 204 double precision minexp(*) ! [in] 205 logical iatype_pt_chg(*) ! [in] 206c 207#include "bas.fh" 208c 209c Distance Squared between Sampling Points and Centers 210c 211 double precision r_q0 ! min dist 212 integer iniz(mcenters),ifin(mcenters) 213 double precision zprim(*) 214 integer n1, icset, ictr, nprim, ncontr, 215 & isphere, l, iprimo,npt,nshbf 216 double precision zmin,acc_loc,x,y,z 217 integer n,ic1,ic2,m,ibfhi,ibflo 218 double precision alpha,logeps 219 double precision gaussian_range,r2,r_arg 220 logical veryfar 221c 222 gaussian_range(n,logeps,alpha) = 223 $ (n*log(-logeps) - n*log(alpha) - 4.0d0*logeps) / 224 $ sqrt(-16.0d0*alpha*logeps) 225c 226 call ifill(mcenters,0,iniz,1) 227 call ifill(mcenters,0,ifin,1) 228 if(acc_gauss.gt.25d0.or.acc_gauss.lt.3d0) call errquit( 229 ' ' grid_nbfm: silly accgauss ',nint(acc_gauss), UNKNOWN_ERR) 230 acc_loc=-acc_gauss 231c 232 n1 = 0 233 npt=0 234 do 400 ictr=1,mcenters ! Loop over all atoms 235 if(iatype_pt_chg(ictr)) goto 400 ! Cycle if this atom is a point charge 236 if (.not.bas_ce2cnr(basis_hand,ictr,ic1,ic2)) 237 & call errquit('Exiting in xc_signf.',11, BASIS_ERR) 238 239 r2=0d0 240 if (ictr.ne.ctr_quad) then ! Calculate r-squared if this atom isn't the center of the quadrature 241 x = xyz(1,ctr_quad) - xyz(1,ictr) 242 y = xyz(2,ctr_quad) - xyz(2,ictr) 243 z = xyz(3,ctr_quad) - xyz(3,ictr) 244 r2 = sqrt(x*x + y*y + z*z) 245 endif 246 r_arg=0d0 247 if (rad_quad.lt.r2) r_arg = (r2-rad_quad) 248c 249c check on most diffuse fn 250c 251 veryfar=r_arg.gt.gaussian_range(ldiff(ictr),acc_loc, 252 . minexp(ictr)) 253 if(veryfar) then 254 if(.not.bas_ce2bfr(basis_hand, ictr, ibflo, ibfhi)) 255 & call errquit('Exiting in grid_nbfm',4, BASIS_ERR) 256 nshbf=ibfhi-ibflo+1 257 npt=npt+nshbf 258 else 259 r_q0=1d10 260 do n=1,nq 261 x = qxyz(1,n) - xyz(1,ictr) 262 y = qxyz(2,n) - xyz(2,ictr) 263 z = qxyz(3,n) - xyz(3,ictr) 264 r_q0 = min(r_q0,sqrt(x*x + y*y + z*z)) 265 enddo 266 do icset = ic1,ic2 267 docset(icset) = .false. 268c 269c get info about current contraction set 270 if (.not. bas_continfo(basis_hand, icset, l ,nprim, 271 & ncontr, isphere)) 272 & call errquit('Exiting in grid_nbfm.',5, BASIS_ERR) 273c 274c get exponents and contraction coefficients for this contraction set 275 if (.not.bas_get_exponent(basis_hand, icset, zprim)) 276 & call errquit('Exiting in grid_nbfm.',7, BASIS_ERR) 277c 278c Determine the minimum Gaussian exponent. 279 zmin = 1.D+06 280 do iprimo = 1,nprim 281 zmin = min(zprim(iprimo),zmin) 282 enddo 283c 284c Only include those basis functions that are "non-zero" for at least 285c one point in the sampling set. 286 nshbf=ncontr*(((l+1)*(l+2))/2) 287 if(isphere.eq.1) then 288 nshbf=ncontr*(2*l+1) 289 endif 290c 291c pre-screening for worst case (max radius) 292c 293 if (r_q0.lt.gaussian_range(l,acc_loc,zmin)) then 294 do m=1,nshbf 295 ibf(n1+m) = npt+m 296 enddo 297 docset(icset) = .true. 298 if (iniz(ictr).eq.0) iniz(ictr)=n1+1 299 n1=n1+nshbf 300 endif 301 npt=npt+nshbf 302 enddo ! icset = ic1, ic2 303 ifin(ictr)= n1 304 endif ! If veryfar 305 400 continue ! End loop over atoms 306c 307 dim_grid_nbfm = n1 308c 309 return 310 end 311c 312c Precompute relevant basis info for XC 313c 314c BGJ - 9/00 315c 316 Subroutine dim_make_basis_info(basis_hand, bas_cent_info, 317 & bas_cset_info, mcenters) 318c 319C$Id: grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $ 320c 321 implicit none 322#include "errquit.fh" 323c 324 integer basis_hand ! [input] 325 integer mcenters ! [input] 326 integer bas_cent_info(3,mcenters) ! [output] 327 integer bas_cset_info(6,*) ! [output] 328c 329#include "bas.fh" 330c 331 integer ictr, icset, ifirst, ilast, nprim, ncontr, l, isphere 332 integer ic1, ic2 333c 334 do ictr = 1,mcenters 335 bas_cent_info(1,ictr) = 0 ! max angular momentum 336 if (.not.bas_ce2cnr(basis_hand,ictr,ic1,ic2)) 337 & call errquit('Exiting in xc_make_basis_info',1, BASIS_ERR) 338 bas_cent_info(2,ictr) = ic1 339 bas_cent_info(3,ictr) = ic2 340c 341 do icset = ic1,ic2 342c 343c Get info about current contraction set - first and last basis function, 344c angular momentum, number of primitives, number of contractions and 345c whether spherical harmomic 346c 347 if (.not. bas_cn2bfr(basis_hand, icset, ifirst, ilast)) 348 & call errquit('Exiting in xc_make_basis_info',1, 349 & BASIS_ERR) 350 if (.not. bas_continfo(basis_hand, icset, l, nprim, 351 & ncontr, isphere)) 352 & call errquit('Exiting in xc_make_basis_info',2, 353 & BASIS_ERR) 354c 355 if (l.lt.0) then 356 call errquit('L code < 0 in xc_make_basis_info',1, 357 & BASIS_ERR) 358 endif 359c 360 bas_cent_info(1,ictr) = max(bas_cent_info(1,ictr),l) 361c 362 bas_cset_info(1,icset) = ifirst 363 bas_cset_info(2,icset) = ilast 364 bas_cset_info(3,icset) = l 365 bas_cset_info(4,icset) = nprim 366 bas_cset_info(5,icset) = ncontr 367 bas_cset_info(6,icset) = isphere 368 enddo 369 enddo 370c 371 return 372 end 373