program inttwoelg
c     two-electron integral - (ab|exp(-eps.r12^2)|cd)
      IMPLICIT NONE
      INTEGER*4 nax,nay,naz,nbx,nby,nbz,ncx,ncy,ncz,ndx,ndy,ndz 
      INTEGER*4 nx,ny,nz,ix,iy,iz,jx,jy,jz
      real*8 ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ai,ain
      real*8 a,b,c,d,xmin,xmax,ymin,ymax,zmin,zmax
      real*8 x1,y1,z1,x2,y2,z2,x1a,x1b,y1a,y1b,z1a,z1b,
     1ee,x2c,x2d,y2c,y2d,z2c,z2d,x12,y12,abcdg,
     1e1,e2,e3,e4,e5,px,py,pz,ksi,ni,R(3),P(3),q(3),
     1abx,aby,abz,cdx,cdy,cdz,GAB,GCD,i0,pi,pi3,eps,
     1w,r2,i0n,g1,g2,g3,g4,g5,th2,t2,ef,we

      open(9,file='TWO.TXT',status='old')
      read(9,*)nax,nay,naz,nbx,nby,nbz,ncx,ncy,ncz,ndx,ndy,ndz
      read(9,*)ax,ay,az
      read(9,*)bx,by,bz
      read(9,*)cx,cy,cz
      read(9,*)dx,dy,dz
      read(9,*)a,b,c,d
      read(9,*)eps
      read(9,*)nx,ny,nz
      close(9)

      call mmx(xmin,xmax,ax,a,bx,b,cx,c,dx,d,px,nx)
      call mmx(ymin,ymax,ay,a,by,b,cy,c,dy,d,py,ny)
      call mmx(zmin,zmax,az,a,bz,b,cz,c,dz,d,pz,nz)
      write(6,601)xmin,xmax,ymin,ymax,zmin,zmax
601   format(6f12.3)

      pi=4.0d0*datan(1.0d0)
      pi3=pi**3
      ksi=a+b
      ni =c+d
      P(1)=(a*ax+b*bx)/(a+b)
      P(2)=(a*ay+b*by)/(a+b)
      P(3)=(a*az+b*bz)/(a+b)
      Q(1)=(c*cx+d*dx)/(c+d)
      Q(2)=(c*cy+d*dy)/(c+d)
      Q(3)=(c*cz+d*dz)/(c+d)
      R(1)=Q(1)-P(1)
      R(2)=Q(2)-P(2)
      R(3)=Q(3)-P(3)
      r2=R(1)**2+R(2)**2+R(3)**2
      abx=ax-bx
      aby=ay-by
      abz=az-bz
      ee=a*b/(a+b)*(abx**2+aby**2+abz**2)
      if(ee.lt.20.0d0)then
       GAB=exp(-ee)
      else
       GAB=0.0d0
      endif
      cdx=cx-dx
      cdy=cy-dy
      cdz=cz-dz
      ee=c*d/(c+d)*(cdx**2+cdy**2+cdz**2)
      if(ee.lt.20.0d0)then
       GCD=exp(-ee)
      else
       GCD=0.0d0
      endif
      th2=ksi*ni/(ksi+ni)
      w=th2/(th2+eps)
      we=eps*th2/(th2+eps)
      ee=w*r2*eps
      i0=0.0d0
      if(ee.lt.20.0d0)then
       i0=pi3*GAB*GCD/dsqrt((ksi*ni)**3)*dsqrt(w**3)*exp(-ee)
      endif 

      ain=0.0d0
      i0n=0.0d0

c     r1
      x1=xmin-px/2.0d0
      do 1 ix=1,Nx
      call xxx(x1,px,x1a,x1b,ax,bx,ee,a,b)
      if(ee.lt.20.0d0)then
       g1=exp(-ee)
       e1=x1a**nax*x1b**nbx*g1
       y1=ymin-py/2.0d0
       do 2 iy=1,Ny
       call xxx(y1,py,y1a,y1b,ay,by,ee,a,b)
       if(ee.lt.20.0d0)then
        g2=exp(-ee)
        e2=y1a**nay*y1b**nby*g2
        z1=zmin-pz/2.0d0
        do 3 iz=1,Nz
        call xxx(z1,pz,z1a,z1b,az,bz,ee,a,b)
        if(ee.lt.20.0d0)then
         g3=exp(-ee)
         e3=z1a**naz*z1b**nbz*g3

c        r2
         x2=xmin-px/2.0d0
         do 4 jx=1,Nx
         call xxx(x2,px,x2c,x2d,cx,dx,ee,c,d)
         if(ee.lt.20.0d0)then
          g4=exp(-ee)
          e4=x2c**ncx*x2d**ndx*g4
          x12=x1-x2
          y2=ymin-py/2.0d0
          do 5 jy=1,Ny
          call xxx(y2,py,y2c,y2d,cy,dy,ee,c,d)
          if(ee.lt.20.0d0)then
           ee=exp(-ee)
           g5=g1*g2*g3*g4*ee
           e5=y2c**ncy*y2d**ndy*e1*e2*e3*e4*ee
           y12=y1-y2
           t2=x12**2+y12**2
           z2=zmin-pz/2.0d0
           do 6 jz=1,Nz
           call xxx(z2,pz,z2c,z2d,cz,dz,ee,c,d)
           ef=ee+eps*(t2+(z1-z2)**2)
           if(ef.lt.20.0d0)then
            ef=exp(-ef)
            i0n=i0n+                   g5*ef
            ain =ain+z2c**ncz*z2d**ndz*e5*ef
           endif
6          continue
          endif
5         continue
         endif
4        continue
        endif
3       continue
       endif
2      continue
      endif
1     continue

      ain=ain*(px*py*pz)**2
      i0n=i0n*(px*py*pz)**2

      ai=i0*
     1abcdg(nax,nay,naz,nbx,nby,nbz,ncx,ncy,ncz,ndx,ndy,ndz,
     1ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ksi,ni,we,
     1P(1),P(2),P(3),Q(1),Q(2),Q(3))

      write(6,719)i0n
719   format(' [0] numerical  :',g15.6)
      write(6,709)i0
709   format(' [0] analytical :',g15.6)
      write(6,609)ain
609   format(' I   numerical  :',g15.6)
      write(6,619)ai
619   format(' I   analytical :',g15.6)

      end

      subroutine mmx(xmin,xmax,ax,a,bx,b,
     1cx,c,dx,d,px,nx)
      implicit none
      integer*4 nx
      real*8 xmin,xmax,ax,a,bx,b,dx,d,cx,c,px,t
      t=5.0d0
      xmin=ax-t/dsqrt(a)
      if(bx-t/dsqrt(b).lt.xmin)xmin=bx-t/dsqrt(b)
      if(cx-t/dsqrt(c).lt.xmin)xmin=cx-t/dsqrt(c)
      if(dx-t/dsqrt(d).lt.xmin)xmin=dx-t/dsqrt(d)
      xmax=ax+t/dsqrt(a)
      if(bx+t/dsqrt(b).gt.xmax)xmax=bx+t/dsqrt(b)
      if(cx+t/dsqrt(c).gt.xmax)xmax=cx+t/dsqrt(c)
      if(dx+t/dsqrt(d).gt.xmax)xmax=dx+t/dsqrt(d)
      px=(xmax-xmin)/dble(Nx)
      return
      end

      subroutine xxx(x1,dx,x1a,x1b,ax,bx,eex,a,b)
      implicit none
      real*8 x1,dx,x1a,x1b,ax,bx,eex,a,b
      x1=x1+dx
      x1a=x1-ax
      x1b=x1-bx
      eex=a*x1a**2+b*x1b**2
      return
      end

      function abcdg(nxi,nyi,nzi,nxj,nyj,nzj,nxk,nyk,nzk,nxl,nyl,nzl,
     1ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ksip,ksiq,w,
     1px,py,pz,qx,qy,qz)
c
c     The gaussian integral (ij||kl) (AO) by a reccurence
c     relationship - make several primitive integrals at once
c     [0] not included
c
      implicit none
      integer*4 ns0
      parameter (ns0=1000)
      integer*4 nxi,nyi,nzi,nxj,nyj,nzj,nxk,nyk,nzk,nxl,nyl,nzl,np,nq,
     1abp(12),cdq(12),qq,alp(9,ns0),alq(9,ns0),alpq(9,ns0),
     1i,j,npq
      real*8 ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,abxyz(6),cdxyz(6),
     1cp(ns0),cq(ns0),ksip,ksiq,ssum,w,wr,hermite,Rx,Ry,Rz,wx,wy,wz,
     1abcdg,cpq(ns0),px,py,pz,qx,qy,qz,P(3),Q(3)
c
      Rx=qx-px
      Ry=qy-py
      Rz=qz-pz
      P(1)=px
      P(2)=py
      P(3)=pz
      Q(1)=qx
      Q(2)=qy
      Q(3)=qz
      abxyz( 1)=ax
      abxyz( 2)=ay
      abxyz( 3)=az
      abxyz( 4)=bx
      abxyz( 5)=by
      abxyz( 6)=bz
      cdxyz( 1)=cx
      cdxyz( 2)=cy
      cdxyz( 3)=cz
      cdxyz( 4)=dx
      cdxyz( 5)=dy
      cdxyz( 6)=dz
c
c     initiate (a b 0|:
      np=0
      abp(1)=nxi
      abp(2)=nyi
      abp(3)=nzi
      abp(4)=nxj
      abp(5)=nyj
      abp(6)=nzj
      abp(7)=0
      abp(8)=0
      abp(9)=0
c     just initialize so compiler does not complain:
      cp(1)=1.0d0
      alp(1,1)=0
      call addp(abp,np,1.0d0,cp,alp,ns0)


c     reduce to sum (0 0 p|:
      call reduce(np,cp,alp,ns0,ksip,abxyz,P)

      write(6,606)nxi,nyi,nzi,nxj,nyj,nzj,0,0,0
      do i=1,np
       write(6,605)cp(i),(alp(j,i),j=1,9)
605    format(g15.6,' ',9i1)
      enddo


      
c     initiate [0 d c): 
      nq=0
      cdq(1)=nxk
      cdq(2)=nyk
      cdq(3)=nzk
      cdq(4)=nxl
      cdq(5)=nyl
      cdq(6)=nzl
      cdq(7)=0
      cdq(8)=0
      cdq(9)=0
      cq(1)=1.0d0
      alq(1,1)=0
      call addp(cdq,nq,1.0d0,cq,alq,ns0)

c     reduce to sum  [q 0 0): 
      call reduce(nq,cq,alq,ns0,ksiq,cdxyz,Q)
      write(6,606)nxk,nyk,nzk,nxl,nyl,nzl,0,0,0
606   format(9i1,' reduced:')
      do i=1,nq
       write(6,605)cq(i),(alq(j,i),j=1,9)
      enddo

c     sum(i,j) cpi cqj [p|q] = sum(i) cpqj [q]
c     [p|q] = (-1)^p [p+q]
      npq=0
      cpq(1)=0.0d0
      alpq(1,1)=0
      do 6 i=1,np
      do 6 j=1,nq
6     call addpq(ns0,npq,cpq,i,j,cp,cq,alp,alq,alpq)

c     get the integral
c     [q]=[0] w^[(qx+qy+qz)/2] (-1)^(qx+qy+qz) Hqx(wx)Hqy(wy)H(qz(wz)
c     [0] not included
      ssum=0.0d0
      wr=sqrt(w)
      wx=wr*Rx
      wy=wr*Ry
      wz=wr*Rz
      do 7 i=1,npq
      qq=alpq(1,i)+alpq(2,i)+alpq(3,i)
      write(6,*)i,cpq(i)*dsqrt(w**qq)*(-1.0d0)**qq
     1*hermite(alpq(1,i),wx)*hermite(alpq(2,i),wy)*hermite(alpq(3,i),wz)
      write(6,*)w,qq
      write(6,*)cpq(i),dsqrt(w**qq),(-1.0d0)**qq
      write(6,*)alpq(1,i),wx,hermite(alpq(1,i),wx)
      write(6,*)alpq(2,i),wy,hermite(alpq(2,i),wy)
      write(6,*)alpq(3,i),wz,hermite(alpq(3,i),wz)
7     ssum=ssum+cpq(i)*dsqrt(w**qq)*(-1.0d0)**qq
     1*hermite(alpq(1,i),wx)*hermite(alpq(2,i),wy)*hermite(alpq(3,i),wz)

      abcdg=ssum

      return
      end

      subroutine reduce(np,cp,alp,ns0,ksi,axyz,pxyz)
c     reduce [ab to sum cp [p]
      implicit none
      integer*4 np,alp(9,ns0),ns0,i,j,k,ix,abp(9),px
      real*8 cp(*),axyz(6),pxyz(3),ksi,co


1000  do 4 i=1,np
      do 4 j=1,6
      if(alp(j,i).gt.0)then
c      remember this:
       co=cp(i)
       do 6 k=1,9
6      abp(k)=alp(k,i)
c      delete the original term:
       call dddp(ns0,i,np,cp,alp)
       write(6,*)i,'deleted,',np,' remains'
c      and add three new:
c      (ax b p ] = px (ax-1 b px-1]+(Px-Ax) (ax-1 b px] +1/(2ksi) (ax-1 b px+1]
       abp(j)=abp(j)-1
       if(j.le.3)then
        ix=j
       else
        ix=j-3
       endif
       px=abp(6+ix)
       abp(6+ix)=abp(6+ix)-1
       if(px.gt.0)call addp(abp,np,co*dble(px),cp,alp,ns0)
       write(6,*)'px,',np,' remains'
       abp(6+ix)=abp(6+ix)+1
       call addp(abp,np,co*(pxyz(ix)-axyz(j)),cp,alp,ns0)
       write(6,*)'pa,',np,' remains',co,ix,j,pxyz(ix),axyz(j)
       abp(6+ix)=abp(6+ix)+1
       call addp(abp,np,co*0.5d0/ksi,cp,alp,ns0)
       write(6,*)'ksi,',np,' remains'
c      take a fresh look:
       goto 1000
      endif
4     continue
      return
      end


      subroutine addp(abp,n,am,fi,al,ns0)
c     add new term (a b p| to those in the list
      implicit none
      integer*4 abp(*),n,i,j,ns0,al(9,ns0)
      real*8 am,fi(*)
c
c     is it already in?:
      do 1 i=1,n
      do 2 j=1,9
2     if(al(j,i).ne.abp(j))goto 1
c     the term is in term i, add it there:
      fi(i)=fi(i)+am
c     accidently zero?:
      if(dabs(fi(i)).lt.1.0d-12) call  dddp(ns0,i,n,fi,al)
      return
1     continue

c     not in, add new:
      n=n+1
      if(n.gt.ns0)call report('addp expansion overflow')
      fi(n)=am
      do 3 j=1,9
3     al(j,n)=abp(j)
c     accidently zero?:
      if(dabs(fi(n)).lt.1.0d-12) call  dddp(ns0,n,n,fi,al)
      return
      end


      subroutine addpq(ns0,npq,cpq,i,j,cp,cq,alp,alq,alpq)
      implicit none
      integer*4 ns0,npq,i,j,alp(9,ns0),alq(9,ns0),alpq(9,ns0),iexp,
     1alt(3),it,ix
      real*8 cpq(*),cp(*),cq(*),ssi
c
c     [p|q] = (-1)^p [p+q]:
c     px+py+pz:
      iexp=alp(6+1,i)+alp(6+2,i)+alp(6+3,i)
c     [p+q]:
      do 1 ix=1,3
1     alt(ix)=alq(6+ix,j)+alp(6+ix,i)
      if(mod(iexp,2).eq.0)then
       ssi=1.0d0
      else
       ssi=-1.0d0
      endif
c     see if this term exists:
      do 2 it=1,npq
      if(alpq(1,it).eq.alt(1).and.
     1   alpq(2,it).eq.alt(2).and.alpq(3,it).eq.alt(3))then
       cpq(it)=cpq(it)+ssi*cp(i)*cq(i)
       return
      endif
2     continue
      npq=npq+1
      if(npq.gt.ns0)call report('npq overflow in addpq')
      cpq(npq)=ssi*cp(i)*cq(j)
      alpq(1,npq)=alt(1)
      alpq(2,npq)=alt(2)
      alpq(3,npq)=alt(3)
      return
      end

      function hermite(n,x)
      implicit none
      integer*4 n,i
      real*8 hermite,tx,hm,hmm,hi,x
      hmm=1.0d0
      if(n.eq.0)then
       hermite=hmm
       return
      endif
      hm=x+x
      if(n.eq.1)then
       hermite=hm
       return
      endif
      tx=hm
      do 1 i=2,n
c     Hi = 2x Hi-1 - 2(i-1) Hi-2
      hi=tx*hm-dble(i+i-2)*hmm
      hmm=hm
1     hm=hi
      hermite=hi
      return
      end

      subroutine dddp(ns0,i,n,fi,al)
c     delete ith term (a b p| from the expansion
      implicit none
      integer*4 ns0
      integer*4 n,i,al(9,ns0)
      real*8 fi(*)
      if(i.lt.n)then
       fi(  i)=fi(  n)
       al(1,i)=al(1,n)
       al(2,i)=al(2,n)
       al(3,i)=al(3,n)
       al(4,i)=al(4,n)
       al(5,i)=al(5,n)
       al(6,i)=al(6,n)
       al(7,i)=al(7,n)
       al(8,i)=al(8,n)
       al(9,i)=al(9,n)
      endif
      n=n-1
      return
      end

      SUBROUTINE REPORT(STRING)
      CHARACTER*(*) STRING
      WRITE(6,1000)
1000  FORMAT(80(1H*))
      WRITE(6,*)STRING
      WRITE(6,1000)
      WRITE(6,2000)
2000  FORMAT(' Program terminated')
      STOP
      END