1c $Id$
2      Subroutine hferi_gen(Ep,Eq,PairP,PairQ,R0,IJK,ERI,E3,t1,t2,t3,t4,
3     &    MPP,MPQ,NPP,NPQ,La,Lb,Lc,Ld,La2,Lb2,Lc2,Ld2,Lqmax,Lqmax3,Lr,
4     &    Acoefs,Bcoefs,Ccoefs,Dcoefs,NPA,NPB,NPC,NPD,NCA,NCB,NCC,NCD,
5     &    MXD,canAB,canCD,canPQ)
6
7      Implicit none
8#include "errquit.fh"
9
10      Logical canAB,canCD,canPQ
11
12c--> Number of pairs in pass, total number
13
14      Integer MPP,MPQ, NPP,NPQ
15
16c--> Angular momenta and number of cartesian components
17
18      Integer La,Lb,Lc,Ld, La2,Lb2,Lc2,Ld2, Lqmax,Lqmax3,Lr
19
20c--> Derivative order
21
22      Integer MXD
23
24c--> Number of primitives, number of contracted functions
25
26      Integer NPA,NPB,NPC,NPD, NCA,NCB,NCC,NCD
27
28c--> Hermite Linear Expansion Coefficients
29
30      Double precision Ep(3,NPP,0:MXD,0:(La+Lb),0:La,0:Lb)
31      Double precision Eq(3,NPQ,0:MXD,0:(Lc+Ld),0:Lc,0:Ld)
32
33c--> Pair index arrays
34
35      Integer PairP(2,NPP)
36      Integer PairQ(2,NPQ)
37
38c--> Auxiliary Function Integrals & Index
39
40      Double precision R0(MPQ,MPP,*)
41      Integer IJK(0:Lr,0:Lr,0:Lr)
42
43c--> Contraction coefficients
44
45      Double precision Acoefs(NPA,NCA)
46      Double precision Bcoefs(NPB,NCB)
47      Double precision Ccoefs(NPC,NCC)
48      Double precision Dcoefs(NPD,NCD)
49
50c--> ERI
51
52      Double precision ERI(Ld2,NCD,Lc2,NCC,Lb2*NCB,La2*NCA)
53
54c--> Scratch Space
55
56      Double precision E3(*)
57      Double precision t1(MPQ,Lqmax3,NPB,NCA)  ! 1st 1/4 trans
58      Double precision t2(MPQ,Lqmax3,MPP)      ! 2nd 1/4 trans, primitives
59      Double precision t3(MPQ)                 ! 3rd 1/4 trans
60      Double precision t4(NPD)                 ! 4th 1/4 trans
61      Integer Nxyz(3)
62
63c--> Local variables
64
65      Integer Lq,ninti
66      Integer ma,mb,mc,md, mp,mq, nq,nr
67      Integer Ia,Ja,Ka, Ib,Jb,Kb, Ic,Jc,Kc, Id,Jd,Kd
68      Integer Ip,Jp,Kp, Iq,Jq,Kq, Ir,Jr,Kr
69      Integer ipa,ipb,ipc,ipd, ica,icb,icc,icd
70      Integer ie,icab,NCAB
71      Integer ind_ca,ind_cb
72      Double precision ca,cb
73
74c
75c Compute electron repulsion integrals (ERI).
76c
77c     Formula:
78c
79c               __
80c               \     Ic,Id;n10   Jc,Jd;n11   Kc,Kd;n12
81c     ERI  =    /   Ex          Ey          Ez          SUM
82c               --    Iq          Jq          Kq           Iq,Jq,Kq
83c            Iq,Jq,Kq
84c
85c                            __
86c                           \     Lq   Ia,Ib;n7   Ja,Jb;n8   Ka,Kb;n9
87c         SUM          =    /  (-1)  Ex         Ey         Ez         R
88c            Iq,Jq,Kq       --         Ip         Jp         Kp        Ir,Jr,Kr
89c                        Ip,Jp,Kp
90c
91c                                Ir = (Ip+n1) + (Iq+n4)
92c                        where   Jr = (Jp+n2) + (Jq+n5)
93c                                Kr = (Kp+n3) + (Kq+n6)
94c
95c                         and    Lq = (Iq+n4) + (Jq+n5) + (Kq+n6)
96c
97c N.B.  For simple ERI (i.e., no derivative integrals)  n[1-12] = 0!
98c
99c******************************************************************************
100c General case:  [ab|cd]
101
102      if (canAB .or. canCD .or. canPQ)
103     &    call errquit ('hferi_gen:can''t do canonical integrals',911,
104     &       INT_ERR)
105
106c Initialize the block of ERIs.
107
108      NCAB = NCA*NCB
109      ninti = La2*Lb2*Lc2*Ld2*NCA*NCB*NCC*NCD
110      call dfill (ninti,0.0d0,eri,1)
111
112c Loop over the components of the "A" and "B" shells.
113
114      do ma = 1,La2
115
116        call getNxyz(La,ma,Nxyz)
117        Ia = Nxyz(1)
118        Ja = Nxyz(2)
119        Ka = Nxyz(3)
120
121        do mb = 1,Lb2
122
123          call getNxyz(Lb,mb,Nxyz)
124          Ib = Nxyz(1)
125          Jb = Nxyz(2)
126          Kb = Nxyz(3)
127
128c Sum across (Ip,Jp,Kp) for each value of (Iq,Jq,Kq).
129
130          call dfill(MPP*MPQ*Lqmax3,0.0d00,t2,1)
131
132          do Ip = 0,Ia+Ib
133            do Jp = 0,Ja+Jb
134              do Kp = 0,Ka+Kb
135
136c Define the product of the Hermite expansion coefficients for
137c overlap distribution "P".
138
139                do mp = 1,MPP
140                  E3(mp) = Ep(1,mp,0,Ip,Ia,Ib)
141     &                    *Ep(2,mp,0,Jp,Ja,Jb)
142     &                    *Ep(3,mp,0,Kp,Ka,Kb)
143                end do
144
145                do Iq = 0,Lqmax
146                  do Jq = 0,Lqmax-Iq
147                    do Kq = 0,Lqmax-Iq-Jq
148
149                      nq = IJK(Iq,Jq,Kq)
150                      Ir = Ip + Iq
151                      Jr = Jp + Jq
152                      Kr = Kp + Kq
153                      nr = IJK(Ir,Jr,Kr)
154
155c Include the factor of (-1)**(Iq+Jq+Kq).
156c Sum over Hermite functions
157
158                      Lq = Iq + Jq + Kq
159                      if( mod(Lq,2).eq.0 )then
160                        do mp = 1,MPP
161                          do mq = 1,MPQ
162                            t2(mq,nq,mp) = t2(mq,nq,mp)
163     &                          + E3(mp)*R0(mq,mp,nr)
164                          end do
165                        end do
166                      else
167                        do mp = 1,MPP
168                          do mq = 1,MPQ
169                            t2(mq,nq,mp) = t2(mq,nq,mp)
170     &                          - E3(mp)*R0(mq,mp,nr)
171                          end do
172                        end do
173                      end if
174
175                    end do
176                  end do
177                end do
178
179              end do
180            end do
181          end do
182
183c Contract over shell a
184
185          call dfill(MPQ*Lqmax3*NPB*NCA,0.0d0,t1,1)
186          do ica = 1,NCA
187            do mp = 1,MPP
188              ipa = PairP(1,mp)
189              ipb = PairP(2,mp)
190              ca = Acoefs(ipa,ica)
191              do nq = 1,Lqmax3
192                do mq = 1,MPQ
193                  t1(mq,nq,ipb,ica) = t1(mq,nq,ipb,ica) +
194     &                t2(mq,nq,mp)*ca
195                end do
196              end do
197            end do
198          end do
199
200c Contract over shell b
201
202          call dfill(MPQ*Lqmax3*NCAB,0.0d0,t2,1)
203          icab = 0
204          do ica = 1,NCA
205            do icb = 1,NCB
206              icab = icab+1
207              do ipb = 1,NPB
208                cb = Bcoefs(ipb,icb)
209                do nq = 1,Lqmax3
210                  do mq = 1,MPQ
211                    t2(mq,nq,icab) = t2(mq,nq,icab)
212     &                  +t1(mq,nq,ipb,ica)*cb
213                  end do
214                end do
215              end do
216            end do
217          end do
218
219
220c Loop over the components of the "C" and "D" shells.
221
222          do mc = 1,Lc2
223
224            call getNxyz(Lc,mc,Nxyz)
225            Ic = Nxyz(1)
226            Jc = Nxyz(2)
227            Kc = Nxyz(3)
228
229            do md = 1,Ld2
230
231              call getNxyz(Ld,md,Nxyz)
232              Id = Nxyz(1)
233              Jd = Nxyz(2)
234              Kd = Nxyz(3)
235
236c Define the product of the Hermite expansion coefficients for
237c overlap distribution "Q".
238
239              ie = 0
240              do Iq = 0,Ic+Id
241                do Jq = 0,Jc+Jd
242                  do Kq = 0,Kc+Kd
243                    do mq = 1,MPQ
244                      ie = ie+1
245                      E3(ie) = Eq(1,mq,0,Iq,Ic,Id)
246     &                        *Eq(2,mq,0,Jq,Jc,Jd)
247     &                        *Eq(3,mq,0,Kq,Kc,Kd)
248                    end do
249                  end do
250                end do
251              end do
252
253c Sum across (Iq,Jq,Kq).
254
255              icab = 0
256              ind_ca = ma
257              do ica = 1,NCA
258                ind_cb = mb
259                do icb = 1,NCB
260                  call dfill (MPQ,0.0d0,t3,1)
261                  icab = icab+1
262                  ie = 0
263                  do Iq = 0,Ic+Id
264                    do Jq = 0,Jc+Jd
265                      do Kq = 0,Kc+Kd
266
267                        nq = IJK(Iq,Jq,Kq)
268
269c Sum over Hermite functions to give integrals
270
271                        do mq = 1,MPQ
272                          ie = ie+1
273                          t3(mq) = t3(mq) + E3(ie)*t2(mq,nq,icab)
274                        end do
275
276                      end do
277                    end do
278                  end do
279
280c Contract over C shell
281
282                  do icc = 1,NCC
283                    call dfill (NPD,0.0d0,t4,1)
284                    do mq = 1,MPQ
285                      ipc = PairQ(1,mq)
286                      ipd = PairQ(2,mq)
287                      t4(ipd) = t4(ipd) + t3(mq)*Ccoefs(ipc,icc)
288                    end do
289
290c Contract over D shell
291
292                    do icd = 1,NCD
293                      do ipd = 1,NPD
294                        ERI(md,icd,mc,icc,ind_cb,ind_ca) =
295     &                      ERI(md,icd,mc,icc,ind_cb,ind_ca) +
296     &                      t4(ipd)*Dcoefs(ipd,icd)
297                      end do
298                    end do
299                  end do
300
301                  ind_cb = ind_cb+Lb2
302                end do ! icb
303                ind_ca = ind_ca+La2
304              end do ! ica
305
306            end do
307          end do
308
309        end do
310      end do
311
312      end
313