1      Subroutine hferi(Ep,Eq,R0,IJK,ERI,E3,sum,MPP,MPQ,
2     &                 NPP,NPQ,Nint,La,Lb,Lc,Ld,Lr,MXD,
3     &                 canAB,canCD,canPQ)
4c $Id$
5
6      Implicit none
7#include "sh_order.fh"
8      integer mpp,mpq,npp,npq,nint
9      integer la,lb,lc,ld,lr,mxd
10      integer ld2,lqmax,lqmax3
11      integer iq,jq,nn,ia,ja,ka,mb_limit
12      integer ij,mb,ib,jb,kb,ip,jp,kp,mp,ir,jr,lq,nq
13      integer kr,nr,mq,mc,ic,jc,kc,md_limit,kl,md,id,jd,kd
14      integer lb2,lc2,kq,ma,nc,la6,lb6,lc6,ld6
15      integer nq1,iq1,iqfin,iqlast
16      Logical canAB,canCD,canPQ
17
18c--> Hermite Linear Expansion Coefficients
19
20      Double precision Ep(3,NPP,0:MXD,0:(La+Lb),0:La,0:Lb)
21      Double precision Eq(3,NPQ,0:MXD,0:(Lc+Ld),0:Lc,0:Ld)
22
23c--> Auxiliary Function Integrals & Index
24
25      Double precision R0(MPQ,MPP,*)
26      Integer IJK(0:Lr,0:Lr,0:Lr)
27
28c--> ERI
29
30      Double precision ERI(Nint)
31
32c--> Scratch Space
33
34      Double precision E3(*),sum(MPQ,*),erinn,zot,
35     ,     zot1,zot2
36      Integer Nxyz(3),nqmx
37c
38c Compute electron repulsion integrals (ERI).
39c
40c     Formula:
41c
42c               __
43c               \     Ic,Id;n10   Jc,Jd;n11   Kc,Kd;n12
44c     ERI  =    /   Ex          Ey          Ez          SUM
45c               --    Iq          Jq          Kq           Iq,Jq,Kq
46c            Iq,Jq,Kq
47c
48c                            __
49c                           \     Lq   Ia,Ib;n7   Ja,Jb;n8   Ka,Kb;n9
50c         SUM          =    /  (-1)  Ex         Ey         Ez         R
51c            Iq,Jq,Kq       --         Ip         Jp         Kp        Ir,Jr,Kr
52c                        Ip,Jp,Kp
53c
54c                                Ir = (Ip+n1) + (Iq+n4)
55c                        where   Jr = (Jp+n2) + (Jq+n5)
56c                                Kr = (Kp+n3) + (Kq+n6)
57c
58c                         and    Lq = (Iq+n4) + (Jq+n5) + (Kq+n6)
59c
60c N.B.  For simple ERI (i.e., no derivative integrals)  n[1-12] = 0!
61c
62c******************************************************************************
63*rak:      integer num_pg, num_tot, num_qg
64*rak:      save num_pg, num_tot, num_qg
65*rak:      data num_pg /0/
66*rak:      data num_qg /0/
67*rak:      data num_tot /0/
68*rak:
69*rak:      if(MPP.ge.MPQ) then
70*rak:        num_pg = num_pg + 1
71*rak:      else
72*rak:        num_qg = num_qg + 1
73*rak:      endif
74*rak:      num_tot = num_tot + 1
75*rak:      if (num_tot.lt.3200.or.mod(num_tot,1000).eq.0) then
76*rak:        write(6,*)' num_pg  ',num_pg
77*rak:        write(6,*)' num_qg  ',num_qg
78*rak:        write(6,*)' num_tot ',num_tot
79*rak:      endif
80c
81c General case:  [ab|cd]
82
83c Define the number of shell components on each center.
84
85      Lb2 = ((Lb+1)*(Lb+2))/2
86      Lc2 = ((Lc+1)*(Lc+2))/2
87      Ld2 = ((Ld+1)*(Ld+2))/2
88      la6=(la*(la+1)*(la+2))/6
89      lb6=(lb*(lb+1)*(lb+2))/6
90      lc6=(lc*(lc+1)*(lc+2))/6
91      ld6=(ld*(ld+1)*(ld+2))/6
92      md_limit = Ld2
93
94c Initialize the block of ERIs.
95
96      Lqmax  = Lc + Ld
97      Lqmax3 = ((Lqmax+1)*(Lqmax+2)*(Lqmax+3))/6
98
99c Loop over the components of the "A" and "B" shells.
100      nqmx=0
101      do Iq = 0,Lqmax
102         do Jq = 0,Lqmax-Iq
103!DEC$ LOOP COUNT MAX=10, MIN=1
104            do Kq = 0,Lqmax-Iq-Jq
105               nqmx = max(IJK(Iq,Jq,Kq),nqmx)
106            enddo
107         enddo
108      enddo
109      nn = 0
110      do ma = 1,((La+1)*(La+2))/2
111
112       nc =  la6 + ma
113       ia = Ixyz(1,nc)
114       ja = Ixyz(2,nc)
115       ka = Ixyz(3,nc)
116
117
118        if( canAB )then
119           mb_limit = ma
120           ij = (ma*(ma-1))/2
121        else
122           mb_limit = Lb2
123           ij = (ma-1)*Lb2
124        end if
125
126        do mb = 1,mb_limit
127           ij=ij+1
128
129
130       nc =  lb6 + mb
131       ib = Ixyz(1,nc)
132       jb = Ixyz(2,nc)
133       kb = Ixyz(3,nc)
134
135c Sum across (Ip,Jp,Kp) for each value of (Iq,Jq,Kq).
136          call dcopy(mpq*nqmx,0d0,0,sum,1)
137
138          do Ip = 0,Ia+Ib
139             do Jp = 0,Ja+Jb
140                do Kp = 0,Ka+Kb
141                   if(MPP.eq.1) then
142                      E3(1) = Ep(1,1,0,Ip,Ia,Ib)*
143     &                     Ep(2,1,0,Jp,Ja,Jb)*
144     &                     Ep(3,1,0,Kp,Ka,Kb)
145                   else
146
147c Define the product of the Hermite expansions coefficients for
148c overlap distribution "P".
149!DEC$ LOOP COUNT MAX=30, MIN=1
150                      do mp = 1,MPP
151                         E3(mp) = Ep(1,mp,0,Ip,Ia,Ib)*
152     &                        Ep(2,mp,0,Jp,Ja,Jb)*
153     &                        Ep(3,mp,0,Kp,Ka,Kb)
154                      end do
155                   endif
156
157                do Iq = 0,Lqmax
158                   Ir = Ip + Iq
159                  do Jq = 0,Lqmax-Iq
160                     Jr = Jp + Jq
161                      Lq = Iq + Jq -1
162                    do Kq = 0,Lqmax-Iq-Jq
163
164                      nq = IJK(Iq,Jq,Kq)
165
166                      Kr = Kp + Kq
167
168                      nr = IJK(Ir,Jr,Kr)
169
170c Include the factor of (-1)**(Iq+Jq+Kq).
171                      Lq=Lq+1
172#if defined(GCC4) || defined(PGLINUX)
173                      if(IAND(Lq,1).eq.1)then
174#else
175                      if(AND(Lq,1).eq.1)then
176#endif
177                        if(MPQ.eq.1) then
178!DEC$ LOOP COUNT MAX=30, MIN=1
179                           do mp = 1,MPP
180                              sum(1,nq) = sum(1,nq)-E3(mp)*R0(1,mp,nr)
181                           end do
182                        else
183                        do mp = 1,MPP
184                           zot=-E3(mp)
185!DEC$ LOOP COUNT MAX=30, MIN=1
186                          do mq = 1,MPQ
187                            sum(mq,nq) = sum(mq,nq)+zot*R0(mq,mp,nr)
188                          end do
189                        end do
190                        endif
191                      else
192                        if(MPQ.eq.1) then
193!DEC$ LOOP COUNT MAX=30, MIN=1
194                           do mp = 1,MPP
195                              sum(1,nq) = sum(1,nq)+E3(mp)*R0(1,mp,nr)
196                           end do
197                        else
198                        do mp = 1,MPP
199                           zot=E3(mp)
200!DEC$ LOOP COUNT MAX=30, MIN=1
201                          do mq = 1,MPQ
202                            sum(mq,nq) = sum(mq,nq)+zot*R0(mq,mp,nr)
203                          end do
204                        end do
205                        endif
206                      end if
207
208                    end do
209                  end do
210                end do
211
212              end do
213            end do
214          end do
215
216c Loop over the components of the "C" and "D" shells.
217
218          do mc = 1,Lc2
219
220             nc =  lc6 + mc
221             ic = Ixyz(1,nc)
222             jc = Ixyz(2,nc)
223             kc = Ixyz(3,nc)
224
225            if( canCD ) md_limit = mc
226
227
228            if( canAB )then
229               kl = (mc*(mc-1))/2
230            else
231               kl = (mc-1)*Ld2
232            end if
233            do md = 1,md_limit
234
235              if( canPQ )then
236                 kl=kl+1
237                 if( kl.gt.ij ) go to 480
238              end if
239
240             nc =  ld6 + md
241             id = Ixyz(1,nc)
242             jd = Ixyz(2,nc)
243             kd = Ixyz(3,nc)
244
245              nn = nn + 1
246
247c Sum across (Iq,Jq,Kq).
248c Define the product of the Hermite expansion coefficients for
249c overlap distribution "Q" and calculate eri
250              erinn=0d0
251
252              if(MPQ.eq.1) then
253                 do Iq = 0,Ic+Id
254                    do Jq = 0,Jc+Jd
255!DEC$ LOOP COUNT MAX=30, MIN=1
256                       do Kq = 0,Kc+Kd
257                          nq = IJK(Iq,Jq,Kq)
258                          ERInn = ERInn +
259     &                         Eq(1,1,0,Iq,Ic,Id)*
260     &                         Eq(2,1,0,Jq,Jc,Jd)*
261     &                         Eq(3,1,0,Kq,Kc,Kd)*sum(1,nq)
262                       enddo
263                    enddo
264                 enddo
265              else
266                 iqfin=0
267                 iqlast=ic+id
268                 if(iqlast.gt.0) then
269                    iqfin=-1
270#if defined(GCC4) || defined(PGLINUX)
271                    if(iand(iqlast,1).eq.0) then
272#else
273                    if(and(iqlast,1).eq.0) then
274#endif
275                       iqfin=ic+id
276                       iqlast=iqlast-2
277                    endif
278                    do Iq = 0,iqlast,2
279                       iq1=iq+1
280                       do Jq = 0,Jc+Jd
281                          do Kq = 0,Kc+Kd
282                             nq = IJK(Iq,Jq,Kq)
283                             nq1 = IJK(Iq1,Jq,Kq)
284!DEC$ LOOP COUNT MAX=30, MIN=1
285                             do mq = 1,MPQ
286                                ERInn = ERInn +
287     &                               (Eq(1,mq,0,Iq,Ic,Id)*sum(mq,nq)+
288     &                               Eq(1,mq,0,Iq1,Ic,Id)*sum(mq,nq1))*
289     &                               Eq(2,mq,0,Jq,Jc,Jd)*
290     &                               Eq(3,mq,0,Kq,Kc,Kd)
291                             end do
292                          end do
293                       enddo
294                    enddo
295                 endif
296                 if(iqfin.ne.-1) then
297                    do Jq = 0,Jc+Jd
298                       do Kq = 0,Kc+Kd
299                          nq = IJK(Iqfin,Jq,Kq)
300!DEC$ LOOP COUNT MAX=30, MIN=1
301                          do mq = 1,MPQ
302                             ERInn = ERInn +
303     &                            Eq(1,mq,0,Iqfin,Ic,Id)*
304     &                            Eq(2,mq,0,Jq,Jc,Jd)*
305     &                            Eq(3,mq,0,Kq,Kc,Kd)*sum(mq,nq)
306                          end do
307                       end do
308                    enddo
309                 endif
310              endif
311              eri(nn)=erinn
312
313            end do              ! md
314  480       continue
315          end do                ! mc
316        end do                  ! mb
317      end do                    ! ma
318
319      end
320