1      Subroutine hf1d_cosmo(
2     &       Axyz,Aprims,Acoefs,NPA,NCA,La,ictrA,
3     &       Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,ictrB,
4     &       Cxyz,zan,
5     &       bNAI,Nint,
6     &       canAB,dryrun,W0,maxW0)
7c $Id$
8c
9c     Calculate derivatives of an nuclear attraction type integral.
10c     The integral is:
11c
12c       E = int(dr;Chi_A(r)Chi_B(r)Q_C/|r-R_C|)
13c
14c     We are interested in the derivatives with respect to the centers
15c     A, B, and C. So upon return we want
16c
17c       bNAI(1:nint,1:3) = dE/dA                  (1)
18c
19c       bNAI(1:nint,4:6) = dE/dB                  (2)
20c
21c       bNAI(1:nint,7:9) = dE/dC = -dE/dA-dE/dB   (3)
22c
23c    Typically the calculation of derivative integrals relies on the
24c    fact that the derivative of a Gaussian function with respect to
25c    its center position is again a Gaussian function but of a higher
26c    angular momentum. Hence the same code can be used to calculate
27c    derivative integrals as well as energy integrals (I am not sure
28c    whether the same would be true if directly calculating the
29c    derivative with respect to C as that would change the operator).
30c
31c    So we expect this code to evaluate eqs. (1) and (2) and use
32c    translational invariance to evaluate eq. (3) (i.e. the second
33c    equality).
34c
35c     Implicit real*8 (a-h,o-z)
36c     Implicit integer (i-n)
37      Implicit none
38
39#include "errquit.fh"
40
41      Logical dryrun
42      Logical canAB,GenCon ! what is canAB???
43      Integer Nint, maxW0
44
45c--> Cartesian Coordinates, Primitives & Contraction Coefficients
46
47      Double precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
48      Double precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
49
50c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents
51
52      Double precision Cxyz(3), zan
53
54c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals
55
56      Double precision bNAI(Nint,9)
57
58c--> Derivative Indices
59
60      Integer inder1(3,3),inder2(6,6)
61
62c--> Scratch Space.
63
64      Double precision W0(maxW0)
65
66c--> Other stuff
67
68      Integer lprod, mxd, Lp3, Lp, Nder, NPP
69      Integer NPA, NCA, NPB, NCB, NCP
70      Integer La, Lb, Li, ictrA, ictrB
71      Integer i_VA, i_VB, i_VP, i_VR, i_Ri, i_Rj, i_RS
72      Integer MaxMem, nd
73      Integer i_P, i_R0, i_top, i_pf, i_PC, i_R0C, i_left, i_right
74      Integer i_IPAIRp, i_ESp, i_ALPHAp, i_Ep
75      Integer i_IJK, i_ff
76
77      Data inder1/ 1,0,0,
78     &             0,1,0,
79     &             0,0,1 /
80
81      Data inder2/ 1,0,0,0,0,0,
82     &             0,1,0,0,0,0,
83     &             0,0,1,0,0,0,
84     &             0,0,0,1,0,0,
85     &             0,0,0,0,1,0,
86     &             0,0,0,0,0,1 /
87c
88c Compute gradient of the overlap, kinetic energy, and nuclear attraction
89c integrals for 2 shells of contracted Gaussians functions.
90c
91c The derivative indices specify the order of differentiation for each
92c coordinate:
93c
94c      inder1(i,n), i=1-3:   d/dRx, d/dRy, d/dRz
95c
96c      inder2(i,n), i=1-3:   d/dPx, d/dPy, d/dPz
97c                   i=4-6:   d/dRx, d/dRy, d/dRz
98c
99c******************************************************************************
100
101#if defined(INTDEBUG)
102      call hf_print_set(1)
103      call hf_print('hf1d: a shell',axyz,aprims,acoefs,npa,nca,la)
104      call hf_print('hf1d: b shell',bxyz,bprims,bcoefs,npb,ncb,lb)
105      call hf_print_set(0)
106#endif
107      MXD = 1
108      Li  = 0
109
110c Determine whether general or segmented contraction is used.
111
112      NCP = NCA*NCB
113
114      GenCon = NCP.ne.1
115
116      if( GenCon )then
117        write(*,*) 'HF1D: Not yet ready for general contraction.'
118        call errquit('HF1D: Not yet ready for general contraction.',
119     &               0,UERR)
120      end if
121
122c Define the angular momentum of the overlap distribution.
123
124      Lp = La + Lb
125
126c Increment "Lp" to account for the order of differentiation.
127
128      Lp = Lp + MXD
129
130c Define the accumulated number of angular momentum functions <= Lp.
131
132      Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6
133
134c Define the prefactor of the overlap distribution "P".
135
136c Assign pointers to scratch space.
137
138      i_ALPHAp = 1
139      i_IPAIRp = i_ALPHAp + 2*(NPA*NPB)
140      i_left   = i_IPAIRp + 2*(NPA*NPB) - 1
141
142      i_ESp   = (maxW0+1) - 3*(NPA*NPB)
143      i_right = i_ESp
144
145      if( i_left.ge.i_right )then
146
147       write(*,*) 'HF1D 1:Insufficient scratch space.'
148       write(*,*) '        needed    ',i_left + (maxW0-(i_right-1))
149       write(*,*) '        allocated ',maxW0
150
151       write(*,*) 'From the left '
152       write(*,*) 'ALPHAp:  ',i_ALPHAp
153       write(*,*) 'IPAIRp:  ',i_IPAIRp
154       write(*,*) 'From the right '
155       write(*,*) 'ESp   :  ',i_ESp
156
157       call errquit('HF1D:  Insufficient scratch space.',i_right-i_left,
158     &              UERR)
159
160      end if
161
162      MaxMem = 1    ! take care of compiler warnings
163      if (dryrun) then
164        MaxMem = i_left + (maxW0 - (i_right-1))
165        NPP = NPA*NPB
166      else
167        call hfset(Axyz,Aprims,Acoefs,NPA,NCA,
168     &             Bxyz,Bprims,Bcoefs,NPB,NCB,
169     &             GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP)
170      endif
171
172c Define the Hermite linear expansion coefficients.
173
174c Assign pointers to scratch space.
175
176      lprod = ((La+Li)+(Lb+Li)+1)*((La+Li)+1)*((Lb+Li)+1)
177
178      i_Ep   = i_IPAIRp + 2*(NPA*NPB)
179      i_pf   = i_Ep     + 3*NPP*(MXD+1)*lprod
180      i_left = i_pf     + 2*NPP - 1
181
182      if( i_left.ge.i_right )then
183
184       write(*,*) 'HF1D 2:Insufficient scratch space.'
185       write(*,*) '        needed    ',i_left + (maxW0-(i_right-1))
186       write(*,*) '        allocated ',maxW0
187
188       write(*,*) 'From the right '
189       write(*,*) 'ALPHAp:  ',i_ALPHAp
190       write(*,*) 'IPAIRp:  ',i_IPAIRp
191       write(*,*) 'Ep    :  ',i_Ep
192       write(*,*) 'pf    :  ',i_pf, npp,mxd,lprod
193       write(*,*) 'From the left '
194       write(*,*) 'ESp   :  ',i_ESp
195
196       call errquit('HF1D:  Insufficient scratch space.',i_right-i_left,
197     &              UERR)
198
199      end if
200
201      if (dryrun) then
202        MaxMem = max(MaxMem, (i_left+(maxW0-(i_right-1))))
203      else
204        do 100 nd = 0,MXD
205
206          call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep),
207     &           W0(i_pf),nd,NPP,MXD,La+Li,Lb+Li)
208
209100     continue
210      endif
211
212
213
214c Compute nuclear attraction integrals, <a|V|b>.
215
216c     if( NAI )then
217
218c Define the auxiliary function integrals.
219
220c Assign scratch space.
221
222       i_R0  = i_Ep  + 3*NPP*(MXD+1)*lprod
223       i_IJK = i_R0  + NPP*Lp3
224       i_R0C = i_IJK + (Lp+1)**3
225       i_P   = i_R0C + NPP*Lp3*1
226       i_RS  = i_P   + NPP*3
227       i_PC  = i_RS  + NPP
228       i_ff  = i_PC  + NPP*3
229       i_Rj  = i_ff  + NPP*2
230       i_top = i_Rj  + NPP*(Lp+1)*Lp3 - 1
231
232       if( i_top.gt.maxW0 )then
233
234        write(*,*) 'HF1D 3:Insufficient scratch space.'
235        write(*,*) '        needed    ',i_top
236        write(*,*) '        allocated ',maxW0
237
238        write(*,*) 'ALPHAp:  ',i_ALPHAp
239        write(*,*) 'IPAIRp:  ',i_IPAIRp
240        write(*,*) 'Ep    :  ',i_Ep
241        write(*,*) 'R0    :  ',i_R0
242        write(*,*) 'IJK   :  ',i_IJK
243        write(*,*) 'R0C   :  ',i_R0C
244        write(*,*) 'P     :  ',i_P
245        write(*,*) 'RS    :  ',i_RS
246        write(*,*) 'PC    :  ',i_PC
247        write(*,*) 'ff    :  ',i_ff
248        write(*,*) 'Rj    :  ',i_Rj
249
250        call errquit('HF1D:  Insufficient scratch space.',maxW0-i_top,
251     &               UERR)
252
253       end if
254
255       if (dryrun) then
256         MaxMem = max(MaxMem, i_top)
257       else
258         call hf1mkr(Axyz,Bxyz,Cxyz,zan,0.0d0,1,
259     &          W0(i_ALPHAp),W0(i_P),W0(i_RS),W0(i_PC),W0(i_ff),
260     &          W0(i_Rj),W0(i_R0),W0(i_R0C),W0(i_IJK),
261     &          NPP,Lp,Lp3,.TRUE.)
262c
263c        Compute the 1-electron - nuclear attraction contributions of
264c        the gradient wrt the nuclear coordinates
265c
266c        call hfefi(W0(i_Ep),W0(i_R0C),W0(i_IJK),bNAI,
267c    &          NPP,Nint,La,Lb,Li,Lp,Lp3,1,
268c    &          MXD,canAB,ictra,ictrb)
269
270       endif
271c Compute compute the derivative wrt to the centers of the Gaussian functions.
272
273c Assign scratch space.
274
275       i_VP  = i_IJK + (Lp+1)**3
276       i_VR  = i_VP  + NPP*(Nint*3)
277       i_VA  = i_VR  + NPP*(Nint*3)
278       i_VB  = i_VA  + (Nint*3)
279       i_ff  = i_VB  + (Nint*3)
280       i_top = i_ff  + NPP*2 - 1
281
282       if( i_top.gt.maxW0 )then
283
284        write(*,*) 'HF1D 4:Insufficient scratch space.'
285        write(*,*) '        needed    ',i_top
286        write(*,*) '        allocated ',maxW0
287
288        write(*,*) 'ALPHAp:  ',i_ALPHAp
289        write(*,*) 'IPAIRp:  ',i_IPAIRp
290        write(*,*) 'Ep    :  ',i_Ep
291        write(*,*) 'R0    :  ',i_R0
292        write(*,*) 'IJK   :  ',i_IJK
293        write(*,*) 'VP    :  ',i_VP
294        write(*,*) 'VR    :  ',i_VR
295        write(*,*) 'VA    :  ',i_VA
296        write(*,*) 'VB    :  ',i_VB
297        write(*,*) 'ff    :  ',i_ff
298
299        call errquit('HF1D:  Insufficient scratch space.',maxW0-i_top,
300     &               UERR)
301
302       end if
303
304       if (dryrun) then
305         MaxMem = max(MaxMem, i_top)
306       else
307c Compute the derivatives wrt to (P,R).
308
309c Compute the derivatives of the primitive integrals.
310
311         Nder = 6
312
313         call hfdnai(W0(i_Ep),W0(i_R0),W0(i_IJK),W0(i_VP),
314     &          NPP,Nint,La,Lb,Li,Lp,Lp3,1,
315     &          MXD,inder2,Nder,canAB)
316
317c Transform to obtain derivatives wrt (A,B).
318
319         call hf1PRtoAB(W0(i_VP),W0(i_VR),W0(i_VA),W0(i_VB),
320     &          W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ff),NPP,Nint*3,
321     &          ictrA,ictrB)
322
323c Combine the two components of the derivatives of NAIs.
324c
325c   (1) the derivative wrt to each of the nuclear attraction centers,
326c   (2) the derivative wrt to the centers of the Gaussian functions.
327
328         if( ictrA.eq.ictrB )then
329           call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,1),1)
330           call daxpy(Nint*3,-1.D0,W0(i_VA),1,bNAI(1,7),1)
331         else
332           call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,1),1)
333           call daxpy(Nint*3,1.D0,W0(i_VB),1,bNAI(1,4),1)
334           call daxpy(Nint*3,-1.D0,W0(i_VA),1,bNAI(1,7),1)
335           call daxpy(Nint*3,-1.D0,W0(i_VB),1,bNAI(1,7),1)
336         end if
337
338       endif
339c     end if
340
341c Return the maximum amount of scratch space required by a "dry run".
342
343      if( DryRun ) maxW0 = MaxMem
344c
345      end
346