From cca5ebc404ab3c6689a677b895afbe0b0ee80329 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 7 Apr 2015 14:24:05 +0200
Subject: [PATCH] Move need.irp.f into MonoInts

---
 src/MonoInts/need.irp.f | 289 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 289 insertions(+)
 create mode 100644 src/MonoInts/need.irp.f

diff --git a/src/MonoInts/need.irp.f b/src/MonoInts/need.irp.f
new file mode 100644
index 00000000..22cb6a48
--- /dev/null
+++ b/src/MonoInts/need.irp.f
@@ -0,0 +1,289 @@
+ 
+      double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB) 
+      implicit double precision(a-h,o-z)
+      dimension nA(3),nB(3)
+      dimension A(3),B(3)
+      gamtot=gamA+gamB
+      SABpartial=1.d0
+
+      l=3
+      u=gamA/gamtot*A(l)+gamB/gamtot*B(l)
+      arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2
+      alpha=dexp(arg)
+     &/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0)
+      wA=dsqrt(gamtot)*(u-A(l))
+      wB=dsqrt(gamtot)*(u-B(l))
+      boundA=dsqrt(gamtot)*(zA-u)
+      boundB=dsqrt(gamtot)*(zB-u)
+
+      accu=0.d0
+      do n=0,nA(l)
+       do m=0,nB(l)
+        integ=nA(l)+nB(l)-n-m
+        accu=accu
+     & +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m)
+     & *(rinteg(integ,boundB)-rinteg(integ,boundA))
+       enddo
+      enddo
+      SABpartial=SABpartial*accu*alpha
+      end
+
+      double precision function rintgauss(n)
+      implicit double precision(a-h,o-z)
+      rintgauss=dsqrt(dacos(-1.d0))
+      if(n.eq.0)return
+      if(n.eq.1)then
+       rintgauss=0.d0
+       return
+      endif
+      if(iand(n,1).eq.1)then
+       rintgauss=0.d0
+       return
+      endif
+      rintgauss=rintgauss/2.d0**(n/2)
+      rintgauss=rintgauss*ddfact2(n-1)
+      end
+
+      double precision function rinteg(n,u)
+      implicit double precision(a-h,o-z)
+      include 'constants.F'
+!     pi=dacos(-1.d0)
+      ichange=1
+      factor=1.d0
+      if(u.lt.0.d0)then
+       u=-u
+       factor=(-1.d0)**(n+1)
+       ichange=-1
+      endif
+      if(iand(n,1).eq.0)then
+       rinteg=0.d0
+       do l=0,n-2,2
+        prod=b_coef(l,u)
+        do k=l+2,n-2,2
+         prod=prod*a_coef(k)
+        enddo
+        rinteg=rinteg+prod
+       enddo
+       prod=dsqrt(pi)/2.d0*erf0(u)
+       do k=0,n-2,2
+        prod=prod*a_coef(k)
+       enddo
+       rinteg=rinteg+prod
+      endif
+
+      if(iand(n,1).eq.1)then
+       rinteg=0.d0
+       do l=1,n-2,2
+        prod=b_coef(l,u)
+        do k=l+2,n-2,2
+         prod=prod*a_coef(k)
+        enddo
+        rinteg=rinteg+prod
+       enddo
+       prod=0.5d0*(1.d0-dexp(-u**2))
+       do k=1,n-2,2
+        prod=prod*a_coef(k)
+       enddo
+       rinteg=rinteg+prod
+      endif
+
+      rinteg=rinteg*factor
+
+      if(ichange.eq.-1)u=-u
+
+      end
+
+!<function type="double precision function" name="erf0">
+!  <arg name="x"
+!       doc ="" />
+!
+!  <doc>
+!
+!  </doc>
+!
+! <fortran>
+      double precision function erf0(x)
+      implicit double precision (a-h,o-z)
+      if(x.lt.0.d0)then
+        erf0=-gammp(0.5d0,x**2)
+      else
+        erf0=gammp(0.5d0,x**2)
+      endif
+      end
+
+
+! </fortran>
+!</function>
+!<function type="double precision function" name="gammp">
+!  <arg name="a"
+!       doc ="" />
+!  <arg name="x"
+!       doc ="" />
+!
+!  <doc>
+!
+!  </doc>
+!
+!  <calls>
+!    gcf
+!    gser
+!  </calls>
+!
+! <fortran>
+      double precision function gammp(a,x)
+      implicit double precision (a-h,o-z)
+      if(x.lt.0..or.a.le.0.)stop 'error in gammp'
+      if(x.lt.a+1.)then
+        call gser(gammp,a,x,gln)
+      else
+        call gcf(gammcf,a,x,gln)
+        gammp=1.-gammcf
+      endif
+      return
+      end
+! </fortran>
+!</function>
+
+
+!<function type="subroutine" name="gser">
+!  <arg name="gamser"
+!       doc ="" />
+!  <arg name="a"
+!       doc ="" />
+!  <arg name="x"
+!       doc ="" />
+!  <arg name="gln"
+!       doc ="" />
+!
+!  <doc>
+!
+!  </doc>
+!
+!  <calledBy>
+!    gammp
+!  </calledBy>
+!
+! <fortran>
+      subroutine gser(gamser,a,x,gln)
+      implicit double precision (a-h,o-z)
+      parameter (itmax=100,eps=3.e-7)
+      gln=gammln(a)
+      if(x.le.0.)then
+        if(x.lt.0.) stop 'error in gser'
+        gamser=0.
+        return
+      endif
+      ap=a
+      sum=1./a
+      del=sum
+      do 11 n=1,itmax
+        ap=ap+1.
+        del=del*x/ap
+        sum=sum+del
+        if(abs(del).lt.abs(sum)*eps)go to 1
+11    continue
+      stop 'a too large, itmax too small'
+1     gamser=sum*exp(-x+a*log(x)-gln)
+      return
+      end
+! </fortran>
+
+!</function>
+!<function type="subroutine" name="gcf">
+!  <arg name="gammcf"
+!       doc ="" />
+!  <arg name="a"
+!       doc ="" />
+!  <arg name="x"
+!       doc ="" />
+!  <arg name="gln"
+!       doc ="" />
+!
+!  <doc>
+!
+!  </doc>
+!
+!  <calledBy>
+!    gammp
+!  </calledBy>
+!
+! <fortran>
+      subroutine gcf(gammcf,a,x,gln)
+      implicit double precision (a-h,o-z)
+      parameter (itmax=100,eps=3.e-7)
+      gln=gammln(a)
+      gold=0.
+      a0=1.
+      a1=x
+      b0=0.
+      b1=1.
+      fac=1.
+      do 11 n=1,itmax
+        an=float(n)
+        ana=an-a
+        a0=(a1+a0*ana)*fac
+        b0=(b1+b0*ana)*fac
+        anf=an*fac
+        a1=x*a0+anf*a1
+        b1=x*b0+anf*b1
+        if(a1.ne.0.)then
+          fac=1./a1
+          g=b1*fac
+          if(abs((g-gold)/g).lt.eps)go to 1
+          gold=g
+        endif
+11    continue
+      stop 'a too large, itmax too small'
+1     gammcf=exp(-x+a*log(x)-gln)*g
+      return
+      end
+
+! </fortran>
+!</function>
+      double precision function ddfact2(n)
+      implicit double precision(a-h,o-z)
+      if(iand(n,1).eq.0)stop 'error in ddfact2'
+      ddfact2=1.d0
+      do i=1,n,2
+       ddfact2=ddfact2*dfloat(i)
+      enddo
+      end
+
+      double precision function a_coef(n)
+      implicit double precision(a-h,o-z)
+      a_coef=dfloat(n+1)/2.d0
+      end
+
+      double precision function b_coef(n,u)
+      implicit double precision(a-h,o-z)
+      b_coef=-0.5d0*u**(n+1)*dexp(-u**2)
+      end
+
+!<function type="double precision function" name="gammln">
+!  <arg name="xx"
+!       doc ="" />
+!
+!  <doc>
+!
+!  </doc>
+!
+! <fortran>
+      double precision function gammln(xx)
+      implicit double precision (a-h,o-z)
+      real*8 cof(6),stp,half,one,fpf,x,tmp,ser
+      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
+     *    -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
+      data half,one,fpf/0.5d0,1.0d0,5.5d0/
+      x=xx-one
+      tmp=x+fpf
+      tmp=(x+half)*log(tmp)-tmp
+      ser=one
+      do 11 j=1,6
+        x=x+one
+        ser=ser+cof(j)/x
+11    continue
+      gammln=tmp+log(stp*ser)
+      return
+      end
+! </fortran>
+!</function>