1      Subroutine hf1(Axyz,Aprims,Acoefs,NPA,NCA,La,
2     &               Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,
3     &               Cxyz,zan,exinv,ncenters,
4     &               bO2I,bKEI,bNAI,Nints,O2I,KEI,NAI,canAB,
5     &               DryRun,W0,maxW0)
6c $Id$
7      Implicit none
8      integer NPA,NCA,La,NPB,NCB,Lb,ncenters,Nints,maxW0
9      logical O2I,KEI,NAI,canAB
10      logical GenCon,DryRun
11
12c--> Cartesian Coordinates, Primitives & Contraction Coefficients
13
14      double precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
15      double precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
16
17c--> Nuclear Cartesian Coordinates, Charges and Inverse Exponents
18
19      double precision Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
20
21c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals
22
23      double precision bO2I(Nints),bKEI(Nints),bNAI(Nints)
24
25c--> Scratch Space.
26
27      double precision W0(maxW0)
28
29c--> Local variables
30
31      integer MXD,NCP,NPP,La2,Lb2,Li,Lp,Lp3,MaxMem,nintlocal,lprod,nd
32      integer i_ALPHAp,i_IPAIRp,i_ESp,i_left,i_right,i_top,i_Ep,i_ptr,
33     &    i_pf,i_prim_ints,i_half_ints,i_Ti,i_R0,i_IJK,i_P,i_RS,i_PC,
34     &    i_ff,i_Rj
35#if defined(INTDEBUG)
36      integer jjjj,iiii
37#endif
38#include "mafdecls.fh"
39#include "stdio.fh"
40#include "errquit.fh"
41c
42c Compute the overlap, kinetic energy, and nuclear attraction integrals for
43c two shells of contracted Gaussians functions. This driver is NOT capable of
44c evaluating integral derivatives.
45c
46c******************************************************************************
47#if defined(INTDEBUG)
48      if (.not.dryrun) then
49        write(LuOut,*)' inside hf1 '
50        write(LuOut,*)' npa,nca,la = ',npa,nca,la
51        write(LuOut,*)' npb,ncb,lb = ',npb,ncb,lb
52        write(LuOut,*)' ncenters   = ',ncenters
53        write(LuOut,*)' NINTS       = ',nints
54        write(LuOut,*)' maxW0      = ',maxw0
55        write(LuOut,*)' <canAB:DryRun>-<',canab,':',dryrun,'>'
56        write(LuOut,*)' <o2i:kei:nai>-<',o2i,':',kei,':',nai,'>'
57        write(6,'(a,3(2x,1pd20.10))')' Axyz =',Axyz
58        write(6,'(a,3(2x,1pd20.10))')' Bxyz =',Bxyz
59        write(6,'(a,100(3(2x,1pd20.10/)))')' Cxyz =',Cxyz
60        do jjjj = 1,nca
61        do iiii = 1,npa
62          write(6,'(a,i3,a,2(2x,1pd20.10))')
63     &        'Aprims:Acoeffs:(',iiii,') =',Aprims(iiii),
64     &        Acoefs(iiii,jjjj)
65        enddo
66        enddo
67        do jjjj = 1,ncb
68        do iiii = 1,npb
69          write(6,'(a,i3,a,2(2x,1pd20.10))')
70     &        'Bprims:Bcoeffs:(',iiii,') =',Bprims(iiii),
71     &        Bcoefs(iiii,jjjj)
72        enddo
73        enddo
74      endif
75      if (.not.dryrun) then
76        write(LuOut,*)' Li/Lj',La,'/',Lb
77        write(LuOut,*)' i_ngen ',nca
78        write(LuOut,*)' j_ngen ',ncb
79        write(LuOut,*)' int_hf1: lstv',Nints
80        write(LuOut,*)' int_hf1: lscr',maxW0
81        call dcopy(maxW0,0d0,0,W0,1)
82        if (O2I) call dcopy(Nints,0d0,0,bO2I,1)
83        if (KEI) call dcopy(Nints,0d0,0,bKEI,1)
84        if (NAI) call dcopy(Nints,0d0,0,bNAI,1)
85*debug_ma:        call MA_summarize_allocated_blocks()
86*debug_ma:        write(LuOut,*)' int_hf1: ma verify 1-b4'
87*debug_ma:        status = ma_verify_allocator_stuff()
88*debug_ma:        write(LuOut,*)' verstat = ',status
89*debug_ma:        write(LuOut,*)' int_hf1: ma verify 1-af'
90        call util_flush(6)
91      endif
92      if (.not.dryrun) then
93        call hf_print_set(1)
94        call hf_print('hf1: a shell',axyz,aprims,acoefs,npa,nca,la)
95        call hf_print('hf1: b shell',bxyz,bprims,bcoefs,npb,ncb,lb)
96        call hf_print_set(0)
97      endif
98#endif
99
100      MXD = 0
101
102c Determine whether general or segmented contraction is used.
103
104      NCP = NCA*NCB
105      GenCon = NCP.ne.1
106
107c To determine all the Hermite expansion coefficients required to evaluate
108c the kinetic energy integrals, increment the angular momenta by one.
109
110      if( KEI )then
111        Li = 1
112      else
113        Li = 0
114      end if
115
116c Define the angular momentum of the overlap distribution.
117
118      Lp = La + Lb
119
120c Increment "Lp" to account for the order of differentiation.
121
122      Lp = Lp + MXD
123
124c Define the accumulated number of angular momentum functions <= Lp.
125
126      Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6
127
128c Define the prefactor of the overlap distribution "P".
129
130c Assign pointers to scratch space.
131
132      i_ALPHAp = 1
133      i_IPAIRp = i_ALPHAp + 2*(NPA*NPB)
134      i_left   = i_IPAIRp + 2*(NPA*NPB) - 1
135
136      i_ESp   = (maxW0+1) - 3*(NPA*NPB)
137      i_right = i_ESp
138
139      MaxMem = 1    ! take care of compiler warnings
140      if( DryRun )then
141
142        MaxMem = i_left + (maxW0 - (i_right-1))
143        NPP = NPA*NPB
144
145      else if (i_left.ge.i_right) then
146
147        write(LuOut,*) 'HF1:  Insufficient scratch space.'
148        write(LuOut,*) '       needed    ',i_left+(maxW0-(i_right-1))
149        write(LuOut,*) '       allocated ',maxW0
150        write(LuOut,*) ' DryRun ',DryRun
151
152        write(LuOut,*) 'From the left '
153        write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
154        write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
155        write(LuOut,*) 'From the right '
156        write(LuOut,*) 'ESp   :  ',i_ESp
157
158        call errquit('hf1: insufficient memory for hfset',911, MEM_ERR)
159
160      else
161
162        call hfset(Axyz,Aprims,Acoefs,NPA,NCA,
163     &      Bxyz,Bprims,Bcoefs,NPB,NCB,
164     &      GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP)
165
166      end if
167
168      La2 = (La+1)*(La+2)/2
169      Lb2 = (Lb+1)*(Lb+2)/2
170
171c Zero out the integrals. Return if screening test gives no pairs
172
173      if (.not.DryRun) then
174        nintlocal = La2*Lb2*nca*ncb
175        if (O2I) call dcopy(nintlocal,0d0,0,bO2I,1)
176        if (KEI) call dcopy(nintlocal,0d0,0,bKEI,1)
177        if (NAI) call dcopy(nintlocal,0d0,0,bNAI,1)
178        if (NPP.eq.0) return
179      end if
180
181c Define the Hermite linear expansion coefficients.
182
183c Assign pointers to scratch space.
184
185      lprod = ((La+Li)+(Lb+Li)+1)*((La+Li)+1)*((Lb+Li)+1)
186
187C      write (LuOut,*) 'before hfmke'
188
189      i_Ep   = i_IPAIRp + 2*(NPA*NPB)
190      i_ptr  = i_Ep     + 3*NPP*(MXD+1)*lprod
191
192      i_prim_ints = i_ptr            ! take care of compiler warnings
193      i_half_ints = i_prim_ints + NPP
194
195      if (gencon) then
196        i_prim_ints = i_ptr
197        i_half_ints = i_prim_ints + NPP
198        i_ptr       = i_half_ints + NPA*NCB
199      end if
200
201      i_pf   = i_ptr
202      i_left = i_pf     + 2*NPP - 1
203
204      if( DryRun )then
205
206        MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) )
207
208      else if( i_left.ge.i_right) then
209
210        write(LuOut,*) 'HF1:  Insufficient scratch space.'
211        write(LuOut,*) '       needed    ',i_left+(maxW0-(i_right-1))
212        write(LuOut,*) '       allocated ',maxW0
213        write(LuOut,*) ' DryRun ',DryRun
214
215        write(LuOut,*) 'From the left '
216        write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
217        write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
218        write(LuOut,*) 'Ep    :  ',i_Ep
219        write(LuOut,*) 'pf    :  ',i_pf
220        write(LuOut,*) 'From the right '
221        write(LuOut,*) 'ESp   :  ',i_ESp
222
223        call errquit('hf1: insufficient memory for hfmke',911, MEM_ERR)
224
225      else
226
227        do nd = 0,MXD
228          call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep),W0(i_pf),
229     &        nd,NPP,MXD,La+Li,Lb+Li)
230        end do
231
232      end if
233
234c Compute the 2-center overlap integrals, <a|S|b>.
235
236      if( O2I )then
237        if( .not. DryRun )then
238          if (gencon) then
239            call hf2oi_gc(W0(i_Ep),bO2I,W0(i_prim_ints),W0(i_half_ints),
240     &          Acoefs,Bcoefs,W0(i_IPAIRp),NPA,NPB,NCA,NCB,NPP,
241     &          La,Lb,La2,Lb2,Li,canAB)
242          else
243            call hf2oi(W0(i_Ep),bO2I,Nints,NPP,La,Lb,Li,canAB)
244          endif
245        end if
246      end if
247
248c Compute kinetic energy integrals, <a|T|b>.
249
250      if (KEI) then
251
252c Assign pointers to scratch space.
253
254        i_Ti  = i_ptr
255        i_top = i_Ti + NPP - 1
256
257        if( DryRun )then
258
259          MaxMem = max( MaxMem, i_top )
260
261        else if( i_top.gt.maxW0 )then
262
263          write(LuOut,*) 'HF1:  Insufficient scratch space.'
264          write(LuOut,*) '       needed    ',i_top
265          write(LuOut,*) '       allocated ',maxW0
266          write(LuOut,*) ' DryRun ',DryRun
267
268          write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
269          write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
270          write(LuOut,*) 'Ep    :  ',i_Ep
271          write(LuOut,*) 'Ti    :  ',i_Ti
272
273          call errquit('hf1: insufficient memory for hfkei',911,
274     &       MEM_ERR)
275
276        else if (gencon) then
277
278          call hfkei_gc(W0(i_ALPHAp),W0(i_Ep),bKEI,
279     &        W0(i_prim_ints),W0(i_half_ints),W0(i_Ti),
280     &        Acoefs,Bcoefs,W0(i_IPAIRp),
281     &        NPA,NPB,NCA,NCB,NPP,La,Lb,La2,Lb2,Li,canAB)
282        else
283
284          call hfkei(W0(i_ALPHAp),W0(i_Ep),bKEI,W0(i_Ti),
285     &        Nints,NPP,La,Lb,Li,canAB)
286
287        end if
288
289      end if
290
291c Compute nuclear attraction integrals, <a|V|b>.
292
293      if( NAI )then
294
295c Define the auxiliary function integrals.
296
297c Assign scratch space.
298
299        i_R0  = i_ptr
300        i_IJK = i_R0  + NPP*Lp3
301        i_P   = i_IJK + (Lp+1)**3
302        i_RS  = i_P   + NPP*3
303        i_PC  = i_RS  + NPP
304        i_ff  = i_PC  + NPP*3
305        i_Rj  = i_ff  + NPP*2
306        i_top = i_Rj  + NPP*(Lp+1)*Lp3 - 1
307
308        if( DryRun )then
309
310          MaxMem = max( MaxMem, i_top )
311C          write (LuOut,*) MaxMem,maxW0
312
313        else if( i_top.gt.maxW0 .and. .not.DryRun )then
314
315          write(LuOut,*) 'HF1:  Insufficient scratch space.'
316          write(LuOut,*) '       needed    ',i_top
317          write(LuOut,*) '       allocated ',maxW0
318          write(LuOut,*) ' DryRun ',DryRun
319
320          write(LuOut,*) 'ALPHAp:  ',i_ALPHAp
321          write(LuOut,*) 'IPAIRp:  ',i_IPAIRp
322          write(LuOut,*) 'Ep    :  ',i_Ep
323          write(LuOut,*) 'R0    :  ',i_R0
324          write(LuOut,*) 'IJK   :  ',i_IJK
325          write(LuOut,*) 'P     :  ',i_P
326          write(LuOut,*) 'RS    :  ',i_RS
327          write(LuOut,*) 'PC    :  ',i_PC
328          write(LuOut,*) 'ff    :  ',i_ff
329          write(LuOut,*) 'Rj    :  ',i_Rj
330
331          call errquit('hf1: insufficient memory for hfmkr',911,
332     &       MEM_ERR)
333
334        else
335
336          call hf1mkr(Axyz,Bxyz,Cxyz,zan,exinv,ncenters,
337     &        W0(i_ALPHAp),W0(i_P),W0(i_RS),W0(i_PC),W0(i_ff),
338     &        W0(i_Rj),W0(i_R0),W0(i_R0),W0(i_IJK),
339     &        NPP,Lp,Lp3,.FALSE.)
340
341          if (gencon) then
342            call hfnai_gc(W0(i_Ep),W0(i_R0),W0(i_IJK),bNAI,
343     &          W0(i_prim_ints),W0(i_half_ints),
344     &          Acoefs,Bcoefs,W0(i_IPAIRp),
345     &          NPA,NPB,NCA,NCB,NPP,La,Lb,La2,Lb2,Li,Lp,Lp3,canAB)
346          else
347            call hfnai(W0(i_Ep),W0(i_R0),W0(i_IJK),bNAI,
348     &          Nints,NPP,La,Lb,Li,Lp,Lp3,canAB)
349          end if
350
351        end if
352
353      end if
354
355c Return the maximum amount of scratch space required by a "dry run".
356
357      if( DryRun ) maxW0 = MaxMem
358c
359      end
360