1 subroutine hnd_oechrg(ixyz,iprim,icoef,i_prim,i_gen,li, 2 & jxyz,jprim,jcoef,j_prim,j_gen,lj, 3 & kxyz,kprim,kcoef,k_prim,k_gen,lk, 4 & lxyz,lprim,lcoef,l_prim,l_gen,ll) 5c 6c $Id$ 7c 8 implicit none 9c 10#include "nwc_const.fh" 11#include "hnd_tol.fh" 12#include "hnd_giao.fh" 13c 14c Input parameters 15c 16 integer li,i_prim,i_gen 17 integer lj,j_prim,j_gen 18 integer lk,k_prim,k_gen 19 integer ll,l_prim,l_gen 20 double precision ixyz(3),iprim(i_prim),icoef(i_prim,i_gen) 21 double precision jxyz(3),jprim(j_prim),jcoef(j_prim,j_gen) 22 double precision kxyz(3),kprim(k_prim),kcoef(k_prim,k_gen) 23 double precision lxyz(3),lprim(l_prim),lcoef(l_prim,l_gen) 24c 25c Local variables 26c 27 integer ia,jb,jbmax 28 double precision ai,aj,arri,aa1,dum 29 double precision daexpa,cci,axi,ayi,azi,dtol,rtol,rri,aa,rrk 30c 31 rtol = 2.30258d00*itol 32 dtol = 1.0d+01**(-itol) 33c 34c ----- -ij- charge distribution ----- 35c 36 lit = li + 1 37 xi = ixyz(1) 38 yi = ixyz(2) 39 zi = ixyz(3) 40 ljt = lj + 1 41 xj = jxyz(1) 42 yj = jxyz(2) 43 zj = jxyz(3) 44 rri=((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 45 expndi = lit.ge.ljt 46 if (lit.ge.ljt) then 47 xc=xi 48 yc=yi 49 zc=zi 50 dxij=xi-xj 51 dyij=yi-yj 52 dzij=zi-zj 53 else 54 xc=xj 55 yc=yj 56 zc=zj 57 dxij=xj-xi 58 dyij=yj-yi 59 dzij=zj-zi 60 endif 61c 62c ----- -giao- factors ----- 63c 64 tijx = xi-xj 65 tijy = yi-yj 66 tijz = zi-zj 67 qijx = yi*zj-zi*yj 68 qijy = zi*xj-xi*zj 69 qijz = xi*yj-yi*xj 70c 71c ----- - i- primitive ----- 72c 73 nij=0 74 do 40 ia=1,i_prim 75 ai=iprim(ia) 76 arri=ai*rri 77 axi=ai*xi 78 ayi=ai*yi 79 azi=ai*zi 80 cci=icoef(ia,i_gen) 81c 82c ----- - j- primitive ----- 83c 84 jbmax=j_prim 85 if(iieqjj) jbmax=ia 86 do 30 jb=1,jbmax 87 aj=jprim(jb) 88 aa=ai+aj 89 aa1=1.0d0/aa 90 dum=aj*arri*aa1 91 if(dum.gt.rtol) go to 30 92 daexpa=cci*jcoef(jb,j_gen)* exp(-dum)*aa1 93 dum= abs(daexpa) 94 if(dum.le.dtol) go to 30 95c 96 nij=nij+1 97 acharg( 1,nij)= daexpa 98 if(iieqjj.and.ia.ne.jb) acharg( 1,nij)=2.0d0*daexpa 99 acharg( 2,nij)= aa 100 acharg( 3,nij)=(axi+aj*xj)*aa1 101 acharg( 4,nij)=(ayi+aj*yj)*aa1 102 acharg( 5,nij)=(azi+aj*zj)*aa1 103c 104 30 continue 105 40 continue 106c 107c ----- -kl- charge distribution ----- 108c 109 lkt = lk + 1 110 xk = kxyz(1) 111 yk = kxyz(2) 112 zk = kxyz(3) 113 lmt = ll + 1 114 xl = lxyz(1) 115 yl = lxyz(2) 116 zl = lxyz(3) 117 rri=((xk-xl)**2+(yk-yl)**2+(zk-zl)**2) 118 expndk = lkt.ge.lmt 119 if (lkt.ge.lmt) then 120 xd=xk 121 yd=yk 122 zd=zk 123 dxkl=xk-xl 124 dykl=yk-yl 125 dzkl=zk-zl 126 else 127 xd=xl 128 yd=yl 129 zd=zl 130 dxkl=xl-xk 131 dykl=yl-yk 132 dzkl=zl-zk 133 endif 134c 135c ----- -giao- factors ----- 136c 137 tklx = xk-xl 138 tkly = yk-yl 139 tklz = zk-zl 140 qklx = yk*zl-zk*yl 141 qkly = zk*xl-xk*zl 142 qklz = xk*yl-yk*xl 143c 144c ----- - k- primitive ----- 145c 146 nkl=0 147 do 60 ia=1,k_prim 148 ai=kprim(ia) 149 arri=ai*rri 150 axi=ai*xk 151 ayi=ai*yk 152 azi=ai*zk 153 cci=kcoef(ia,k_gen) 154c 155c ----- - l- primitive ----- 156c 157 jbmax=l_prim 158 if(kkeqll) jbmax=ia 159 do 50 jb=1,jbmax 160 aj=lprim(jb) 161 aa=ai+aj 162 aa1=1.0d0/aa 163 dum=aj*arri*aa1 164 if(dum.gt.rtol) go to 50 165 daexpa=cci*lcoef(jb,l_gen)* exp(-dum)*aa1 166 dum= abs(daexpa) 167 if(dum.le.dtol) go to 50 168c 169 nkl=nkl+1 170 bcharg( 1,nkl)= daexpa 171 if(kkeqll.and.ia.ne.jb) bcharg( 1,nkl)=2.0d0*daexpa 172 bcharg( 2,nkl)= aa 173 bcharg( 3,nkl)=(axi+aj*xl)*aa1 174 bcharg( 4,nkl)=(ayi+aj*yl)*aa1 175 bcharg( 5,nkl)=(azi+aj*zl)*aa1 176c 177 50 continue 178 60 continue 179 return 180 end 181