diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f
index 83616969..1193e126 100644
--- a/src/utils/cgtos_one_e.irp.f
+++ b/src/utils/cgtos_one_e.irp.f
@@ -21,7 +21,7 @@ complex*16 function overlap_cgaussian_x(Ae_center, Be_center, alpha, beta, power
   integer                :: i, iorder_p
   complex*16             :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p
 
-  complex*16             :: Fc_integral
+  complex*16, external   :: Fc_integral
 
 
   call give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_p, iorder_p, &
diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f
index 075d9bbc..8d18ef46 100644
--- a/src/utils/cgtos_utils.irp.f
+++ b/src/utils/cgtos_utils.irp.f
@@ -421,11 +421,12 @@ complex*16 function Fc_integral(n, inv_sq_p)
   implicit none
   include 'constants.include.F'
 
-  integer,    intent(in) :: n
-  complex*16, intent(in) :: inv_sq_p 
+  integer,    intent(in)     :: n
+  complex*16, intent(in)     :: inv_sq_p 
 
+  complex*16                 :: inv_sq_p2, inv_sq_p3, inv_sq_p4
   ! (n)! 
-  double precision       :: fact
+  double precision, external :: fact
 
   if(n < 0) then
     Fc_integral = (0.d0, 0.d0)
@@ -438,13 +439,29 @@ complex*16 function Fc_integral(n, inv_sq_p)
     return
   endif
 
-  if(n == 0) then
+  select case(n)
+  case(0)
     Fc_integral = sqpi * inv_sq_p
-    return
-  endif
-
-  Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1))
+  case(2)
+    Fc_integral = 0.5d0 * sqpi * inv_sq_p * inv_sq_p * inv_sq_p
+  case(4)
+    inv_sq_p2 = inv_sq_p * inv_sq_p
+    Fc_integral = 0.75d0 * sqpi * inv_sq_p * inv_sq_p2 * inv_sq_p2
+  case(6)
+    inv_sq_p3 = inv_sq_p * inv_sq_p * inv_sq_p
+    Fc_integral = 1.875d0 * sqpi * inv_sq_p * inv_sq_p3 * inv_sq_p3
+  case(8)
+    inv_sq_p3 = inv_sq_p * inv_sq_p * inv_sq_p
+    Fc_integral = 6.5625d0 * sqpi * inv_sq_p3 * inv_sq_p3 * inv_sq_p3
+  case(10)
+    inv_sq_p2 = inv_sq_p * inv_sq_p
+    inv_sq_p4 = inv_sq_p2 * inv_sq_p2
+    Fc_integral = 29.53125d0 * sqpi * inv_sq_p * inv_sq_p2 * inv_sq_p4 * inv_sq_p4
+  case default
+    Fc_integral = sqpi * 0.5d0**n * inv_sq_p**(n+1) * fact(n) / fact(shiftr(n, 1))
+  end select
 
+  return
 end
 
 ! ---