1 subroutine grid_signf(basis, mcenters, 2 & xyz,xyzm,ctr_quad, rad_quad, 3 & zprim, iandex,iatype_pt_chg, 4 & nq,xyzw) 5c 6C$Id$ 7c 8 implicit none 9#include "errquit.fh" 10c 11 integer basis 12 integer ctr_quad ! atom where the grid is centered [input] 13 integer mcenters ! [input] 14 double precision xyz(3,*) ! [input] 15 double precision xyzm(3,*) ! [input/output] 16 double precision rad_quad ! [input] 17 logical iatype_pt_chg(*) ! [input] 18 integer iandex(*) ! [scratch] 19 double precision xyzw(4,*) 20 integer nq 21 22c 23#include "bas.fh" 24c 25c Distance Squared between Sampling Points and Centers 26c 27 double precision zprim(*),acc_sigf,acc_sigf2 28 integer icset, ictr, 29 & ifirst, ilast, nprim, iprimo 30 double precision zmin,r2,x,y,z,r_arg 31 integer mcenters_scr,ctr_out 32 integer ic1,ic2,dum 33 integer n,l 34c 35 logical qpts_in 36c 37 double precision alpha,logeps,bfspread,bfspread2 38 double precision gaussian_range 39 gaussian_range(n,logeps,alpha) = 40 $ (n*log(-logeps) - n*log(alpha) - 4.0d0*logeps) / 41 $ sqrt(-16.0d0*alpha*logeps) 42 43c 44c 45c 46 acc_sigf=log(1d-10) 47 acc_sigf2=log(1d-13) 48 mcenters_scr = 0 49c 50 do ictr=1,mcenters 51 if(iatype_pt_chg(ictr)) goto 2001 52 if (ictr.eq.ctr_quad) then 53 mcenters_scr=mcenters_scr+1 54 iandex(mcenters_scr)=ictr 55 goto 2001 56 endif 57 if (.not.bas_ce2cnr(basis,ictr,ic1,ic2)) 58 & call errquit('Exiting in xc_signf.',11, BASIS_ERR) 59c 60c no bf on center ie but we keep grid 61c 62 if(ic1.eq.0) then 63 mcenters_scr=mcenters_scr+1 64 iandex(mcenters_scr)=ictr 65 goto 2001 66 endif 67 68 x = xyz(1,ctr_quad) - xyz(1,ictr) 69 y = xyz(2,ctr_quad) - xyz(2,ictr) 70 z = xyz(3,ctr_quad) - xyz(3,ictr) 71 r2 = sqrt(x*x + y*y + z*z) 72 r_arg=0d0 73 if (rad_quad.lt.r2) r_arg = (r2-rad_quad) 74 do icset = ic1,ic2 75c 76c get info about current contraction set 77c 78 if (.not. bas_cn2bfr(basis, icset, ifirst, ilast)) 79 & call errquit('Exiting in grid_signf.',4, BASIS_ERR) 80 if (.not. bas_continfo(basis, icset, l ,nprim, 81 & dum, dum)) 82 & call errquit('Exiting in grid_signf.',4, BASIS_ERR) 83c 84c 85c get exponents and contraction coefficients for this contraction set 86c 87 if (.not.bas_get_exponent(basis, icset, zprim)) 88 & call errquit('Exiting in grid_signf.',7, BASIS_ERR) 89c 90c Determine the minimum Gaussian exponent. 91c 92 93 zmin = 1.D+06 94 do iprimo = 1,nprim 95 zmin = min(zprim(iprimo),zmin) 96 enddo 97c 98c Only include those basis functions that are "non-zero" for at least 99c one point in the sampling set. 100c 101 bfspread=gaussian_range(l,acc_sigf,zmin) 102 if (r_arg.lt.bfspread) then 103#ifdef GRID_ASCREEN 104 bfspread2=gaussian_range(l,acc_sigf2,zmin) 105c 106c check if all grid pts are really in the bf spread 107c 108 do iprimo=1,nq 109 x = xyzw(1,iprimo) - xyz(1,ictr) 110 y = xyzw(2,iprimo) - xyz(2,ictr) 111 z = xyzw(3,iprimo) - xyz(3,ictr) 112 r2 = sqrt(x*x + y*y + z*z) 113 if (r2.lt.bfspread2) then 114#endif 115 mcenters_scr=mcenters_scr+1 116 iandex(mcenters_scr)=ictr 117 goto 2001 118#ifdef GRID_ASCREEN 119 endif 120 enddo 121#endif 122 endif 123 124 enddo 125 2001 continue 126 enddo 127 ctr_out = 1 ! take care of compiler warnings 128 do ictr = 1, mcenters_scr 129 ic1=iandex(ictr) 130 xyzm(1,ictr) = xyzm(1,ic1) 131 xyzm(2,ictr) = xyzm(2,ic1) 132 xyzm(3,ictr) = xyzm(3,ic1) 133 if(ic1.eq.ctr_quad) ctr_out=ictr 134 enddo 135 mcenters=mcenters_scr 136 ctr_quad=ctr_out 137 return 138 end 139