1C$PRAGMA SUN OPT=2
2      Subroutine hf1mke3(Axyz,Bxyz,Cxyz,alpha,G,GT,ABC2I,E,
3     &                   NABC,La,Lb,Lc)
4c $Id$
5
6      Implicit real*8 (a-h,o-z)
7      Implicit integer (i-n)
8
9      Dimension Axyz(3),Bxyz(3),Cxyz(3),alpha(4,NABC),G(3,NABC)
10
11      Dimension GT(NABC,3),ABC2I(3*NABC)
12
13      Dimension E(NABC,3,0:(La+Lb+Lc),0:La,0:Lb,0:Lc)
14c
15c Define the Hermite linear expansion coefficients for the product of
16c three monomials. These coefficients are scaled by a factor appropriate
17c for 3-ctr OIs, defined as
18c
19c         /    PI     \ 3/2     /   a b   __2 \      /   (a + b) c   __2 \
20c   ES = | ----------- |    EXP| - -----  AB   | EXP| - -----------  PC   |
21c         \ a + b + c /         \  a + b      /      \   a + b + c       /
22c
23c This scaling factor is passed as the 4th index of the exponents array
24c (i.e., "alpha(4,m)").
25c
26c     Recursion Formulae:
27c
28c       0,0,0
29c     Ex      = ES
30c       0
31c
32c       0,0,0
33c     Ey      = 1.0
34c       0
35c
36c       0,0,0
37c     Ez      = 1.0
38c       0
39c
40c       Ia+1,Ib,Ic           Ia,Ib,Ic         Ia,Ib,Ic            Ia,Ib,Ic
41c     Ex           = ABC2I Ex         + GAx Ex         + (Ip+1) Ex
42c       Ip                   Ip-1             Ip                  Ip+1
43c
44c       Ia,Ib+1,Ic           Ia,Ib,Ic         Ia,Ib,Ic            Ia,Ib,Ic
45c     Ex           = ABC2I Ex         + GBx Ex         + (Ip+1) Ex
46c       Ip                   Ip-1             Ip                  Ip+1
47c
48c       Ia,Ib,Ic+1           Ia,Ib,Ic         Ia,Ib,Ic            Ia,Ib,Ic
49c     Ex           = ABC2I Ex         + GCx Ex         + (Ip+1) Ex
50c       Ip                   Ip-1             Ip                  Ip+1
51c
52c     Indices for E(m,k,Ip,Ia,Ib,Ic):
53c
54c     1 [  m  ] NABC
55c     1 [  k  ] 3     {X,Y,Z}
56c     0 [  Ip ] Ia+Ib+Ic
57c     0 [  Ia ] La
58c     0 [  Ib ] Lb
59c     0 [  Ic ] Lc
60c
61c******************************************************************************
62
63c Initialize the Hermite expansion coefficients.
64
65      isz_e = NABC*3*(La+Lb+Lc+1)*(La+1)*(Lb+1)*(Lc+1)
66      call dcopy(isz_e,0d0,0,E,1)
67c Define E(Ip,Ia,Ib,Ic) for Ip=0, Ia=0, Ib=0, Ic=0.
68
69      do 100 m = 1,NABC
70
71       E(m,1,0,0,0,0) = alpha(4,m)
72       E(m,2,0,0,0,0) = 1.D0
73       E(m,3,0,0,0,0) = 1.D0
74
75       ABC2I(m       ) = 0.5D0/(alpha(1,m) + alpha(2,m) + alpha(3,m))
76       ABC2I(m+1*NABC) = ABC2I(m)
77       ABC2I(m+2*NABC) = ABC2I(m)
78
79       GT(m,1) = G(1,m) - Axyz(1)
80       GT(m,2) = G(2,m) - Axyz(2)
81       GT(m,3) = G(3,m) - Axyz(3)
82
83  100 continue
84
85c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=1,La, Ib=0, Ic=0.
86
87      do 250 Ia = 1,La
88
89c                          ===>   Ip = 0   <===
90
91       do 210 m = 1,3*NABC
92        E(m,1,0,Ia,0,0) =   GT(m,1)*E(m,1,0,Ia-1,0,0)
93     &                    +         E(m,1,1,Ia-1,0,0)
94  210  continue
95
96c                       ===>   Ip = 1,Ia-1   <===
97
98        do 230 Ip = 1,Ia-1
99
100         do 220 m = 1,3*NABC
101          E(m,1,Ip,Ia,0,0) =   ABC2I(m)*E(m,1,Ip-1,Ia-1,0,0)
102     &                       +  GT(m,1)*E(m,1,Ip  ,Ia-1,0,0)
103     &                       +   (Ip+1)*E(m,1,Ip+1,Ia-1,0,0)
104  220   continue
105
106  230  continue
107
108c                         ===>   Ip = Ia   <===
109
110       Ip = Ia
111
112       do 240 m = 1,3*NABC
113        E(m,1,Ip,Ia,0,0) =   ABC2I(m)*E(m,1,Ip-1,Ia-1,0,0)
114     &                     +  GT(m,1)*E(m,1,Ip  ,Ia-1,0,0)
115  240  continue
116
117  250 continue
118
119c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=0,La, Ib=1,Lb, Ic=0.
120
121      do 300 m = 1,NABC
122       GT(m,1) = G(1,m) - Bxyz(1)
123       GT(m,2) = G(2,m) - Bxyz(2)
124       GT(m,3) = G(3,m) - Bxyz(3)
125  300 continue
126
127      do 360 Ib = 1,Lb
128
129       do 350 Ia = 0,La
130
131c                          ===>   Ip = 0   <===
132
133        do 310 m = 1,3*NABC
134         E(m,1,0,Ia,Ib,0) =   GT(m,1)*E(m,1,0,Ia,Ib-1,0)
135     &                      +         E(m,1,1,Ia,Ib-1,0)
136  310   continue
137
138c                       ===>   Ip = 1,Ia+Ib-1   <===
139
140        do 330 Ip = 1,Ia+Ib-1
141
142         do 320 m = 1,3*NABC
143          E(m,1,Ip,Ia,Ib,0) =   ABC2I(m)*E(m,1,Ip-1,Ia,Ib-1,0)
144     &                        +  GT(m,1)*E(m,1,Ip  ,Ia,Ib-1,0)
145     &                        +   (Ip+1)*E(m,1,Ip+1,Ia,Ib-1,0)
146  320    continue
147
148  330   continue
149
150c                        ===>   Ip = Ia+Ib   <===
151
152        Ip = Ia + Ib
153
154        do 340 m = 1,3*NABC
155         E(m,1,Ip,Ia,Ib,0) =   ABC2I(m)*E(m,1,Ip-1,Ia,Ib-1,0)
156     &                       +  GT(m,1)*E(m,1,Ip  ,Ia,Ib-1,0)
157  340   continue
158
159  350  continue
160
161  360 continue
162
163c Define E(Ip,Ia,Ib,Ic) for Ip=0,Ia+Ib+Ic, Ia=0,La, Ib=0,Lb, Ic=1,Lc.
164
165      do 400 m = 1,NABC
166       GT(m,1) = G(1,m) - Cxyz(1)
167       GT(m,2) = G(2,m) - Cxyz(2)
168       GT(m,3) = G(3,m) - Cxyz(3)
169  400 continue
170
171      do 470 Ic = 1,Lc
172
173       do 460 Ib = 0,Lb
174
175        do 450 Ia = 0,La
176
177c                          ===>   Ip = 0   <===
178
179         do 410 m = 1,3*NABC
180          E(m,1,0,Ia,Ib,Ic) =   GT(m,1)*E(m,1,0,Ia,Ib,Ic-1)
181     &                        +         E(m,1,1,Ia,Ib,Ic-1)
182  410    continue
183
184c                       ===>   Ip = 1,Ia+Ib+Ic-1   <===
185
186         do 430 Ip = 1,Ia+Ib+Ic-1
187
188          do 420 m = 1,3*NABC
189           E(m,1,Ip,Ia,Ib,Ic) =   ABC2I(m)*E(m,1,Ip-1,Ia,Ib,Ic-1)
190     &                          +  GT(m,1)*E(m,1,Ip  ,Ia,Ib,Ic-1)
191     &                          +   (Ip+1)*E(m,1,Ip+1,Ia,Ib,Ic-1)
192  420     continue
193
194  430    continue
195
196c                        ===>   Ip = Ia+Ib+Ic   <===
197
198         Ip = Ia + Ib + Ic
199
200         do 440 m = 1,3*NABC
201          E(m,1,Ip,Ia,Ib,Ic) =   ABC2I(m)*E(m,1,Ip-1,Ia,Ib,Ic-1)
202     &                         +  GT(m,1)*E(m,1,Ip  ,Ia,Ib,Ic-1)
203  440    continue
204
205  450   continue
206
207  460  continue
208
209  470 continue
210
211      end
212