program inttwoel
c     two-electron integral
      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,i
      real*8 ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ai,aig
      real*8 a,b,c,d,rr,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,r12,
     1e1,e2,e3,f1,f2,px,py,pz,ksi,ni,R(3),P(3),q(3),
     1abx,aby,abz,cdx,cdy,cdz,GAB,GCD,i0,cw(12),ew(12),pi,pi3,
     1w,r2,i0n,g1,g2,g3,h1,h2,th2,t2,it(12),itn(12),ef,xx

      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,*)rr
      read(9,*)nx,ny,nz
      close(9)

      open(8,file='R.FIT')
      read(8,*)
      do 10 i=1,12
10    read(8,*)cw(i),cw(i),ew(i)
      close(8)



      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)
      i0=0.0d0
      do 7 i=1,12
      w=th2/(th2+ew(i))
      ee=w*r2*ew(i)
      it(i)=0.0d0
      itn(i)=0.0d0
      if(ee.lt.20.0d0)then
       it(i)=pi3*GAB*GCD/dsqrt((ksi*ni)**3)*dsqrt(w**3)*exp(-ee)
       i0=i0+cw(i)*it(i)
      endif 
7     continue
 


      ai=0.0d0
      aig=0.0d0
      i0n=0.0d0

c     r1
      x1=xmin-px/2.0d0
      do 1 ix=1,Nx
      write(6,*)ix
      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
          h1=exp(-ee)
          f1=x2c**ncx*x2d**ndx*h1
          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)
           h2=g1*g2*g3*h1*ee
           f2=y2c**ncy*y2d**ndy*e1*e2*e3*f1*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)
           if(ee.lt.20.0d0)then
            r12=dsqrt(t2+(z1-z2)**2)
            ee=exp(-ee)
            ai=ai+z2c**ncz*z2d**ndz*ee*f2/(r12+rr)**6
            i0n=i0n+ee*h2/(r12+rr)**6

            do 8 i=1,12
            ef=ew(i)*r12**2
            if(ef.lt.100.0d0)then
             itn(i)=itn(i)+ee*h2*exp(-ef)
             aig=aig+cw(i)*ee*h2*exp(-ef)
            endif
8           continue

           endif
6          continue
          endif
5         continue
         endif
4        continue
        endif
3       continue
       endif
2      continue
      endif
1     continue


      ai =ai *(px*py*pz)**2
      aig=aig*(px*py*pz)**2
      i0n=i0n*(px*py*pz)**2
      do 9 i=1,12
      itn(i)=itn(i)*(px*py*pz)**2
9     write(6,603)i,it(i),itn(i)
603   format(i3,2f18.9)
      

      write(6,709)i0
709   format(' [0] analytical :',f12.6)
      write(6,719)i0n
719   format(' [0] numerical  :',f12.6)
      write(6,609)ai
609   format(' I              :',F22.6)
      write(6,109)aig
109   format(' I from G       :',F22.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
      xmin=ax-3.0d0/dsqrt(a)
      if(bx-3.0d0/dsqrt(b).lt.xmin)xmin=bx-3.0d0/dsqrt(b)
      if(cx-3.0d0/dsqrt(c).lt.xmin)xmin=cx-3.0d0/dsqrt(c)
      if(dx-3.0d0/dsqrt(d).lt.xmin)xmin=dx-3.0d0/dsqrt(d)
      xmax=ax+3.0d0/dsqrt(a)
      if(bx+3.0d0/dsqrt(b).gt.xmax)xmax=bx+3.0d0/dsqrt(b)
      if(cx+3.0d0/dsqrt(c).gt.xmax)xmax=cx+3.0d0/dsqrt(c)
      if(dx+3.0d0/dsqrt(d).gt.xmax)xmax=dx+3.0d0/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