1c $Id$
2      Subroutine hf2(Axyz,Aprims,Acoefs,NPA,NCA,La,
3     &               Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,
4     &               Cxyz,Cprims,Ccoefs,NPC,NCC,Lc,
5     &               Dxyz,Dprims,Dcoefs,NPD,NCD,Ld,
6     &               bERI,Nints,canAB,canCD,canPQ,
7     &               DryRun,W0,maxW0)
8
9      Implicit none
10#include "stdio.fh"
11#include "mafdecls.fh"
12#include "errquit.fh"
13#include "case.fh"
14
15      Integer NPA,NCA,La
16      Integer NPB,NCB,Lb
17      Integer NPC,NCC,Lc
18      Integer NPD,NCD,Ld
19      Integer Nints,maxW0
20
21      Logical canAB,canCD,canPQ
22
23      Logical GenCon,DryRun
24
25c--> Cartesian Coordinates, Primitives & Contraction Coefficients
26
27      Double Precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
28      Double Precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
29      Double Precision Cxyz(3),Cprims(NPC),Ccoefs(NPC,NCC)
30      Double Precision Dxyz(3),Dprims(NPD),Dcoefs(NPD,NCD)
31c
32c--> Block of Electron Repulsion Integrals
33
34      Double Precision bERI(Nints)
35
36c--> Scratch Space
37
38      Double Precision W0(maxW0)
39
40c--> Local variables
41
42      Integer NCP,NCQ,MXD,nd
43      Integer La2,Lb2,Lc2,Ld2
44      Integer lp,lq,lr,lp3,lq3,lr3,lpq3
45      Integer i_ALPHAp,i_IPAIRp,i_ESp,i_Ep
46      Integer i_ALPHAq,i_IPAIRq,i_ESq,i_Eq
47      Integer i_left,i_right,i_pf
48      Integer MaxMem,MaxAvail,MemPerPair,MaxPairs,NPass
49      Integer maxpp,maxpq,ipp,ipq
50      Integer MPP,MPQ,MPR,NPP,NPQ,NPR
51      Integer lszp,lszq,ninti
52      Integer i_R0,i_IJK,i_P,i_Q,i_PQ,i_ff,i_Rj
53      Integer i_E3,i_t1,i_t2,i_t3,i_t4,i_eri,i_top
54      Integer i,j,istep
55      integer i_ffscr,l_ffscr,i_rscr,l_rscr,nvals
56
57      Logical doit_PQ
58
59c
60c Compute 4-ctr electron repulsion integrals (ERI) for four shells of
61c contracted Gaussian functions.
62c
63c******************************************************************************
64
65      MXD = 0
66      istep = MA_sizeof(MT_DBL,1,MT_INT)
67      if (istep .gt. 2) call errquit(
68     &    'hf2:Too many integers per real for Ipair addressing',911,
69     &       INT_ERR)
70      istep = 2/istep
71
72c Determine whether general or segmented contraction is used.
73
74      NCP = NCA*NCB
75      NCQ = NCC*NCD
76
77      GenCon = (NCP.ne.1) .or. (NCQ.ne.1)
78
79c Define the angular momentum of the overlap distributions.
80
81      Lp = La + Lb
82      Lq = Lc + Ld
83      Lr = Lp + Lq
84
85c Increment "Lr" to account for the order of differentiation.
86
87      Lr = Lr + MXD
88
89c Define the accumulated number of angular momentum functions <= Lr.
90
91      Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6
92      Lq3 = ((Lq+1)*(Lq+2)*(Lq+3))/6
93      Lpq3 = max(Lp3,Lq3)
94      Lr3 = ((Lr+1)*(Lr+2)*(Lr+3))/6
95
96c Define the prefactor of the overlap distribution "P".
97
98c Assign pointers to scratch space.
99
100      i_ALPHAp = 1
101      i_IPAIRp = i_ALPHAp + 2*(NPA*NPB)
102      i_left   = i_IPAIRp + 2*(NPA*NPB) - 1
103      i_ESp   = (maxW0+1) - 3*(NPA*NPB)
104      i_right = i_ESp
105
106      MaxMem = 1   ! take care of compiler warnings
107      if (DryRun) then
108        MaxMem = i_left + (maxW0 - (i_right-1))
109        NPP = NPA*NPB
110      else
111        if (i_left.ge.i_right) then
112          write(LuOut,*) 'HF2:  Insufficient scratch space.'
113          write(LuOut,*) '      needed    ',i_left + (maxW0-(i_right-1))
114          write(LuOut,*) '      allocated ',maxW0
115          write(LuOut,*) 'From the left '
116          write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
117          write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
118          write(LuOut,*) 'From the right '
119          write(LuOut,*) 'ESp   :  ',i_ESp
120          call errquit('hf2: Not enough scratch space for hfset AB',911,
121     &       MEM_ERR)
122        end if
123        call hfset(Axyz,Aprims,Acoefs,NPA,NCA,
124     &      Bxyz,Bprims,Bcoefs,NPB,NCB,
125     &      GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP)
126      end if
127
128c Define the prefactor of the overlap distribution "Q".
129
130c Assign pointers to scratch space.
131
132      i_ALPHAq = i_IPAIRp + 2*(NPA*NPB)
133      i_IPAIRq = i_ALPHAq + 2*(NPC*NPD)
134      i_left   = i_IPAIRq + 2*(NPC*NPD) - 1
135      i_ESq   = i_right - 3*(NPC*NPD)
136      i_right = i_ESq
137
138      if( DryRun )then
139        MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) )
140        NPQ = NPC*NPD
141      else
142        if (i_left.ge.i_right) then
143          write(LuOut,*) 'HF2:  Insufficient scratch space.'
144          write(LuOut,*) '      needed    ',i_left + (maxW0-(i_right-1))
145          write(LuOut,*) '      allocated ',maxW0
146          write(LuOut,*) 'From the left '
147          write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
148          write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
149          write(LuOut,*) 'ALPHAq:  ',i_ALPHAq
150          write(LuOut,*) 'IPAIRq:  ',i_IPAIRq
151          write(LuOut,*) 'From the right '
152          write(LuOut,*) 'ESp   :  ',i_ESp
153          write(LuOut,*) 'ESq   :  ',i_ESq
154          call errquit('hf2: Not enough scratch space for hfset CD',911,
155     &       MEM_ERR)
156        end if
157        call hfset(Cxyz,Cprims,Ccoefs,NPC,NCC,
158     &      Dxyz,Dprims,Dcoefs,NPD,NCD,
159     &      GenCon,W0(i_ALPHAq),W0(i_IPAIRq),W0(i_ESq),NPQ)
160      end if
161
162c Zero out integral block. Return if NPR = 0, i.e. no integrals
163
164      NPR = NPP*NPQ
165      if (NPR .eq. 0) then
166!     if (.not.DryRun) call dfill(Nints,0.0d00,beri,1)
167      if (.not.DryRun) call dcopy(Nints,0d0,0,beri,1)
168         go to 99
169      endif
170
171c Define cartesian components for shells, bra and ket length and
172c number of integrals.
173
174      La2 = (la+1)*(la+2)/2
175      Lb2 = (lb+1)*(lb+2)/2
176      lszp = La2*Lb2*NCP
177
178      Lc2 = (lc+1)*(lc+2)/2
179      Ld2 = (ld+1)*(ld+2)/2
180      lszq = Lc2*Ld2*NCQ
181
182      ninti = lszp*lszq
183      if (.not.DryRun) call dfill(Ninti,0.0d00,beri,1)
184
185c Define the Hermite linear expansion coefficients.
186
187c Assign pointers to scratch space.
188
189      i_Ep   = i_IPAIRq + 2*(NPC*NPD)
190      i_Eq   = i_Ep     + 3*NPP*(MXD+1)*((Lp+1)*(La+1)*(Lb+1))
191      i_pf   = i_Eq     + 3*NPQ*(MXD+1)*((Lq+1)*(Lc+1)*(Ld+1))
192      i_left = i_pf     + 2*max(NPP,NPQ) - 1
193
194      if( DryRun )then
195        MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) )
196      else
197        if (i_left.ge.i_right) then
198          write(LuOut,*) 'HF2:  Insufficient scratch space.'
199          write(LuOut,*) '      needed    ',i_left + (maxW0-(i_right-1))
200          write(LuOut,*) '      allocated ',maxW0
201          write(LuOut,*) 'From the left '
202          write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
203          write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
204          write(LuOut,*) 'ALPHAq:  ',i_ALPHAq
205          write(LuOut,*) 'IPAIRq:  ',i_IPAIRq
206          write(LuOut,*) 'Ep    :  ',i_Ep
207          write(LuOut,*) 'Eq    :  ',i_Eq
208          write(LuOut,*) 'pf    :  ',i_pf
209          write(LuOut,*) 'From the right '
210          write(LuOut,*) 'ESp   :  ',i_ESp
211          write(LuOut,*) 'ESq   :  ',i_ESq
212          call errquit('hf2: Not enough scratch space for hfmke',911,
213     &       MEM_ERR)
214        end if
215
216        do nd = 0,MXD
217          call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep),W0(i_pf),
218     &        nd,NPP,MXD,La,Lb)
219          call hfmke(Cxyz,Dxyz,W0(i_ALPHAq),W0(i_ESq),W0(i_Eq),W0(i_pf),
220     &        nd,NPQ,MXD,Lc,Ld)
221        end do
222
223      end if
224
225c Evaluate the auxiliary function integrals.
226
227c   Determine memory requirements and assign pointers to scratch space.
228c   If insufficient memory, do multipassing.
229
230      i_IJK = i_Eq  + 3*NPQ*(MXD+1)*((Lq+1)*(Lc+1)*(Ld+1))
231      i_R0  = i_IJK + (Lr+1)**3
232
233      MaxAvail = maxW0-i_R0-3*(NPP+NPQ)+1
234      MemPerPair = (Lr+2)*Lr3 + 5
235      if (MaxAvail-MemPerPair .le. 0) then
236        write(LuOut,*) 'HF2:  Insufficient scratch space.'
237        write(LuOut,*) '      needed    ',i_R0+3*(NPP+NPQ)
238     &      +MemPerPair*NPR
239        write(LuOut,*) '      allocated ',maxW0
240        write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
241        write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
242        write(LuOut,*) 'ALPHAq:  ',i_ALPHAq
243        write(LuOut,*) 'IPAIRq:  ',i_IPAIRq
244        write(LuOut,*) 'Ep    :  ',i_Ep
245        write(LuOut,*) 'Eq    :  ',i_Eq
246        write(LuOut,*) 'R0    :  ',i_R0
247        write(LuOut,*) 'IJK   :  ',i_IJK
248        call errquit('hf2: Not enough scratch space for hf2mkr',911,
249     &       MEM_ERR)
250      end if
251
252      MaxPairs = Maxavail/MemPerPair
253
254      if (MaxPairs .ge. NPR) then
255        maxpp = NPP
256        maxpq = NPQ
257      else if (NPP .ge. NPQ) then
258        if (MaxPairs .ge. NPP) then
259          maxpp = NPP
260          maxpq = MaxPairs/NPP
261        else if (MaxPairs .ge. NPQ) then
262          maxpq = NPQ
263          maxpp = MaxPairs/NPQ
264        else
265          NPass = (NPP-1)/MaxPairs+1
266          maxpp = min((NPP-1)/NPass+1,MaxPairs)
267          maxpq = 1
268        end if
269      else
270        if (MaxPairs .ge. NPQ) then
271          maxpq = NPQ
272          maxpp = MaxPairs/NPQ
273        else if (MaxPairs .ge. NPP) then
274          maxpp = NPP
275          maxpq = MaxPairs/NPP
276        else
277          NPass = (NPP-1)/MaxPairs+1
278          maxpp = min((NPP-1)/NPass+1,MaxPairs)
279          maxpq = 1
280        end if
281      end if
282
283      do ipp = 0,NPP-1,maxpp
284        MPP = min(maxpp,NPP-ipp)
285        do ipq = 0,NPQ-1,maxpq
286          MPQ = min(maxpq,NPQ-ipq)
287          MPR = MPP*MPQ
288          i_P   = i_R0  + MPR*Lr3
289          i_Q   = i_P   + 3*MPP
290          i_PQ  = i_Q   + 3*MPQ
291          i_ff  = i_PQ  + MPR*3
292          i_Rj  = i_ff  + 2*MPR
293          i_top = i_Rj  + MPR*(Lr+1)*Lr3 - 1
294
295          if (Lp.eq.Lq) then
296            doit_PQ = MPP.ge.MPQ
297          else if (abs((MPP-MPQ)).ge.4) then
298            doit_PQ = MPP.ge.MPQ
299          else if (Lp.lt.Lq) then
300            doit_PQ = .true.
301          else
302            doit_PQ = .false.
303          endif
304
305          if( DryRun )then
306            MaxMem = max( MaxMem, i_top )
307          else
308             if(doscreen) then
309                nvals=2*MPP*MPQ
310                if (.not.MA_Push_get(mt_dbl,nvals,'ffscr',
311     C               l_ffscr,i_ffscr)) call errquit(
312     E               'hf2: ran out of MA memory',0,MEM_ERR)
313                nvals=2*MPP*MPQ*lr3*(lr+1)
314                if (.not.MA_Push_get(mt_dbl,nvals,'rscr',
315     C               l_rscr,i_rscr)) call errquit(
316     E               'hf2: ran out of MA memory',1,MEM_ERR)
317                else
318                   nvals=1
319                   i_rscr=0
320                   i_ffscr=0
321                endif
322            if (doit_PQ) then
323              call hf2mkr(Axyz,Bxyz,Cxyz,Dxyz,
324     &            W0(i_ALPHAp+2*ipp),W0(i_ALPHAq+2*ipq),
325     &            W0(i_R0),W0(i_IJK),W0(i_P),W0(i_Q),W0(i_PQ),
326     &            W0(i_ff),W0(i_Rj),
327     &            dbl_mb(i_ffscr),dbl_mb(i_rscr),
328     &            MPP,MPQ,Lr,Lr3)
329            else
330              call hf2mkr(Cxyz,Dxyz,Axyz,Bxyz,
331     &            W0(i_ALPHAq+2*ipq),W0(i_ALPHAp+2*ipp),
332     &            W0(i_R0),W0(i_IJK),W0(i_Q),W0(i_P),W0(i_PQ),
333     &            W0(i_ff),W0(i_Rj),
334     &            dbl_mb(i_ffscr),dbl_mb(i_rscr),
335     &            MPQ,MPP,Lr,Lr3)
336            end if
337            if(doscreen) then
338               if (.not.ma_chop_stack(l_ffscr)) call errquit(
339     E              'hf2: machopstack failed',0,MEM_ERR)
340            endif
341          end if
342
343c Compute the ERI.
344
345c Assign pointers to scratch space.
346
347          i_eri = i_P
348          i_t2 = i_eri + ninti
349          if (gencon) then
350            if (doit_PQ) then
351              i_E3 = i_t2 + MPQ*Lq3*max(MPP,NCP)
352              i_t1 = i_E3
353              i_t3 = i_E3 + MPQ*Lq3
354              i_t4 = i_t3 + MPQ
355              i_top = max(i_t1+MPQ*Lq3*NPB*NCA,i_E3+MPP,i_t4+NPD)-1
356            else
357              i_E3 = i_t2 + MPP*Lp3*max(MPQ,NCQ)
358              i_t1 = i_E3
359              i_t3 = i_E3 + MPP*Lp3
360              i_t4 = i_t3 + MPP
361              i_top = max(i_t1+MPP*Lp3*NPD*NCC,i_E3+MPQ,i_t4+NPB)-1
362            end if
363          else
364            i_E3 = i_t2 + MPR*Lpq3
365            i_top  = i_E3 + max(MPP,MPQ)-1
366          end if
367          if( DryRun )then
368            MaxMem = max( MaxMem, i_top )
369          else if (i_top.gt.maxW0) then
370            write(LuOut,*) 'HF2:  Insufficient scratch space.'
371            write(LuOut,*) '      needed    ',i_top
372            write(LuOut,*) '      allocated ',maxW0
373            write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
374            write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
375            write(LuOut,*) 'ALPHAq:  ',i_ALPHAq
376            write(LuOut,*) 'IPAIRq:  ',i_IPAIRq
377            write(LuOut,*) 'Ep    :  ',i_Ep
378            write(LuOut,*) 'Eq    :  ',i_Eq
379            write(LuOut,*) 'IJK   :  ',i_IJK
380            write(LuOut,*) 'R0    :  ',i_R0
381            write(LuOut,*) 'eri   :  ',i_eri
382            write(LuOut,*) 't2    :  ',i_t2
383            write(LuOut,*) 'E3    :  ',i_E3
384            if (gencon) then
385              write(LuOut,*) 't1  :  ',i_t1
386              write(LuOut,*) 't3  :  ',i_t3
387              write(LuOut,*) 't4  :  ',i_t4
388            end if
389            call errquit('hf2: Not enough scratch space for hferi',911,
390     &       MEM_ERR)
391          else if (doit_PQ) then
392            if (GenCon) then
393              call hferi_gen(W0(i_Ep+3*ipp),W0(i_Eq+3*ipq),
394     &            W0(i_IPAIRp+istep*ipp),W0(i_IPAIRq+istep*ipq),
395     &            W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),
396     &            W0(i_t1),W0(i_t2),W0(i_t3),W0(i_t4),
397     &            MPP,MPQ,NPP,NPQ,La,Lb,Lc,Ld,La2,Lb2,Lc2,Ld2,
398     &            Lq,Lq3,Lr,Acoefs,Bcoefs,Ccoefs,Dcoefs,
399     &            NPA,NPB,NPC,NPD,NCA,NCB,NCC,NCD,
400     &            MXD,canAB,canCD,canPQ)
401            else
402              call hferi(W0(i_Ep+3*ipp),W0(i_Eq+3*ipq),
403     &            W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),W0(i_t2),
404     &            MPP,MPQ,NPP,NPQ,Nints,La,Lb,Lc,Ld,Lr,MXD,
405     &            canAB,canCD,canPQ)
406            end if
407            call daxpy(ninti,1.0D0,W0(i_eri),1,beri,1)
408          else
409            if (GenCon) then
410              call hferi_gen(W0(i_Eq+3*ipq),W0(i_Ep+3*ipp),
411     &            W0(i_IPAIRq+istep*ipq),W0(i_IPAIRp+istep*ipp),
412     &            W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),
413     &            W0(i_t1),W0(i_t2),W0(i_t3),W0(i_t4),
414     &            MPQ,MPP,NPQ,NPP,Lc,Ld,La,Lb,Lc2,Ld2,La2,Lb2,
415     &            Lp,Lp3,Lr,Ccoefs,Dcoefs,Acoefs,Bcoefs,
416     &            NPC,NPD,NPA,NPB,NCC,NCD,NCA,NCB,
417     &            MXD,canCD,canAB,canPQ)
418            else
419              call hferi(W0(i_Eq+3*ipq),W0(i_Ep+3*ipp),
420     &            W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),W0(i_t2),
421     &            MPQ,MPP,NPQ,NPP,Ninti,Lc,Ld,La,Lb,Lr,MXD,
422     &            canCD,canAB,canPQ)
423            end if
424            j = 1
425            do i = 0,lszp-1
426              call daxpy(lszq,1.0d0,W0(i_eri+i),lszp,beri(j),1)
427              j = j+lszq
428            end do
429          end if
430
431        end do
432      end do
433
434c Return the maximum amount of scratch space required by a "dry run".
435
436   99 if( DryRun ) maxW0 = MaxMem
437c
438      end
439
440