1      subroutine dim_elfder(xyzi,expi,coefi,i_nprim,i_ngen, Li,
2     1  xyzj,expj,coefj,j_nprim,j_ngen,Lj,nder,nint,elfder,scr,lscr,
3     2  xyzpt, npt)
4c
5c $Id: hnd_elfder.F 21348 2011-10-31 22:51:25Z d3p852 $
6c
7      implicit double precision (a-h,o-z)
8#include "hnd_pointers.fh"
9      dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
10      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
11      dimension xyzpt(3,npt)
12      dimension scr(lscr)
13      dimension elfder(*)
14c
15c
16c     ----- Wrapper routine that sets the sizes of scratch blocks -----
17c
18      idim = max(nder*3,1)
19c
20      call dim_elfder1(xyzi,expi,coefi,i_nprim,i_ngen, Li, xyzj,expj,
21     1 coefj, j_nprim, j_ngen, Lj, nder, nint, elfder, xyzpt, npt, idim,
22     2 scr(elpt(1)),scr(elpt(2)),scr(elpt(3)),scr(elpt(4)),scr(elpt(5)),
23     3 scr(elpt(6)),scr(elpt(7)),scr(elpt(8)),scr(elpt(9)))
24c
25      return
26      end
27c
28      subroutine dim_elfder1(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,
29     1  expj,coefj, j_nprim, j_ngen, Lj, nder, nint, elfder, xyzpt, npt,
30     2  idim,xv,yv,zv,dxv,dyv,dzv,ddxv,ddyv,ddzv)
31c
32      implicit double precision (a-h,o-z)
33#include "nwc_const.fh"
34#include "hnd_rys.fh"
35#include "stdio.fh"
36#include "hnd_tol.fh"
37#include "dimqm.fh"
38      common/hnd_xyzder/xint,yint,zint,tx,x0,y0,z0,xi,yi,zi,xj,yj,zj,
39     1                  ni,nj,cx,cy,cz
40      dimension w2(maxrys),w4(maxrys)
41      dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
42      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
43      dimension xyzpt(3,npt),elfder(nint,idim,npt)
44c
45      dimension xv(Li+1,Lj+1,*)
46      dimension yv(Li+1,Lj+1,*)
47      dimension zv(Li+1,Lj+1,*)
48      dimension dxv(Li+1,Lj+1,*)
49      dimension dyv(Li+1,Lj+1,*)
50      dimension dzv(Li+1,Lj+1,*)
51      dimension ddxv(Li+1,Lj+1,*)
52      dimension ddyv(Li+1,Lj+1,*)
53      dimension ddzv(Li+1,Lj+1,*)
54c
55c     Routine calculates the following integrals for a given shell:
56c
57c     nder =-2 : Electronic wave function
58c     nder =-1 : Electronic density
59c     nder = 0 : Electrostatic potential
60c     nder = 1 : Electric field
61c     nder = 2 : Electric field gradient
62c
63c     Maximum scratch size needs to be Max(nder*3,1))*(Li+1)*(Lj+1)*((Li+Lj)/2+1)
64c     Data format elfder will be (nints,ipts,Max(nder*3,1))
65c
66      data rln10  /2.30258d+00/
67      data pi212  /1.1283791670955d+00/  ! 2/sqrt(pi)
68c
69      dtol=rln10*itol
70c
71c     Zero integral array
72c
73      call dcopy(nint*max(nder*3,1)*npt,0.0d0,0,elfder,1)
74      xintmax = 0.0d0
75      yintmax = 0.0d0
76      zintmax = 0.0d0
77c
78c     ----- ishell -----
79c
80      xi=xyzi(1)
81      yi=xyzi(2)
82      zi=xyzi(3)
83      lit = Li + 1
84      maxi = lit*(lit+1)/2
85c
86c     ----- jshell -----
87c
88      xj=xyzj(1)
89      yj=xyzj(2)
90      zj=xyzj(3)
91      ljt = Lj + 1
92      maxj = ljt*(ljt+1)/2
93c      write(luout,*) "i shell:", xi, yi, zi
94c      write(luout,*) "j shell:", xj, yj, zj
95c
96      rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
97      nroots=(lit+ljt+nder-2)/2+1
98
99      if(nroots.gt.maxrys) then
100         write(luout,9997) maxrys,lit,ljt,nroots
101         call errquit('hnd_elfgrd: need higher Rys root',nroots,INT_ERR)
102      endif
103c
104c     ----- i primitive -----
105c
106      do 7000 ig=1,i_nprim
107      ai=expi(ig)
108      arri=ai*rr
109      axi=ai*xi
110      ayi=ai*yi
111      azi=ai*zi
112      csi=coefi(ig,i_ngen)
113c      write(luout,*) "Li:", Li
114c      write(luout,*) "ai:", ai
115c
116c     ----- j primitive -----
117c
118      do 6000 jg=1,j_nprim
119      aj=expj(jg)
120c      write(luout,*) "Lj:", Lj
121c      write(luout,*) "aj:", aj
122      aa=ai+aj
123c      write(luout,*) "aa:", aa
124      aa1=1.0d0/aa
125      dum=aj*arri*aa1
126      if(dum.gt.dtol) go to 6000
127      fac= exp(-dum)
128      csj=coefj(jg,j_ngen)
129      ax=(axi+aj*xj)*aa1
130      ay=(ayi+aj*yj)*aa1
131      az=(azi+aj*zj)*aa1
132c
133c     ----- density factor -----
134c
135      dum1 = csi * fac
136      dij = dum1 * csj
137c
138c     Electronic density integrals and electronic wave
139c     function are done differently
140c     then the electric field, electric field gradient, and
141c     electrostatic potential integrals
142c     Switch between the two classes here:
143c
144      if (nder.ge.0) goto 399
145c
146c     ----- electronic wave function / density integral -----
147c
148c     For wave function we simply set j to 0.0d0 and reset the
149c     density factor appropriately
150c
151      if (nder.eq.-2) then
152         dij = csi
153         ax = xi
154         ay = yi
155         az = zi
156         aa = ai
157         ljt = 1
158         maxj = 1
159      endif
160c
161      do 300 ipt=1,npt
162         x0=xyzpt(1,ipt)
163         y0=xyzpt(2,ipt)
164         z0=xyzpt(3,ipt)
165         dum = aa*((x0-ax)**2+(y0-ay)**2+(z0-az)**2)
166         if(dum.gt.dtol) go to 300
167         fac = exp(-dum)
168c
169c        ----- density values -----
170c
171         do 270 j=1,ljt
172            nj=j
173            do 270 i=1,lit
174               ni=i
175               call hnd_denxyz
176               xv(i,j,1)=xint
177               yv(i,j,1)=yint
178               zv(i,j,1)=zint
179  270    continue
180c
181c        ----- combining the pieces together -----
182c
183         ij=0
184         do 290 j=1,maxj
185            call getNxyz(Lj,j,Nxyz)
186            jx = Nxyz(1) + 1
187            jy = Nxyz(2) + 1
188            jz = Nxyz(3) + 1
189            do 290 i=1,maxi
190               call getNxyz(Li,i,Nxyz)
191               ix = Nxyz(1) + 1
192               iy = Nxyz(2) + 1
193               iz = Nxyz(3) + 1
194               ij=ij+1
195               elfder(ij,1,ipt)=elfder(ij,1,ipt)+fac*dij*xv(ix,jx,1)*
196     1                          yv(iy,jy,1)*zv(iz,jz,1)
197  290    continue
198  300 continue
199c
200      goto 6000
201c
202c     ----- electric field (gradient) term -----
203c
204  399 dij = dij * pi212 * aa1
205      aax=aa*ax
206      aay=aa*ay
207      aaz=aa*az
208c      write(luout,*) "dij:", dij
209c      aa = aa/erf(aa)
210      do 500 ipt=1,npt
211         znuc=1.0d0
212         cx=xyzpt(1,ipt)
213         cy=xyzpt(2,ipt)
214         cz=xyzpt(3,ipt)
215         cu = 1.0d0/(2.7306542620197547)**2
216c         xci = exp(-ai*(xi-cx)**2/(ai + 1))
217c         yci = exp(-ai*(yi-cy)**2/(ai + 1))
218c         zci = exp(-ai*(zi-cz)**2/(ai + 1))
219c         xcj = exp(-aj*(xj-cx)**2/(aj + 1))
220c         ycj = exp(-aj*(yj-cy)**2/(aj + 1))
221c         zcj = exp(-aj*(zj-cz)**2/(aj + 1))
222c              xci = exp(-(aa*cu/(aa+cu))*(ax-cx)**2)
223c              yci = exp(-(aa*cu/(aa+cu))*(ay-cy)**2)
224c              zci = exp(-(aa*cu/(aa+cu))*(az-cz)**2)
225c            temp = xci*yci*zci
226c            damp = erfc(temp)
227c            write(luout,*) xci, yci, zci
228c            write(luout,*) xci*yci*zci
229c             if(temp .ge. 1.0d-2) then
230c               write(luout,*) "overlap:", temp
231c               write(luout,*) "scale: " , erfc(temp)
232c             end if
233
234c         write(luout,*) ipt, 'ci'
235c         write(luout,*) xci, yci, zci
236c         write(luout,*) ipt, 'cj'
237c         write(luout,*) xcj, ycj, zcj
238c
239c        Testing DIM atom as Gaussian
240c         rr = (ax-cx)**2+(ay-cy)**2+(az-cz)**2
241c         aac = aa+ac
242c         dum = aa*rr*ac
243c         dijc = dij * exp(-dum/aac)
244
245c         temp = sqrt((ax-cx)**2+(ay-cy)**2+(az-cz)**2)
246c         if(temp .le. 10.0) write(luout,*) temp
247         yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2)
248c         if(yy .le. 20.0) write(luout,*) "yy:", yy
249c         if(yy .le. 20.0) write(luout,*) "aa:", aa
250c         if(yy .le. 20.0) write(luout,*) "temp:", temp
251c         if(yy .le. 20.0) write(luout,*) "nroots:", nroots
252         call hnd_droot
253         do 420 iroot=1,nroots
254            uu=u9(iroot)*aa
255c            write(luout,*) "uu:", uu
256c            uu = erfc(uu) * uu
257c             uu = 1 - exp(-uu)
258c            ddfac = 1.0d0
259c            if(uu .ge. 0.1d0) then
260c              write(luout,*) "uu:", uu
261c              ddfac = 0.0d0
262c            end if
263            u4=uu*uu
264            ww=w9(iroot)*znuc
265            w2(iroot)=ww*uu*2.0d0
266            w4(iroot)=ww*u4*4.0d0
267            w9(iroot)=ww
268            tt=1.0d0/(aa+uu)
269            tx= sqrt(tt)
270            x0=(aax+uu*cx)*tt
271            y0=(aay+uu*cy)*tt
272            z0=(aaz+uu*cz)*tt
273c              xci = exp(-(aa*uu*tt)*(ax-cx)**2)
274c              yci = exp(-(aa*uu*tt)*(ay-cy)**2)
275c              zci = exp(-(aa*uu*tt)*(az-cz)**2)
276c              temp = erfc(xci*yci*zci)
277c              if(temp .le. 0.5) then
278c                write(luout,*) "overlap:", xci*yci*zci
279c              end if
280c            write(luout,*) x0, y0, z0
281c            temp = sqrt((x0-xi)**2+(y0-yi)**2+(z0-zi)**2)
282c            if(temp .le. 10.0) write(luout,*) "Ri0:", temp
283c            temp = sqrt((x0-xj)**2+(y0-yj)**2+(z0-zj)**2)
284c            if(temp .le. 10.0) write(luout,*) "Rj0:", temp
285c            write(luout,*) temp
286
287            do 410 i=1,lit
288               ni=i
289               do 410 j=1,ljt
290                  nj=j
291                  goto (402,401) nder+1
292                  call dim_dervxyz(2)
293                  ddxv(i,j,iroot)=xint
294                  ddyv(i,j,iroot)=yint
295                  ddzv(i,j,iroot)=zint
296  401             call dim_dervxyz(1)
297                  dxv(i,j,iroot)=xint
298                  dyv(i,j,iroot)=yint
299                  dzv(i,j,iroot)=zint
300  402             call dim_sxyz
301                  xv(i,j,iroot)=xint
302                  yv(i,j,iroot)=yint
303                  zv(i,j,iroot)=zint
304c                  write(luout,*) xint, yint, zint
305  410       continue
306  420    continue
307c
308c        ----- combining the pieces together -----
309c
310         ij=0
311         do 440 j=1,maxj
312            call getNxyz(Lj,j,Nxyz)
313            jx = Nxyz(1) + 1
314            jy = Nxyz(2) + 1
315            jz = Nxyz(3) + 1
316            do 450 i=1,maxi
317            call getNxyz(Li,i,Nxyz)
318            ix = Nxyz(1) + 1
319            iy = Nxyz(2) + 1
320            iz = Nxyz(3) + 1
321            ij=ij+1
322            if (nder.eq.2) then
323               dumxx=0.0d0
324               dumyy=0.0d0
325               dumzz=0.0d0
326               dumxy=0.0d0
327               dumxz=0.0d0
328               dumyz=0.0d0
329               do 430 iroot=1,nroots
330                  dum=xv(ix,jx,iroot)*yv(iy,jy,iroot)*
331     1                zv(iz,jz,iroot)*w2(iroot)
332                  dumxx=dumxx - dum + ddxv(ix,jx,iroot)*
333     1                  yv(iy,jy,iroot)*zv(iz,jz,iroot)*w4(iroot)
334                  dumyy=dumyy - dum + ddyv(iy,jy,iroot)*
335     1                  xv(ix,jx,iroot)*zv(iz,jz,iroot)*w4(iroot)
336                  dumzz=dumzz - dum + ddzv(iz,jz,iroot)*
337     1                  xv(ix,jx,iroot)*yv(iy,jy,iroot)*w4(iroot)
338                  dumxy=dumxy + dxv(ix,jx,iroot)*dyv(iy,jy,iroot)*
339     1                  zv(iz,jz,iroot)*w4(iroot)
340                  dumxz=dumxz + dxv(ix,jx,iroot)*dzv(iz,jz,iroot)*
341     1                  yv(iy,jy,iroot)*w4(iroot)
342                  dumyz=dumyz + dyv(iy,jy,iroot)*dzv(iz,jz,iroot)*
343     1                  xv(ix,jx,iroot)*w4(iroot)
344  430          continue
345
346               elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumxx*dij
347               elfder(ij,2,ipt) = elfder(ij,2,ipt) + dumyy*dij
348               elfder(ij,3,ipt) = elfder(ij,3,ipt) + dumzz*dij
349               elfder(ij,4,ipt) = elfder(ij,4,ipt) + dumxy*dij
350               elfder(ij,5,ipt) = elfder(ij,5,ipt) + dumxz*dij
351               elfder(ij,6,ipt) = elfder(ij,6,ipt) + dumyz*dij
352            elseif (nder.eq.1) then
353               dumx=0.0d0
354               dumy=0.0d0
355               dumz=0.0d0
356               do 431 iroot=1,nroots
357                  dumx=dumx + dxv(ix,jx,iroot)*yv(iy,jy,iroot)*
358     1                 zv(iz,jz,iroot)*w2(iroot)
359                  dumy=dumy + xv(ix,jx,iroot)*dyv(iy,jy,iroot)*
360     1                 zv(iz,jz,iroot)*w2(iroot)
361                  dumz=dumz + xv(ix,jx,iroot)*yv(iy,jy,iroot)*
362     1                 dzv(iz,jz,iroot)*w2(iroot)
363  431          continue
364               elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumx*dij
365               elfder(ij,2,ipt) = elfder(ij,2,ipt) + dumy*dij
366               elfder(ij,3,ipt) = elfder(ij,3,ipt) + dumz*dij
367            else
368               dumx=0.0d0
369               do 432 iroot=1,nroots
370                  dumx=dumx + xv(ix,jx,iroot)*yv(iy,jy,iroot)*
371     1                 zv(iz,jz,iroot)*w9(iroot)
372  432          continue
373               elfder(ij,1,ipt) = elfder(ij,1,ipt) + dumx*dij
374            endif
375  450       continue
376  440    continue
377c
378  500 continue
379c
380 6000 continue
381 7000 continue
382c
383      return
384 9997 format(' in -elfgrd- , the rys quadrature is not implemented',
385     1       ' beyond -nroots- = ',i3,/,
386     2       ' lit,ljt,nroots = ',3i3)
387      end
388