1      subroutine hf3pot(
2     &    Axyz,Aprims,Acoef,NPA,NCA,La,
3     &    Bxyz,Bprims,Bcoef,NPB,NCB,Lb,
4     &    Cxyz,Cprims,Ccoef,NPC,NCC,Lc,
5     &    b3nai,b3nai_sz,nint,
6     &    DryRun,Scr,lscr)
7c
8c $Id$
9c
10      implicit none
11c:: includes
12#include "stdio.fh"
13#include "mafdecls.fh"
14#include "errquit.fh"
15c:: passed
16      integer La  ! [input] angular momentum for center a
17      integer Lb  ! [input] angular momentum for center b
18      integer Lc  ! [input] angular momentum for center c
19      integer NPA ! [input] number of primitives on center a
20      integer NPB ! [input] number of primitives on center b
21      integer NPC ! [input] number of primitives on center c
22      integer NCA ! [input] number of general contractions on center a
23      integer NCB ! [input] number of general contractions on center b
24      integer NCC ! [input] number of general contractions on center c
25      double precision Axyz(3) ! [input] coordinates for center a
26      double precision Bxyz(3) ! [input] coordinates for center b
27      double precision Cxyz(3) ! [input] coordinates for center c
28      double precision Aprims(NPA) ! [input] primitive exponents on center a
29      double precision Bprims(NPB) ! [input] primitive exponents on center b
30      double precision Cprims(NPC) ! [input] primitive exponents on center c
31      double precision Acoef(NPA,NCA) ! [input] primitive coefficients on center a
32      double precision Bcoef(NPB,NCB) ! [input] primitive coefficients on center b
33      double precision Ccoef(NPC,NCC) ! [input] primitive coefficients on center c
34      integer b3nai_sz                 ! [input] input size of 3 center nai buffer
35      double precision b3nai(b3nai_sz) ! [output] 3 center nai buffer
36      integer nint                     ! [output] actual number of integrals computed
37      logical DryRun                   ! [input] run routine only to test memory size needed
38      integer lscr                     ! [input] size of scratch array
39      double precision scr(lscr)       ! [scratch] scratch array for computations
40c::local
41      logical gencon ! general contraction flag
42      Integer Nabc   ! product of primitives
43      Integer maxmem ! memory returned in lscr for dryrun
44      Integer Lg     ! total angular momentum
45      Integer Lg3    ! l3 of total angular momentum
46c:: scratch pointers
47      Integer i_top   ! top of scratch
48      Integer i_save  ! point at which scratch can be reused
49      Integer i_alpha ! scr ptr for exponents and Hermite expansion coef prefactor
50      Integer i_E     ! Hermite expansion coefficients
51      Integer i_G1    ! center G coordinates (first use;through making of E)
52      Integer i_G2    ! center G coordinates (second use; making of R)
53      Integer i_Gt    ! Temp array for G coorinates
54      Integer i_ABC2I ! Temp array for prefactors of E
55      Integer i_R0    ! the auxillary hermite integrals required for sum
56      Integer i_ijk   ! pointer for auxillary hermite integrals
57      Integer i_RS    ! prefactors for auxillary hermite integrals
58      Integer i_GC    ! center for auxillary hermite integral recursion
59      Integer i_r     ! the auxillary hermite integrals required for recursion
60      Integer i_ff    ! prefactor array for forming recusion set of R's
61
62#if defined(INTDEBUG)
63      write(luout,*)' hf3pot debug'
64      if (.not.dryrun) then
65        write(luout,*)'b3nai_sz',b3nai_sz
66        write(luout,*)'dryrun',dryrun
67        write(luout,*)'lscr',lscr
68        call hf_print_set(1)
69        call hf_print('hf3pot a',Axyz,Aprims,Acoef,NPA,NCA,La)
70        call hf_print('hf3pot b',Bxyz,Bprims,Bcoef,NPB,NCB,Lb)
71        call hf_print('hf3pot c',Cxyz,Cprims,Ccoef,NPC,NCC,Lc)
72        call hf_print_set(0)
73      endif
74#endif
75c
76c Compute 3-center nuclear attraction integrals for ECPs
77c
78      gencon = nca.gt.1.or.ncb.gt.1.or.ncc.gt.1
79      if (gencon) call errquit
80     &    (' hf3pot not ready for general contractions ',911,
81     &       INT_ERR)
82c
83      i_alpha = 1
84      i_top   = i_alpha + (NPA*NPB*NPC)*4 - 1
85
86#if defined(INTDEBUG)
87      write(luout,*)'i_alpha :',i_alpha
88      write(luout,*)'i_top   :',i_top
89      write(luout,*)npa,npb,npc,(NPA*NPB*NPC),(NPA*NPB*NPC)*4
90#endif
91      if (i_top.gt.lscr) then
92        write(luout,*)' hf3pot: insufficient scratch space '
93        write(luout,*)'       : needed    :',i_top
94        write(luout,*)'       : allocated :',lscr
95        call errquit('hf3pot: fatal error ',1, INT_ERR)
96      endif
97c
98c.. determine actual product pairs to be kept.
99c   The following line is to take care of compiler warnings.
100c
101      maxmem = i_top
102      if (DryRun) then
103        maxmem = i_top
104        Nabc = NPA*NPB*NPC
105      else
106        call hf1set3(
107     &      Axyz,Aprims,Acoef,NPA,
108     &      Bxyz,Bprims,Bcoef,NPB,
109     &      Cxyz,Cprims,Ccoef,NPC,
110     &      scr(i_alpha),Nabc)
111      endif
112#if defined(INTDEBUG)
113      write(luout,*)
114      write(luout,*)' nabc = ',nabc
115      write(luout,*)
116#endif
117
118c Define the center of the charge distribution.
119c A+B -> P; P + C -> G
120
121      Lg = La+Lb+Lc ! total angular momentum
122
123c
124      i_E    = i_top + 1
125      i_G1   = i_E   + Nabc*3*(Lg+1)*(La+1)*(Lb+1)*(Lc+1)
126      i_save = i_G1  ! i_alpha and I_E saved
127      i_top  = i_G1  + Nabc*3 - 1
128
129#if defined(INTDEBUG)
130      write(luout,*)
131      write(luout,*)'i_alpha :',i_alpha
132      write(luout,*)'i_E     :',i_E
133      write(luout,*)'i_G1    :',i_G1
134      write(luout,*)'i_save  :',i_save
135      write(luout,*)'i_top   :',i_top
136      write(luout,*)
137#endif
138
139      if (i_top.gt.lscr) then
140        write(luout,*)' hf3pot: insufficient scratch space '
141        write(luout,*)'       : needed    :',i_top
142        write(luout,*)'       : allocated :',lscr
143
144        write(luout,*)'i_alpha :',i_alpha
145        write(luout,*)'i_E     :',i_E
146        write(luout,*)'i_G1    :',i_G1
147        call errquit('hf3pot: fatal error ',2, INT_ERR)
148      endif
149c
150      if (DryRun) then
151        maxmem = max(maxmem, i_top)
152      else
153        call hfctr3(Axyz,Bxyz,Cxyz,scr(i_alpha),scr(i_G1),Nabc)
154      endif
155c
156c Define the Hermite linear expansion coefficients
157c
158      i_Gt    = i_top   + 1
159      i_ABC2I = i_Gt    + Nabc*3
160      i_top   = i_ABC2I + Nabc*3 - 1
161
162#if defined(INTDEBUG)
163      write(luout,*)
164      write(luout,*)'i_save  :',i_save
165      write(luout,*)'i_alpha :',i_alpha
166      write(luout,*)'i_E     :',i_E
167      write(luout,*)'i_G1    :',i_G1
168      write(luout,*)'i_Gt    :',i_Gt
169      write(luout,*)'i_ABC2I :',i_ABC2I
170      write(luout,*)'i_top   :',i_top
171      write(luout,*)
172#endif
173
174      if (i_top.gt.lscr) then
175        write(luout,*)' hf3pot: insufficient scratch space '
176        write(luout,*)'       : needed    :',i_top
177        write(luout,*)'       : allocated :',lscr
178
179        write(luout,*)'i_alpha :',i_alpha
180        write(luout,*)'i_E     :',i_E
181        write(luout,*)'i_G1    :',i_G1
182        write(luout,*)'i_Gt    :',i_Gt
183        write(luout,*)'i_ABC2I :',i_ABC2I
184        call errquit('hf3pot: fatal error ',3, INT_ERR)
185      endif
186c
187      if (DryRun) then
188        maxmem = max(maxmem, i_top)
189      else
190        call hf1mke3(Axyz, Bxyz, Cxyz,
191     &      scr(i_alpha), scr(i_G1), scr(i_Gt), scr(i_ABC2I),
192     &      scr(i_E), Nabc, La, Lb, Lc)
193      endif
194c
195      Lg3   = ((Lg+1)*(Lg+2)*(Lg+3))/6
196
197#if defined(INTDEBUG)
198      write(luout,*)
199      write(luout,*)'lg,lg3',lg,lg3
200      write(luout,*)
201#endif
202
203      i_R0  = i_save
204      i_IJK = i_R0   + Nabc*Lg3
205      i_RS  = i_IJK  + (Lg+1)*(Lg+1)*(Lg+1)
206      i_GC  = i_RS   + Nabc
207      i_G2  = i_GC   + Nabc*3
208      i_ff  = i_G2   + Nabc*3
209      i_R   = i_ff   + 2*Nabc
210      i_top = i_R    + Nabc*(Lg+1)*Lg3 - 1
211
212#if defined(INTDEBUG)
213      write(luout,*)
214      write(luout,*)'i_save  :',i_save
215      write(luout,*)'i_alpha :',i_alpha
216      write(luout,*)'i_E     :',i_E
217      write(luout,*)'i_R0    :',i_R0
218      write(luout,*)'i_IJK   :',i_IJK
219      write(luout,*)'i_RS    :',i_RS
220      write(luout,*)'i_GC    :',i_GC
221      write(luout,*)'i_G2    :',i_G2
222      write(luout,*)'i_ff    :',i_ff
223      write(luout,*)'i_R     :',i_R
224      write(luout,*)'i_top   :',i_top
225      write(luout,*)
226#endif
227
228      if (i_top.gt.lscr) then
229        write(luout,*)' hf3pot: insufficient scratch space '
230        write(luout,*)'       : needed    :',i_top
231        write(luout,*)'       : allocated :',lscr
232
233        write(luout,*)'i_alpha :',i_alpha
234        write(luout,*)'i_E     :',i_E
235        write(luout,*)'i_R0    :',i_R0
236        write(luout,*)'i_IJK   :',i_IJK
237        write(luout,*)'i_RS    :',i_RS
238        write(luout,*)'i_GC    :',i_GC
239        write(luout,*)'i_G2    :',i_G2
240        write(luout,*)'i_ff    :',i_ff
241        write(luout,*)'i_R     :',i_R
242        call errquit('hf3pot: fatal error ',3, INT_ERR)
243      endif
244c
245      if (DryRun) then
246        maxmem = max(maxmem, i_top)
247      else
248        call hf3mkr(Axyz, Bxyz, Cxyz, scr(i_alpha), scr(i_G2),
249     &      scr(i_RS), scr(i_GC), scr(i_ff), scr(i_R), scr(i_R0),
250     &      scr(i_IJK), Nabc, Lg, Lg3)
251      endif
252c
253      if (DryRun) then
254        lscr = maxmem
255        return
256      endif
257c
258      call hf3PEabc(b3nai, scr(i_E), scr(i_R0), scr(i_IJK),
259     &    Nabc, La, Lb, Lc, Lg, Lg3, nint, b3nai_sz)
260c
261      end
262