From 41e69c56e50ff68ec1d3b2e0b1ec0adb6c2de3d7 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 30 Sep 2016 15:34:59 +0200
Subject: [PATCH] HF_bitmask in cas definition

---
 src/Determinants/psi_cas.irp.f       |  4 ++--
 src/Integrals_Monoelec/pseudopot.f90 | 13 +++++-------
 src/Utils/util.irp.f                 | 31 ++++++++++++++++++++--------
 3 files changed, 29 insertions(+), 19 deletions(-)

diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f
index 304a2370..968ced57 100644
--- a/src/Determinants/psi_cas.irp.f
+++ b/src/Determinants/psi_cas.irp.f
@@ -21,9 +21,9 @@ use bitmasks
       do k=1,N_int
         good = good .and. (                                          &
             iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
-            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,1,l)), hf_bitmask(k,1)) ) .and. (  &
             iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
-            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+            iand(not(cas_bitmask(k,2,l)), hf_bitmask(k,2)) )
       enddo
       if (good) then
         exit
diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90
index c89bb019..a21ee836 100644
--- a/src/Integrals_Monoelec/pseudopot.f90
+++ b/src/Integrals_Monoelec/pseudopot.f90
@@ -136,7 +136,7 @@ end
 
          if(l.eq.2.and.m.eq.2)ylm_real=dsqrt(15.d0/16.d0/pi)*(xchap*xchap-ychap*ychap)
          if(l.eq.2.and.m.eq.1)ylm_real=dsqrt(15.d0/fourpi)*xchap*zchap
-         if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(-xchap*xchap-ychap*ychap+2.d0*zchap*zchap)
+         if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(2.d0*zchap*zchap-xchap*xchap-ychap*ychap)
          if(l.eq.2.and.m.eq.-1)ylm_real=dsqrt(15.d0/fourpi)*ychap*zchap
          if(l.eq.2.and.m.eq.-2)ylm_real=dsqrt(15.d0/fourpi)*xchap*ychap
 
@@ -279,8 +279,7 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
    do m=-l,l
     prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
     prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
-
-      accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
+    accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
 
    enddo
   enddo
@@ -1783,10 +1782,8 @@ double precision function coef_nk(n,k)
   double precision gam,dble_fact,fact
 
   gam=dble_fact(n+n+k+k+1)
-
-!  coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
-  
-  coef_nk=1.d0/(2.d0**k*fact(k)*gam)
+!  coef_nk=1.d0/(2.d0**k*fact(k)*gam)
+  coef_nk=1.d0/(dble(ibset(0,k))*fact(k)*gam)
   
   return
 
@@ -1815,7 +1812,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
   logical done
 
-  u=(a+b)/(2.d0*dsqrt(gam))
+  u=(a+b)*0.5d0/dsqrt(gam)
   freal=dexp(-arg)
 
   if(a.eq.0.d0.and.b.eq.0.d0)then
diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f
index 751610a5..8d375ac6 100644
--- a/src/Utils/util.irp.f
+++ b/src/Utils/util.irp.f
@@ -84,10 +84,8 @@ double precision function fact(n)
     memo(i) = memo(i-1)*dble(i)
   enddo
   memomax = min(n,100)
-  fact = memo(memomax)
-  do i=101,n
-    fact = fact*dble(i)
-  enddo
+  double precision :: logfact
+  fact = dexp(logfact(n))
 end function
 
 double precision function logfact(n)
@@ -164,12 +162,27 @@ double precision function dble_fact_even(n) result(fact2)
 
   ASSERT (iand(n,1) /= 1)
 
-  prod=1.d0
-  do k=2,n,2
-   prod=prod*dfloat(k)
+  if (n <= memomax) then
+    if (n < 2) then
+      fact2 = 1.d0
+    else
+      fact2 = memo(n)
+    endif
+    return
+  endif
+
+  integer                        :: i
+  memo(1)=1.d0
+  do i=memomax+2,min(n,100),2
+    memo(i) = memo(i-2)* dble(i)
   enddo
-  fact2=prod
-  return
+  memomax = min(n,100)
+  fact2 = memo(memomax)
+  
+  if (n > 100) then
+    double precision :: dble_logfact
+    fact2 = dexp(dble_logfact(n))
+  endif
 
 end function