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