From dc20ea7c07f9a198927800decb26c085f4dfd39a Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Mon, 4 May 2015 20:26:15 +0200
Subject: [PATCH 01/18] ezfio.pseudo_intergrals => ezfio.pseudo_integrals

---
 scripts/pseudo/put_pseudo_in_ezfio.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py
index 87db7845..1f8594c5 100755
--- a/scripts/pseudo/put_pseudo_in_ezfio.py
+++ b/scripts/pseudo/put_pseudo_in_ezfio.py
@@ -311,11 +311,11 @@ if __name__ == "__main__":
     # ~#~#~#~#~ #
 
     klocmax = max([len(i) for i in v_k])
-    ezfio.pseudo_intergrals_klocmax = klocmax
+    ezfio.pseudo_integrals_klocmax = klocmax
 
-    ezfio.pseudo_intergrals_v_k = zip(*v_k)
-    ezfio.pseudo_intergrals_n_k = zip(*n_k)
-    ezfio.pseudo_intergrals_dz_k = zip(*dz_k)
+    ezfio.pseudo_integrals_v_k = zip(*v_k)
+    ezfio.pseudo_integrals_n_k = zip(*n_k)
+    ezfio.pseudo_integrals_dz_k = zip(*dz_k)
 
     # ~#~#~#~#~#~#~#~#~ #
     # N o n _ L o c a l #
@@ -324,15 +324,15 @@ if __name__ == "__main__":
     lmax = max([len(i) for i in v_kl])
     kmax = max([len(sublist) for list_ in v_kl for sublist in list_])
 
-    ezfio.pseudo_intergrals_lmaxpo = lmax
-    ezfio.pseudo_intergrals_kmax = kmax
+    ezfio.pseudo_integrals_lmaxpo = lmax
+    ezfio.pseudo_integrals_kmax = kmax
 
     v_kl = make_it_square(v_kl, [lmax, kmax])
     n_kl = make_it_square(n_kl, [lmax, kmax], int)
     dz_kl = make_it_square(dz_kl, [lmax, kmax])
 
-    ezfio.pseudo_intergrals_v_kl = zip(*v_kl)
-    ezfio.pseudo_intergrals_n_kl = zip(*n_kl)
-    ezfio.pseudo_intergrals_dz_kl = zip(*dz_kl)
+    ezfio.pseudo_integrals_v_kl = zip(*v_kl)
+    ezfio.pseudo_integrals_n_kl = zip(*n_kl)
+    ezfio.pseudo_integrals_dz_kl = zip(*dz_kl)
 
-    ezfio.pseudo_intergrals_do_pseudo = True
+    ezfio.pseudo_integrals_do_pseudo = True

From ad3d6317b230ef1b287b2100431af38ece9737e3 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 10:28:04 +0200
Subject: [PATCH 02/18] Add Pseudo Unit Test for SO2

---
 src/FCIdump/README.rst       |  3 ++
 tests/SO2.xyz                |  5 +++
 tests/unit_test/unit_test.py | 86 +++++++++++++++++++++++++++---------
 3 files changed, 72 insertions(+), 22 deletions(-)
 create mode 100644 tests/SO2.xyz

diff --git a/src/FCIdump/README.rst b/src/FCIdump/README.rst
index faf69203..8a467b16 100644
--- a/src/FCIdump/README.rst
+++ b/src/FCIdump/README.rst
@@ -10,6 +10,9 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`fcidump <http://github.com/LCPQ/quantum_package/tree/master/src/FCIdump/fcidump.irp.f#L1>`_
+  Undocumented
+
 
 
 Needed Modules
diff --git a/tests/SO2.xyz b/tests/SO2.xyz
new file mode 100644
index 00000000..5b6fca03
--- /dev/null
+++ b/tests/SO2.xyz
@@ -0,0 +1,5 @@
+3
+SO2 Geo: Experiment Mult: 1 symmetry: 32
+O    0.0    1.2371    0.7215
+O    0.0   -1.2371    0.7215
+S    0.0    0.0       0.0
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index d9498700..19d129d1 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -12,6 +12,9 @@ sys.path = [EZFIO + "/Python"] + sys.path
 
 from ezfio import ezfio
 from collections import defaultdict
+from collections import namedtuple
+
+Energy = namedtuple('Energy', ['without_pseudo', 'with_pseudo'])
 
 # ~#~#~ #
 # O p t #
@@ -34,7 +37,7 @@ has_hf_alredy = False
 global filename_check
 
 
-def init_folder(geo, basis, mult=1, ezfio_name=None):
+def init_folder(geo, basis, mult=1, pseudo=False, ezfio_name=None):
     '''
     Take a geo in arg (aka a existing geo.xyz in test/)
     And create the geo.ezfio with the adeguate basis and multipliciti
@@ -47,15 +50,21 @@ def init_folder(geo, basis, mult=1, ezfio_name=None):
     if not ezfio_name:
         ezfio_name = geo
 
-    cmd = "qp_create_ezfio_from_xyz -b {0} -m {1} {2}.xyz -o {3}.ezfio"
+    if pseudo:
+        cmd = "qp_create_ezfio_from_xyz -b {0} -m {1} {2}.xyz -p -o {3}.ezfio"
+    else:
+        cmd = "qp_create_ezfio_from_xyz -b {0} -m {1} {2}.xyz -o {3}.ezfio"
+
     subprocess.check_call([cmd.format(basis, mult, geo, ezfio_name)],
                           shell=True)
 
+    subprocess.call(["rm {0}.xyz".format(geo)], shell=True)
+
 
 def get_error_message(l_exepected, l_cur):
     l_msg = ["Need {0} get {1} error is {2}".format(i, j, abs(i - j))
              for i, j in zip(l_exepected, l_cur)]
-    return "\n".join(l_msg)
+    return "\n" + "\n".join(l_msg)
 
 
 #  _
@@ -146,18 +155,20 @@ def check_mo_guess(geo, basis, mult=1):
 # /  |_   _   _ |        _. |      _   _
 # \_ | | (/_ (_ |<   \/ (_| | |_| (/_ _>
 #
-def run_hf(geo, basis, mult=1):
+def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
     """
     Run a simle by default hf
     EZFIO path = geo.ezfio
     """
+
     # ~#~#~#~#~#~#~#~#~#~ #
     # R e f _ e n e r g y #
     # ~#~#~#~#~#~#~#~#~#~ #
 
-    ref_energy = defaultdict(dict)
+    ref_energy = defaultdict(defaultdict)
 
-    ref_energy["sto-3g"]["methane"] = -39.7267433402
+    ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
+    ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #
@@ -170,7 +181,7 @@ def run_hf(geo, basis, mult=1):
     # I n i t #
     # ~#~#~#~ #
 
-    init_folder(geo, basis, mult)
+    init_folder(geo, basis, mult, pseudo)
     ezfio.set_file("{0}.ezfio".format(geo))
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~ #
@@ -187,7 +198,7 @@ def run_hf(geo, basis, mult=1):
     ezfio.hartree_fock_thresh_scf = 1.e-10
     ezfio.hartree_fock_n_it_scf_max = 100
 
-    ezfio.pseudo_integrals_do_pseudo = False
+    ezfio.pseudo_integrals_do_pseudo = pseudo
 
     # ~#~#~ #
     # R u n #
@@ -201,15 +212,25 @@ def run_hf(geo, basis, mult=1):
     # ~#~#~#~#~ #
 
     cur_e = ezfio.get_hartree_fock_energy()
+
     ref_e = ref_energy[basis][geo]
+    if pseudo:
+        ref_e = ref_e.with_pseudo
+    else:
+        ref_e = ref_e.without_pseudo
 
     if abs(cur_e - ref_e) <= precision:
+
+        if remove_after_sucess:
+            subprocess.call(["rm -R {0}.ezfio".format(geo)], shell=True)
+
         return True
+
     else:
         raise ValueError(get_error_message([ref_e], [cur_e]))
 
 
-def run_full_ci_10k_pt2_end(geo, basis):
+def run_full_ci_10k_pt2_end(geo, basis, pseudo):
     """
     Run a Full_ci with 10k with the TruePT2
     EZFIO path = geo.ezfio
@@ -222,8 +243,8 @@ def run_full_ci_10k_pt2_end(geo, basis):
     ref_energy_var = defaultdict(dict)
     ref_energy_pt2 = defaultdict(dict)
 
-    ref_energy_var["sto-3g"]["methane"] = -0.398058753535695E+02
-    ref_energy_pt2["sto-3g"]["methane"] = -0.398059182483741E+02
+    ref_energy_var["sto-3g"]["methane"] = Energy(-0.398058753535695E+02, None)
+    ref_energy_pt2["sto-3g"]["methane"] = Energy(-0.398059182483741E+02, None)
 
     # ~#~#~#~ #
     # I n i t #
@@ -256,6 +277,13 @@ def run_full_ci_10k_pt2_end(geo, basis):
     ref_var = ref_energy_var[basis][geo]
     ref_pt2 = ref_energy_pt2[basis][geo]
 
+    if pseudo:
+        ref_var = ref_var.with_pseudo
+        ref_pt2 = ref_pt2.with_pseudo
+    else:
+        ref_var = ref_var.without_pseudo
+        ref_pt2 = ref_pt2.without_pseudo
+
     t = [abs(cur_var - ref_var) <= precision,
          abs(cur_pt2 - ref_pt2) <= precision]
 
@@ -266,15 +294,16 @@ def run_full_ci_10k_pt2_end(geo, basis):
                                            [cur_var, cur_pt2]))
 
 
-def hf_then_10k_test(geo, basis):
-    if not has_hf_alredy:
-        run_hf(geo, basis)
+def hf_then_10k_test(geo, basis, mult=1, pseudo=False):
+
+    run_hf(geo, basis, mult, pseudo, remove_after_sucess=False)
 
     try:
-        run_full_ci_10k_pt2_end(geo, basis)
-        return_code = True
+        run_full_ci_10k_pt2_end(geo, basis, pseudo)
     except:
-        return_code = False
+        raise
+    else:
+        return_code = True
 
     # ~#~#~#~#~#~#~#~ #
     # F i n a l i z e #
@@ -335,15 +364,27 @@ def check_convert(path_out):
     else:
         raise ValueError(get_error_message([ref_e], [cur_e]))
 
+
 # ___
 #  |  _   _ _|_
 #  | (/_ _>  |_
 #
 class ValueTest(unittest.TestCase):
 
-    def test_full_ci_10k_pt2_end(self):
-        self.assertTrue(hf_then_10k_test("methane", "sto-3g"))
+    def test_hf_then_full_ci_10k_pt2_end(self):
+        self.assertTrue(hf_then_10k_test(geo="methane",
+                                         basis="sto-3g",
+                                         mult=1,
+                                         pseudo=False))
 
+    def test_hf(self):
+        self.assertTrue(run_hf(geo="SO2",
+                               basis="vdz",
+                               mult=1,
+                               pseudo=True))
+
+
+class ConvertTest(unittest.TestCase):
     def test_check_convert_hf_energy(self):
         self.assertTrue(check_convert("HBO.out"))
 
@@ -351,11 +392,12 @@ class ValueTest(unittest.TestCase):
 class InputTest(unittest.TestCase):
 
     def test_check_disk_acess(self):
-        self.assertTrue(check_disk_acess("methane", "un-ccemd-ref"))
+        self.assertTrue(check_disk_acess(geo="methane",
+                                         basis="un-ccemd-ref"))
 
     def test_check_mo_guess(self):
-        self.assertTrue(check_mo_guess("methane", "maug-cc-pVDZ"))
-
+        self.assertTrue(check_mo_guess(geo="methane",
+                                       basis="maug-cc-pVDZ"))
 
 if __name__ == '__main__':
     unittest.main()

From ec26b457b305acd278a62731417e981dd6f3b33d Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 10:33:22 +0200
Subject: [PATCH 03/18] Change SO2 for HBO in unitest (test need to be fast)

---
 tests/HBO.xyz                | 5 +++++
 tests/unit_test/unit_test.py | 3 ++-
 2 files changed, 7 insertions(+), 1 deletion(-)
 create mode 100644 tests/HBO.xyz

diff --git a/tests/HBO.xyz b/tests/HBO.xyz
new file mode 100644
index 00000000..796a1e19
--- /dev/null
+++ b/tests/HBO.xyz
@@ -0,0 +1,5 @@
+3
+HBO Geo: Experiment Mult: 1 symmetry: 14
+B    0.0    0.0    1.166
+H    0.0    0.0    0.0
+O    0.0    0.0    2.366
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index 19d129d1..7ce1c0f6 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -169,6 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
 
     ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
     ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
+    ref_energy["vdz"]["HBO"] = Energy(None, -19.11982537379352)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #
@@ -378,7 +379,7 @@ class ValueTest(unittest.TestCase):
                                          pseudo=False))
 
     def test_hf(self):
-        self.assertTrue(run_hf(geo="SO2",
+        self.assertTrue(run_hf(geo="HBO",
                                basis="vdz",
                                mult=1,
                                pseudo=True))

From 7ccfee69a5f99e248f00851da3bbc6aba7036039 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 10:45:44 +0200
Subject: [PATCH 04/18] Add qp_utils.py for cache decorator, remove pickl in
 ei_handler.py

---
 scripts/ezfio_interface/ei_handler.py | 21 +++++----------------
 scripts/module/module_handler.py      | 19 +------------------
 scripts/qp_utils.py                   | 17 +++++++++++++++++
 src/Makefile                          |  4 ----
 4 files changed, 23 insertions(+), 38 deletions(-)
 create mode 100644 scripts/qp_utils.py

diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py
index 6d18d071..4667478f 100755
--- a/scripts/ezfio_interface/ei_handler.py
+++ b/scripts/ezfio_interface/ei_handler.py
@@ -62,6 +62,9 @@ import ConfigParser
 from collections import defaultdict
 from collections import namedtuple
 
+
+from qp_utils import cache
+
 Type = namedtuple('Type', 'fancy ocaml fortran')
 
 
@@ -78,6 +81,7 @@ def is_bool(str_):
         raise TypeError
 
 
+@cache
 def get_type_dict():
     """
     This function makes the correspondance between the type of value read in
@@ -89,17 +93,7 @@ def get_type_dict():
     # ~#~#~#~#~ #
     # P i c l e #
     # ~#~#~#~#~ #
-
-    import cPickle as pickle
-
-    from os import listdir
-
     qpackage_root = os.environ['QPACKAGE_ROOT']
-    fancy_type_pickle = qpackage_root + "/scripts/ezfio_interface/fancy_type.p"
-
-    if fancy_type_pickle in listdir(os.getcwd()):
-        fancy_type = pickle.load(open(fancy_type_pickle, "rb"))
-        return fancy_type
 
     # ~#~#~#~ #
     # I n i t #
@@ -148,9 +142,7 @@ def get_type_dict():
         b = r.find('let untouched = "')
         e = r.find(';;', b)
 
-        l_un = [
-            i for i in r[
-                b:e].splitlines() if i.strip().startswith("module")]
+        l_un = [i for i in r[b:e].splitlines() if i.strip().startswith("module")]
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # q p _ t y p e s _ g e n e r a t e #
@@ -174,9 +166,6 @@ def get_type_dict():
     # ~#~#~#~#~#~#~#~ #
     # F i n a l i z e #
     # ~#~#~#~#~#~#~#~ #
-
-    pickle.dump(dict(fancy_type), open(fancy_type_pickle, "wb"))
-
     return dict(fancy_type)
 
 
diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py
index 3bc1d0e4..a35dbe2c 100755
--- a/scripts/module/module_handler.py
+++ b/scripts/module/module_handler.py
@@ -24,24 +24,7 @@ from docopt import docopt
 import os
 import sys
 import os.path
-from functools import wraps
-
-
-def cache(func):
-    """
-    A decorator for lazy evaluation off true function
-    """
-    saved = {}
-
-    @wraps(func)
-    def newfunc(*args):
-        if args in saved:
-            return saved[args]
-
-        result = func(*args)
-        saved[args] = result
-        return result
-    return newfunc
+from qp_utils import cache
 
 
 @cache
diff --git a/scripts/qp_utils.py b/scripts/qp_utils.py
new file mode 100644
index 00000000..9be4ea2c
--- /dev/null
+++ b/scripts/qp_utils.py
@@ -0,0 +1,17 @@
+from functools import wraps
+
+def cache(func):
+    """
+    A decorator for lazy evaluation off true function
+    """
+    saved = {}
+
+    @wraps(func)
+    def newfunc(*args):
+        if args in saved:
+            return saved[args]
+
+        result = func(*args)
+        saved[args] = result
+        return result
+    return newfunc
diff --git a/src/Makefile b/src/Makefile
index 93b5234f..1bb5aa87 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -16,10 +16,6 @@ default: ezfio
 
 veryclean:
 	$(QPACKAGE_ROOT)/scripts/module/clean_modules.sh $(ALL_MODULES)
-	# Define the dict [type in EZFIO.cfg] = ocaml type , f90 type
-	# If you change the qptypes_generator.ml, you need to rm this
-	# For simplicity add this to the veryclean rule
-	rm -f $(QPACKAGE_ROOT)/scripts/ezfio_interface/fancy_type.p
 
 $(ALL_MODULES): ezfio
 	$(QPACKAGE_ROOT)/scripts/module/build_modules.sh $@ 

From 8dc7ea8481ba79094d2f073df9d0fd1f68f20d2e Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 11:31:05 +0200
Subject: [PATCH 05/18] Touch EZFIO.cfg and ezfio.config when make veryclean 
 for forcing the recompilation of the EZFIO

---
 scripts/module/clean_modules.sh | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/scripts/module/clean_modules.sh b/scripts/module/clean_modules.sh
index 608e0161..a11deb70 100755
--- a/scripts/module/clean_modules.sh
+++ b/scripts/module/clean_modules.sh
@@ -16,6 +16,8 @@ function do_clean()
     IRPF90_temp IRPF90_man Makefile.depend \
     $(module_handler.py print_genealogy) include \
     ezfio_interface.irp.f irpf90.make irpf90_entities tags $(ls_exe) *.mod
+
+  touch -c EZFIO.cfg *.ezfio_config
 }
 
 if [[ -z $1 ]]

From eeb60de46d861f0233177c9f906b28c127d3b378 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 11:56:37 +0200
Subject: [PATCH 06/18] solve overflow error in int.f90 for the pseudo

---
 src/Pseudo_integrals/int.f90 | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index be806b3f..46e204eb 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1879,13 +1879,18 @@ end
       RETURN
       END
 
-      double precision function coef_nk(n,k)
-      implicit none
-      integer n,k
-      double precision gam,dblefact,fact
-      gam=dblefact(2*(n+k)+1)
-      coef_nk=1.d0/(2.d0**k*fact(k)*gam)
-      end
+double precision function coef_nk(n,k)
+  implicit none
+  integer n,k
+  double precision gam,dblefact,fact
+
+  if(k.GE.80) then
+    coef_nk = 0.d0
+  else
+    gam=dblefact(2*(n+k)+1)
+    coef_nk=1.d0/(2.d0**k*fact(k)*gam)
+  endif
+end
 
 !!    Calculation of
 !!

From db8f905b98565a180d26b92244f8a7312ccdf827 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 13:39:11 +0200
Subject: [PATCH 07/18] Integration with the new code

---
 src/Pseudo_integrals/int.f90 | 165 +++++++++++++++++++----------------
 tests/unit_test/unit_test.py |   2 +-
 2 files changed, 90 insertions(+), 77 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 46e204eb..f7b35817 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -195,7 +195,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
 
 double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm
 double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big
-double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num_soph_p
+double precision :: areal,freal,breal,t1,t2,int_prod_bessel
 double precision :: arg
 
 integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot
@@ -276,7 +276,7 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
     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_num_soph_p(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
@@ -309,7 +309,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
    do lambdap=0,lmax+ntotB
         do k=1,kmax
           do l=0,lmax
-          array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
+          array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
           enddo
         enddo
    enddo
@@ -431,7 +431,7 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
     do k=1,kmax
       do l=0,lmax
 
-        array_R(ktot,k,l,0,lambdap)=  int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
+        array_R(ktot,k,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
        enddo
     enddo
   enddo
@@ -517,7 +517,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
     do k=1,kmax
       do l=0,lmax
 
-      array_R(ktot,k,l,lambda,0)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
+      array_R(ktot,k,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
       enddo
     enddo
   enddo
@@ -1898,14 +1898,32 @@ end
 !!
 !!    M_n(x) modified spherical bessel function 
 !!
-      double precision function int_prod_bessel(l,gam,n,m,a,b)
+      double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
       implicit none
       integer n,k,m,q,l,kcp
       double precision gam,dblefact,fact,pi,a,b
-      double precision int,intold,sum,coef_nk,crochet
+      double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
       logical done
 
+      u=(a+b)/(2.d0*dsqrt(gam))
+      freal=dexp(-arg)
+
+      if(a.eq.0.d0.and.b.eq.0.d0)then
+       if(n.ne.0.or.m.ne.0)then
+        int_prod_bessel=0.d0
+        return
+       endif
+       int_prod_bessel=crochet(l,gam)*freal
+       return
+      endif
+
+      if(u.gt.6.d0)then
+       int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+      return
+      endif
+
       if(a.ne.0.d0.and.b.ne.0.d0)then
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1925,17 +1943,12 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
-       if(kcp.gt.100) then
-          print*,"l",l
-          print*, "gam", gam
-          print*, "n", n
-          print*, "m", m
-          print*, "a", a
-          print*, "b", b
-          print*, "kcp", kcp
-          print*,'**WARNING** bad convergence in int_prod_bessel'
-        endif
+       int_prod_bessel=int*freal
+       if(kcp.gt.100)then
+        print*,'**WARNING** bad convergence in int_prod_bessel'
+!        else
+!        print*,'kcp=',kcp
+       endif
        return
       endif
 
@@ -1944,6 +1957,7 @@ end
         int_prod_bessel=0.d0
         return
        endif
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1959,7 +1973,7 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
+       int_prod_bessel=int*freal
        if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
        return
       endif
@@ -1969,6 +1983,7 @@ end
         int_prod_bessel=0.d0
         return
        endif
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1984,76 +1999,74 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
+       int_prod_bessel=int*freal
        if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
        return
-       endif
-
-      if(a.eq.0.d0.and.b.eq.0.d0)then
-       if(n.ne.0.or.m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-       int_prod_bessel=crochet(l,gam)
-       return
       endif
 
       stop 'pb in int_prod_bessel!!'
       end
 
-
-double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg)
+double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg)
       implicit none
-      integer :: n,m,l
-      double precision :: gam,a,b,arg,arg_new
-      double precision :: bessel_mod,factor
-      logical :: not_done
-      double precision :: bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal
-      integer :: iter, i, nI, n0
-      double precision :: eps
+      integer n,m,i,npts,l
+      double precision gam,a,b
+      double precision sum,x,bessel_mod,u,factor,arg
+      double precision xq(100),wq(100)
 
       u=(a+b)/(2.d0*dsqrt(gam))
-      arg_new=u**2-arg
-      freal=dexp(arg_new)
-      v=u/dsqrt(gam)
+      factor=dexp(u**2-arg)/dsqrt(gam)
 
-  bigA=v+dsqrt(-dlog(1.d-15)/gam)
-  n0=5
-  accu=0.d0
-  dx=bigA/(float(n0)-1.d0)
-  iter=0
-  do i=1,n0
-   x=(float(i)-1.d0)*dx
-   accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
-  enddo
+xq(1)= 5.38748089001123    
+xq(2)= 4.60368244955074   
+xq(3)= 3.94476404011563  
+xq(4)= 3.34785456738322 
+xq(5)= 2.78880605842813
+xq(6)= 2.25497400208928       
+xq(7)= 1.73853771211659      
+xq(8)= 1.23407621539532     
+xq(9)= 0.737473728545394   
+xq(10)= 0.245340708300901 
+xq(11)=-0.245340708300901
+xq(12)=-0.737473728545394    
+xq(13)=-1.23407621539532    
+xq(14)=-1.73853771211659   
+xq(15)=-2.25497400208928  
+xq(16)=-2.78880605842813 
+xq(17)=-3.34785456738322     
+xq(18)=-3.94476404011563    
+xq(19)=-4.60368244955074   
+xq(20)=-5.38748089001123  
+wq(1)=  2.229393645534151E-013
+wq(2)=  4.399340992273176E-010
+wq(3)=  1.086069370769280E-007
+wq(4)=  7.802556478532063E-006
+wq(5)=  2.283386360163528E-004
+wq(6)=  3.243773342237853E-003
+wq(7)=  2.481052088746362E-002
+wq(8)=  0.109017206020022
+wq(9)=  0.286675505362834
+wq(10)= 0.462243669600610
+wq(11)= 0.462243669600610
+wq(12)= 0.286675505362834
+wq(13)= 0.109017206020022
+wq(14)= 2.481052088746362E-002
+wq(15)= 3.243773342237853E-003
+wq(16)= 2.283386360163528E-004
+wq(17)= 7.802556478532063E-006
+wq(18)= 1.086069370769280E-007
+wq(19)= 4.399340992273176E-010
+wq(20)= 2.229393645534151E-013
 
-accu=accu*freal
-intold=accu*dx
+      npts=20
+!      call gauher(xq,wq,npts)
 
-eps=1.d-08
-nI=n0-1
-dx=dx/2.d0
-not_done=.true.
-
-do while(not_done)
- iter=iter+1
- accu=0.d0
- do i=1,nI
-  x=dx+(float(i)-1.d0)*2.d0*dx
-  accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
- enddo
- accu=accu*freal
- intnew=intold/2.d0+accu
- if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then
-  not_done=.false.
- else
-  intold2=intold
-  intold=intnew
-  dx=dx/2.d0
-  nI=2*nI
- endif
-enddo
-int_prod_bessel_num_soph_p=intold
+      sum=0.d0
+      do i=1,npts
+       x=(xq(i)+u)/dsqrt(gam)
+       sum=sum+wq(i)*x**l*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
+      enddo
+      int_prod_bessel_large=sum*factor
 end
 
 !!    Calculation of
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index 7ce1c0f6..b2881c92 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
 
     ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
     ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
-    ref_energy["vdz"]["HBO"] = Energy(None, -19.11982537379352)
+    ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #

From a8eaf1cfb43c62e4242288746fe19a35a7ff4bdc Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 14:22:38 +0200
Subject: [PATCH 08/18] Add indent

---
 src/Pseudo_integrals/int.f90 | 209 +++++++++++++++++++----------------
 tests/unit_test/unit_test.py |   2 +-
 2 files changed, 115 insertions(+), 96 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index f7b35817..4dffef09 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1898,114 +1898,133 @@ end
 !!
 !!    M_n(x) modified spherical bessel function 
 !!
-      double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
-      implicit none
-      integer n,k,m,q,l,kcp
-      double precision gam,dblefact,fact,pi,a,b
-      double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
-      logical done
 
-      u=(a+b)/(2.d0*dsqrt(gam))
-      freal=dexp(-arg)
+double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
-      if(a.eq.0.d0.and.b.eq.0.d0)then
-       if(n.ne.0.or.m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-       int_prod_bessel=crochet(l,gam)*freal
-       return
-      endif
+  implicit none
+  integer n,k,m,q,l,kcp
+  double precision gam,dblefact,fact,pi,a,b
+  double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
+  logical done
 
-      if(u.gt.6.d0)then
-       int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+  u=(a+b)/(2.d0*dsqrt(gam))
+  freal=dexp(-arg)
+
+  if(a.eq.0.d0.and.b.eq.0.d0)then
+    if(n.ne.0.or.m.ne.0)then
+      int_prod_bessel=0.d0
       return
+    endif
+    
+    int_prod_bessel=crochet(l,gam)*freal
+    return
+  endif
+
+  if(u.gt.6.d0)then
+   int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+    return
+  endif
+
+  if(a.ne.0.d0.and.b.ne.0.d0)then
+
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    do while (.not.done)  
+    
+      kcp=kcp+1
+      sum=0.d0
+     
+      do k=0,q
+        sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q)
+      enddo
+     
+     int=int+sum*crochet(2*q+n+m+l,gam)
+     
+     if(dabs(int-intold).lt.1d-15)then
+       done=.true.
+      else
+        
+        q=q+1
+        intold=int
       endif
+    enddo
+    
+    int_prod_bessel=int*freal
+    
+    if(kcp.gt.100)then
+      print*,'**WARNING** bad convergence in int_prod_bessel'
+    endif
+    
+    return
+  endif
 
-      if(a.ne.0.d0.and.b.ne.0.d0)then
+  if(a.eq.0.d0.and.b.ne.0.d0)then
+       
+    if(n.ne.0)then
+      int_prod_bessel=0.d0
+      return
+    endif
 
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)  
-        kcp=kcp+1
-        sum=0.d0
-        do k=0,q
-         sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q)
-        enddo
-        int=int+sum*crochet(2*q+n+m+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    do while (.not.done)
+        
+      kcp=kcp+1
+      int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
+      
+      if(dabs(int-intold).lt.1d-15)then
+        done=.true.
+      else
+        q=q+1
+        intold=int
+      endif
+    
+    enddo
+
+    int_prod_bessel=int*freal
+    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
+    return
+  endif
+
+  if(a.ne.0.d0.and.b.eq.0.d0)then
+    if(m.ne.0)then
+      int_prod_bessel=0.d0
+      return
+    endif
+
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    do while (.not.done)
+      kcp=kcp+1
+      int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
+      if(dabs(int-intold).lt.1d-15)then
          done=.true.
-        else
+      else
          q=q+1
          intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)then
-        print*,'**WARNING** bad convergence in int_prod_bessel'
-!        else
-!        print*,'kcp=',kcp
-       endif
-       return
       endif
+    enddo
+    
+    int_prod_bessel=int*freal
+    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
+    return
 
-      if(a.eq.0.d0.and.b.ne.0.d0)then
-       if(n.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
+  endif
 
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)
-        kcp=kcp+1
-        int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
-         done=.true.
-        else
-         q=q+1
-         intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
-       return
-      endif
-
-      if(a.ne.0.d0.and.b.eq.0.d0)then
-       if(m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)
-        kcp=kcp+1
-        int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
-         done=.true.
-        else
-         q=q+1
-         intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
-       return
-      endif
-
-      stop 'pb in int_prod_bessel!!'
-      end
+  stop 'pb in int_prod_bessel!!'
+end
 
 double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg)
       implicit none
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index b2881c92..23e831ec 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
 
     ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
     ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
-    ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160)
+    ref_energy["vdz"]["HBO"] = Energy(None, -19.11982540413317)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #

From 69c872333b2294560bfb39ada664ee93d7ca3dbf Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 15:19:17 +0200
Subject: [PATCH 09/18] >50% acceleration of int.f90

---
 src/Pseudo_integrals/int.f90 | 20 +++++++++++++++---
 tests/unit_test/unit_test.py | 40 ++++++++++++++++++------------------
 2 files changed, 37 insertions(+), 23 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 4dffef09..1fdb2c4d 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1905,6 +1905,10 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
   integer n,k,m,q,l,kcp
   double precision gam,dblefact,fact,pi,a,b
   double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
+  double precision dump
+
+  double precision, allocatable :: temp_array(:)
+
   logical done
 
   u=(a+b)/(2.d0*dsqrt(gam))
@@ -1933,13 +1937,20 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
     done=.false.
     kcp=0
 
+    allocate( temp_array(0:101))
+
     do while (.not.done)  
-    
+
       kcp=kcp+1
       sum=0.d0
-     
+
+      dump = a**n
+      do k=q,q+1
+         temp_array(k) = coef_nk(n,k)*dump*(a**(2*k))
+      enddo
+
       do k=0,q
-        sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q)
+        sum=sum+temp_array(k)*coef_nk(m,q-k)*b**(m-2*k+2*q)
       enddo
      
      int=int+sum*crochet(2*q+n+m+l,gam)
@@ -1951,8 +1962,11 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
         q=q+1
         intold=int
       endif
+
     enddo
     
+    deallocate(temp_array)
+
     int_prod_bessel=int*freal
     
     if(kcp.gt.100)then
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index 23e831ec..76dc6fd3 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -20,7 +20,7 @@ Energy = namedtuple('Energy', ['without_pseudo', 'with_pseudo'])
 # O p t #
 # ~#~#~ #
 
-precision = 1.e-7
+precision = 1.e-8
 
 # A test get a geo file and a basis file.
 # A global dict containt the result for this test
@@ -372,11 +372,11 @@ def check_convert(path_out):
 #
 class ValueTest(unittest.TestCase):
 
-    def test_hf_then_full_ci_10k_pt2_end(self):
-        self.assertTrue(hf_then_10k_test(geo="methane",
-                                         basis="sto-3g",
-                                         mult=1,
-                                         pseudo=False))
+#    def test_hf_then_full_ci_10k_pt2_end(self):
+#        self.assertTrue(hf_then_10k_test(geo="methane",
+#                                         basis="sto-3g",
+#                                         mult=1,
+#                                         pseudo=False))
 
     def test_hf(self):
         self.assertTrue(run_hf(geo="HBO",
@@ -385,20 +385,20 @@ class ValueTest(unittest.TestCase):
                                pseudo=True))
 
 
-class ConvertTest(unittest.TestCase):
-    def test_check_convert_hf_energy(self):
-        self.assertTrue(check_convert("HBO.out"))
-
-
-class InputTest(unittest.TestCase):
-
-    def test_check_disk_acess(self):
-        self.assertTrue(check_disk_acess(geo="methane",
-                                         basis="un-ccemd-ref"))
-
-    def test_check_mo_guess(self):
-        self.assertTrue(check_mo_guess(geo="methane",
-                                       basis="maug-cc-pVDZ"))
+#class ConvertTest(unittest.TestCase):
+#    def test_check_convert_hf_energy(self):
+#        self.assertTrue(check_convert("HBO.out"))
+#
+#
+#class InputTest(unittest.TestCase):
+#
+#    def test_check_disk_acess(self):
+#        self.assertTrue(check_disk_acess(geo="methane",
+#                                         basis="un-ccemd-ref"))
+#
+#    def test_check_mo_guess(self):
+#        self.assertTrue(check_mo_guess(geo="methane",
+#                                       basis="maug-cc-pVDZ"))
 
 if __name__ == '__main__':
     unittest.main()

From 9fb1161a557d9b83d8b04c52dd400b3a9cd722bc Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 15:19:17 +0200
Subject: [PATCH 10/18] >50% acceleration of int.f90

---
 src/Pseudo_integrals/int.f90 | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 4dffef09..1fdb2c4d 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1905,6 +1905,10 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
   integer n,k,m,q,l,kcp
   double precision gam,dblefact,fact,pi,a,b
   double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
+  double precision dump
+
+  double precision, allocatable :: temp_array(:)
+
   logical done
 
   u=(a+b)/(2.d0*dsqrt(gam))
@@ -1933,13 +1937,20 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
     done=.false.
     kcp=0
 
+    allocate( temp_array(0:101))
+
     do while (.not.done)  
-    
+
       kcp=kcp+1
       sum=0.d0
-     
+
+      dump = a**n
+      do k=q,q+1
+         temp_array(k) = coef_nk(n,k)*dump*(a**(2*k))
+      enddo
+
       do k=0,q
-        sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q)
+        sum=sum+temp_array(k)*coef_nk(m,q-k)*b**(m-2*k+2*q)
       enddo
      
      int=int+sum*crochet(2*q+n+m+l,gam)
@@ -1951,8 +1962,11 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
         q=q+1
         intold=int
       endif
+
     enddo
     
+    deallocate(temp_array)
+
     int_prod_bessel=int*freal
     
     if(kcp.gt.100)then

From c1c684b0e731cb5a8445735c70bea8159ea39947 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 14:22:38 +0200
Subject: [PATCH 11/18] >50% acceleration of int.f90

---
 src/Pseudo_integrals/int.f90 | 223 ++++++++++++++++++++---------------
 tests/unit_test/unit_test.py |   2 +-
 2 files changed, 129 insertions(+), 96 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index f7b35817..1fdb2c4d 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1898,114 +1898,147 @@ end
 !!
 !!    M_n(x) modified spherical bessel function 
 !!
-      double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
-      implicit none
-      integer n,k,m,q,l,kcp
-      double precision gam,dblefact,fact,pi,a,b
-      double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
-      logical done
 
-      u=(a+b)/(2.d0*dsqrt(gam))
-      freal=dexp(-arg)
+double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
-      if(a.eq.0.d0.and.b.eq.0.d0)then
-       if(n.ne.0.or.m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-       int_prod_bessel=crochet(l,gam)*freal
-       return
-      endif
+  implicit none
+  integer n,k,m,q,l,kcp
+  double precision gam,dblefact,fact,pi,a,b
+  double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
+  double precision dump
 
-      if(u.gt.6.d0)then
-       int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+  double precision, allocatable :: temp_array(:)
+
+  logical done
+
+  u=(a+b)/(2.d0*dsqrt(gam))
+  freal=dexp(-arg)
+
+  if(a.eq.0.d0.and.b.eq.0.d0)then
+    if(n.ne.0.or.m.ne.0)then
+      int_prod_bessel=0.d0
       return
+    endif
+    
+    int_prod_bessel=crochet(l,gam)*freal
+    return
+  endif
+
+  if(u.gt.6.d0)then
+   int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+    return
+  endif
+
+  if(a.ne.0.d0.and.b.ne.0.d0)then
+
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    allocate( temp_array(0:101))
+
+    do while (.not.done)  
+
+      kcp=kcp+1
+      sum=0.d0
+
+      dump = a**n
+      do k=q,q+1
+         temp_array(k) = coef_nk(n,k)*dump*(a**(2*k))
+      enddo
+
+      do k=0,q
+        sum=sum+temp_array(k)*coef_nk(m,q-k)*b**(m-2*k+2*q)
+      enddo
+     
+     int=int+sum*crochet(2*q+n+m+l,gam)
+     
+     if(dabs(int-intold).lt.1d-15)then
+       done=.true.
+      else
+        
+        q=q+1
+        intold=int
       endif
 
-      if(a.ne.0.d0.and.b.ne.0.d0)then
+    enddo
+    
+    deallocate(temp_array)
 
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)  
-        kcp=kcp+1
-        sum=0.d0
-        do k=0,q
-         sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q)
-        enddo
-        int=int+sum*crochet(2*q+n+m+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
+    int_prod_bessel=int*freal
+    
+    if(kcp.gt.100)then
+      print*,'**WARNING** bad convergence in int_prod_bessel'
+    endif
+    
+    return
+  endif
+
+  if(a.eq.0.d0.and.b.ne.0.d0)then
+       
+    if(n.ne.0)then
+      int_prod_bessel=0.d0
+      return
+    endif
+
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    do while (.not.done)
+        
+      kcp=kcp+1
+      int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
+      
+      if(dabs(int-intold).lt.1d-15)then
+        done=.true.
+      else
+        q=q+1
+        intold=int
+      endif
+    
+    enddo
+
+    int_prod_bessel=int*freal
+    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
+    return
+  endif
+
+  if(a.ne.0.d0.and.b.eq.0.d0)then
+    if(m.ne.0)then
+      int_prod_bessel=0.d0
+      return
+    endif
+
+    q=0
+    intold=-1.d0
+    int=0.d0
+    done=.false.
+    kcp=0
+
+    do while (.not.done)
+      kcp=kcp+1
+      int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
+      if(dabs(int-intold).lt.1d-15)then
          done=.true.
-        else
+      else
          q=q+1
          intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)then
-        print*,'**WARNING** bad convergence in int_prod_bessel'
-!        else
-!        print*,'kcp=',kcp
-       endif
-       return
       endif
+    enddo
+    
+    int_prod_bessel=int*freal
+    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
+    return
 
-      if(a.eq.0.d0.and.b.ne.0.d0)then
-       if(n.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
+  endif
 
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)
-        kcp=kcp+1
-        int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
-         done=.true.
-        else
-         q=q+1
-         intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
-       return
-      endif
-
-      if(a.ne.0.d0.and.b.eq.0.d0)then
-       if(m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-
-       q=0
-       intold=-1.d0
-       int=0.d0
-       done=.false.
-       kcp=0
-       do while (.not.done)
-        kcp=kcp+1
-        int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
-        if(dabs(int-intold).lt.1d-15)then
-         done=.true.
-        else
-         q=q+1
-         intold=int
-        endif
-       enddo
-       int_prod_bessel=int*freal
-       if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
-       return
-      endif
-
-      stop 'pb in int_prod_bessel!!'
-      end
+  stop 'pb in int_prod_bessel!!'
+end
 
 double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg)
       implicit none
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index b2881c92..23e831ec 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
 
     ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
     ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
-    ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160)
+    ref_energy["vdz"]["HBO"] = Energy(None, -19.11982540413317)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #

From 0f26522eedf16eb414bbaf5d73ca7474a3cce14b Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 15:41:10 +0200
Subject: [PATCH 12/18] Restore unit test

---
 tests/unit_test/unit_test.py | 40 ++++++++++++++++++------------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index 76dc6fd3..23e831ec 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -20,7 +20,7 @@ Energy = namedtuple('Energy', ['without_pseudo', 'with_pseudo'])
 # O p t #
 # ~#~#~ #
 
-precision = 1.e-8
+precision = 1.e-7
 
 # A test get a geo file and a basis file.
 # A global dict containt the result for this test
@@ -372,11 +372,11 @@ def check_convert(path_out):
 #
 class ValueTest(unittest.TestCase):
 
-#    def test_hf_then_full_ci_10k_pt2_end(self):
-#        self.assertTrue(hf_then_10k_test(geo="methane",
-#                                         basis="sto-3g",
-#                                         mult=1,
-#                                         pseudo=False))
+    def test_hf_then_full_ci_10k_pt2_end(self):
+        self.assertTrue(hf_then_10k_test(geo="methane",
+                                         basis="sto-3g",
+                                         mult=1,
+                                         pseudo=False))
 
     def test_hf(self):
         self.assertTrue(run_hf(geo="HBO",
@@ -385,20 +385,20 @@ class ValueTest(unittest.TestCase):
                                pseudo=True))
 
 
-#class ConvertTest(unittest.TestCase):
-#    def test_check_convert_hf_energy(self):
-#        self.assertTrue(check_convert("HBO.out"))
-#
-#
-#class InputTest(unittest.TestCase):
-#
-#    def test_check_disk_acess(self):
-#        self.assertTrue(check_disk_acess(geo="methane",
-#                                         basis="un-ccemd-ref"))
-#
-#    def test_check_mo_guess(self):
-#        self.assertTrue(check_mo_guess(geo="methane",
-#                                       basis="maug-cc-pVDZ"))
+class ConvertTest(unittest.TestCase):
+    def test_check_convert_hf_energy(self):
+        self.assertTrue(check_convert("HBO.out"))
+
+
+class InputTest(unittest.TestCase):
+
+    def test_check_disk_acess(self):
+        self.assertTrue(check_disk_acess(geo="methane",
+                                         basis="un-ccemd-ref"))
+
+    def test_check_mo_guess(self):
+        self.assertTrue(check_mo_guess(geo="methane",
+                                       basis="maug-cc-pVDZ"))
 
 if __name__ == '__main__':
     unittest.main()

From 92258bb7a3226d5cd37fefc763ab673a5779db6e Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 16:00:35 +0200
Subject: [PATCH 13/18] Remove test_int useless function

---
 src/Pseudo_integrals/int.f90 | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 1fdb2c4d..f634c510 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1835,24 +1835,6 @@ double precision function binom_gen(alpha,n)
  enddo
 end
 
-double precision function test_int(g_a,g_b,g_c,ac,bc)
-implicit none
-double precision factor,g_a,g_b,g_c,ac,bc,x,dx,sum,alpha,beta,pi
-integer i,npts
-pi=dacos(-1.d0)
-factor=0.5d0*pi/(g_a*g_b*ac*bc*dsqrt(g_a+g_b+g_c))*dexp(-g_a*ac**2-g_b*bc**2)
-npts=2000
-dx=20.d0/npts
-sum=0.d0
-alpha=(2.d0*g_a*ac+2.d0*g_b*bc)/dsqrt(g_c+g_a+g_b)
-beta=(2.d0*g_b*bc-2.d0*g_b*bc)/dsqrt(g_c+g_a+g_b)
-do i=1,npts
- x=(i-1)*dx+0.5d0*dx
- sum=sum+dx*dexp(-x**2)*(dcosh(alpha*x)-dcosh(beta*x))
-enddo
-test_int=factor*sum
-end
-
 recursive function fact1(n,a) result(x)
 implicit none
 integer n

From 17155b5ec19100ccdfa1aa190eefa55706abea97 Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 19:13:32 +0200
Subject: [PATCH 14/18] Optimise a litle the pseudo(10-20%)

---
 src/Pseudo_integrals/int.f90 | 124 ++++++++++++++++++++---------------
 src/Utils/util.irp.f         |  41 +++++++++++-
 2 files changed, 110 insertions(+), 55 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index f634c510..49677e02 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -800,13 +800,13 @@ end
 double precision function bigA(i,j,k)
 implicit none
 integer i,j,k
-double precision fourpi,dblefact
+double precision fourpi,dble_fact
 fourpi=4.d0*dacos(-1.d0)
 bigA=0.d0
 if(mod(i,2).eq.1)return
 if(mod(j,2).eq.1)return
 if(mod(k,2).eq.1)return
-bigA=fourpi*dblefact(i-1)*dblefact(j-1)*dblefact(k-1)/dblefact(i+j+k+1)
+bigA=fourpi*dble_fact(i-1)*dble_fact(j-1)*dble_fact(k-1)/dble_fact(i+j+k+1)
 end
 !!
 !! I_{lambda,mu,l,m}^{k1,k2,k3} = /int dOmega  Y_{lambda mu} xchap^k1 ychap^k2 zchap^k3  Y_{lm}
@@ -996,10 +996,10 @@ end
 double precision function crochet(n,g)
 implicit none
 integer n
-double precision g,pi,dblefact,expo
+double precision g,pi,dble_fact,expo
 pi=dacos(-1.d0)
 expo=0.5d0*dfloat(n+1)
-crochet=dblefact(n-1)/(2.d0*g)**expo
+crochet=dble_fact(n-1)/(2.d0*g)**expo
 if(mod(n,2).eq.0)crochet=crochet*dsqrt(pi/2.d0)
 end
 
@@ -1431,30 +1431,6 @@ end
        stop
       end
 
-double precision function dblefact(n)
-implicit none
-integer :: n,k
-double precision prod
-dblefact=1.d0
-
-if(n.lt.0)return
-if(mod(n,2).eq.1)then
-  prod=1.d0
-  do k=1,n,2
-   prod=prod*dfloat(k)
-  enddo
-  dblefact=prod
-  return
- endif
- if(mod(n,2).eq.0)then
-  prod=1.d0
-  do k=2,n,2
-   prod=prod*dfloat(k)
-  enddo
-  dblefact=prod
-  return
- endif
-end
 !!
 !! R_{lambda,lamba',N}= exp(-ga_a AC**2 -g_b BC**2) /int_{0}{+infty} r**(2+n) exp(-(g_a+g_b+g_k)r**2)
 !!                       * M_{lambda}( 2g_a ac r) M_{lambda'}(2g_b bc r)
@@ -1510,10 +1486,10 @@ end
        double precision function bessel_mod_exp(n,x)
        implicit none
        integer n,k
-       double precision x,coef,accu,fact,dblefact
+       double precision x,coef,accu,fact,dble_fact
        accu=0.d0
        do k=0,10
-        coef=1.d0/fact(k)/dblefact(2*(n+k)+1)
+        coef=1.d0/fact(k)/dble_fact(2*(n+k)+1)
         accu=accu+(x**2/2.d0)**k*coef
        enddo
        bessel_mod_exp=x**n*accu
@@ -1835,22 +1811,6 @@ double precision function binom_gen(alpha,n)
  enddo
 end
 
-recursive function fact1(n,a) result(x)
-implicit none
-integer n
-double precision a,x,erf
-if(n.eq.0)then
-x=dsqrt(dacos(-1.d0))/2.d0*erf(a)
-return
-endif
-if(n.eq.1)then
-x=1.d0-dexp(-a**2)
-return
-endif
-if(mod(n,2).eq.0)x=0.5d0*dfloat(n-1)*fact1(n-2,a)+a**n*dexp(-a**2)
-if(mod(n,2).eq.1)x=0.5d0*dfloat(n-1)*fact1(n-2,a)+0.5d0*a**(n-1)*dexp(-a**2)
-end
-
       double precision FUNCTION ERF(X)
       implicit double precision(a-h,o-z)
       IF(X.LT.0.d0)THEN
@@ -1863,15 +1823,71 @@ end
 
 double precision function coef_nk(n,k)
   implicit none
-  integer n,k
-  double precision gam,dblefact,fact
+  integer n,k, ISHFT
+
+  double precision gam,dble_fact,fact
+
+  double precision, save :: store_result(0:2,0:10)
+  integer, save :: ifirst
+
+  if (ifirst == 0) then
+    ifirst =1
 
-  if(k.GE.80) then
-    coef_nk = 0.d0
-  else
-    gam=dblefact(2*(n+k)+1)
-    coef_nk=1.d0/(2.d0**k*fact(k)*gam)
   endif
+
+  store_result(0, 0) =  1.00000000000000d0
+  store_result(0, 1) =  0.166666666666667d0
+  store_result(0, 2) =  8.333333333333333d-3
+  store_result(0, 3) =  1.984126984126984d-4
+  store_result(0, 4) =  2.755731922398589d-6
+  store_result(0, 5) =  2.505210838544172d-8
+  store_result(0, 6) =  1.605904383682161d-10
+  store_result(0, 7) =  7.647163731819816d-13
+  store_result(0, 8) =  2.811457254345521d-15
+  store_result(0, 9) =  8.220635246624329d-018
+  store_result(0,10) =  1.957294106339126d-020
+
+  store_result(1, 0) =  0.333333333333333d0     
+  store_result(1, 1) =  3.333333333333333d-002   
+  store_result(1, 2) =  1.190476190476191d-003
+  store_result(1, 3) =  2.204585537918871d-005   
+  store_result(1, 4) =  2.505210838544172d-007   
+  store_result(1, 5) =  1.927085260418594d-009   
+  store_result(1, 6) =  1.070602922454774d-011   
+  store_result(1, 7) =  4.498331606952833d-014   
+  store_result(1, 8) =  1.479714344392379d-016
+  store_result(1, 9) =  3.914588212678252d-019
+  store_result(1,10) =  8.509974375387505d-022
+
+  store_result(2, 0) =  6.666666666666667d-002
+  store_result(2, 1) =  4.761904761904762d-003
+  store_result(2, 2) =  1.322751322751323d-004
+  store_result(2, 3) =  2.004168670835337d-006
+  store_result(2, 4) =  1.927085260418594d-008
+  store_result(2, 5) =  1.284723506945729d-010
+  store_result(2, 6) =  6.297664249733966d-013
+  store_result(2, 7) =  2.367542951027807d-015
+  store_result(2, 8) =  7.046258782820854d-018
+  store_result(2, 9) =  1.701994875077501d-020
+  store_result(2,10) =  3.403989750155002d-023
+
+  if (n.LE.2) then
+    if (k.LE.10) then
+      coef_nk = store_result(n,k)
+      return 
+    endif
+  endif
+
+  if ((n+k).GE.30) then
+    coef_nk = 0.d0
+    return
+  endif
+
+  gam=dble_fact(2*(n+k)+1)
+  coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
+  
+  return
+
 end
 
 !!    Calculation of
@@ -1885,7 +1901,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
   implicit none
   integer n,k,m,q,l,kcp
-  double precision gam,dblefact,fact,pi,a,b
+  double precision gam,dble_fact,fact,pi,a,b
   double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
   double precision dump
 
diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f
index 1186d789..7a78ad59 100644
--- a/src/Utils/util.irp.f
+++ b/src/Utils/util.irp.f
@@ -134,7 +134,46 @@ BEGIN_PROVIDER [ double precision, fact_inv, (128) ]
   enddo
 END_PROVIDER
 
-double precision function dble_fact(n) result(fact2)
+
+double precision function dble_fact(n)
+  implicit none
+  integer :: n
+  double precision :: dble_fact_peer, dble_fact_odd
+
+  dble_fact = 1.d0
+
+  if(n.lt.0) return
+
+  if(iand(n,1).eq.0)then
+    dble_fact = dble_fact_peer(n)
+  else
+    dble_fact= dble_fact_odd(n)
+  endif
+
+end function
+
+double precision function dble_fact_peer(n) result(fact2)
+  implicit none
+  BEGIN_DOC
+  ! n!!
+  END_DOC
+  integer                        :: n,k
+  double precision, save         :: memo(1:100)
+  integer, save                  :: memomax = 2
+  double precision               :: prod
+
+  ASSERT (iand(n,1) /= 1)
+
+  prod=1.d0
+  do k=2,n,2
+   prod=prod*dfloat(k)
+  enddo
+  fact2=prod
+  return
+
+end function
+
+double precision function dble_fact_odd(n) result(fact2)
   implicit none
   BEGIN_DOC
   ! n!!

From bd95f3f86ea5a0c26d7bdd559d901e7478f7b2bf Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 19:58:33 +0200
Subject: [PATCH 15/18] Gain a Major 50% in int integral

---
 src/Pseudo_integrals/int.f90 | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 49677e02..6f0d4988 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1905,7 +1905,8 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
   double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
   double precision dump
 
-  double precision, allocatable :: temp_array(:)
+  double precision, allocatable :: temp_array_a(:)
+  double precision, allocatable :: temp_array_b(:)
 
   logical done
 
@@ -1935,22 +1936,21 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
     done=.false.
     kcp=0
 
-    allocate( temp_array(0:101))
+    allocate( temp_array_a(0:100))
+    allocate( temp_array_b(0+m:200+m))
 
     do while (.not.done)  
 
       kcp=kcp+1
       sum=0.d0
 
-      dump = a**n
-      do k=q,q+1
-         temp_array(k) = coef_nk(n,k)*dump*(a**(2*k))
-      enddo
+      temp_array_a(q) = coef_nk(n,q)*(a**(n+2*q))
+      temp_array_b(m+2*q) = coef_nk(m,q)*b**(m+2*q)
 
       do k=0,q
-        sum=sum+temp_array(k)*coef_nk(m,q-k)*b**(m-2*k+2*q)
+        sum=sum+temp_array_a(k)*temp_array_b(m-2*k+2*q)
       enddo
-     
+
      int=int+sum*crochet(2*q+n+m+l,gam)
      
      if(dabs(int-intold).lt.1d-15)then
@@ -1963,14 +1963,15 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
     enddo
     
-    deallocate(temp_array)
+    deallocate(temp_array_a)
+    deallocate(temp_array_b)
 
     int_prod_bessel=int*freal
     
     if(kcp.gt.100)then
       print*,'**WARNING** bad convergence in int_prod_bessel'
     endif
-    
+
     return
   endif
 

From 9452013782815966c6a6f3c0504e878ea1ea37eb Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Tue, 5 May 2015 20:58:41 +0200
Subject: [PATCH 16/18] More clarity

---
 src/Pseudo_integrals/int.f90 | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 6f0d4988..bade7eb2 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1828,12 +1828,6 @@ double precision function coef_nk(n,k)
   double precision gam,dble_fact,fact
 
   double precision, save :: store_result(0:2,0:10)
-  integer, save :: ifirst
-
-  if (ifirst == 0) then
-    ifirst =1
-
-  endif
 
   store_result(0, 0) =  1.00000000000000d0
   store_result(0, 1) =  0.166666666666667d0
@@ -1878,7 +1872,7 @@ double precision function coef_nk(n,k)
     endif
   endif
 
-  if ((n+k).GE.30) then
+  if ((n+k).GE.60) then
     coef_nk = 0.d0
     return
   endif
@@ -1945,10 +1939,10 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
       sum=0.d0
 
       temp_array_a(q) = coef_nk(n,q)*(a**(n+2*q))
-      temp_array_b(m+2*q) = coef_nk(m,q)*b**(m+2*q)
+      temp_array_b(q) = coef_nk(m,q)*b**(m+2*q)
 
       do k=0,q
-        sum=sum+temp_array_a(k)*temp_array_b(m-2*k+2*q)
+        sum=sum+temp_array_a(k)*temp_array_b(q-2*k)
       enddo
 
      int=int+sum*crochet(2*q+n+m+l,gam)

From 2ec05fd08dafee17af4bf701c93cbd0b76019c1d Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Thu, 7 May 2015 11:15:14 +0200
Subject: [PATCH 17/18] Remove overflow and 400% acceleration

---
 src/Pseudo_integrals/int.f90 | 233 +++++++++++++----------------------
 1 file changed, 88 insertions(+), 145 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index bade7eb2..7db0dcd8 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -430,7 +430,6 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
   do lambdap=0,lmax+ntotB
     do k=1,kmax
       do l=0,lmax
-
         array_R(ktot,k,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
        enddo
     enddo
@@ -1827,58 +1826,11 @@ double precision function coef_nk(n,k)
 
   double precision gam,dble_fact,fact
 
-  double precision, save :: store_result(0:2,0:10)
-
-  store_result(0, 0) =  1.00000000000000d0
-  store_result(0, 1) =  0.166666666666667d0
-  store_result(0, 2) =  8.333333333333333d-3
-  store_result(0, 3) =  1.984126984126984d-4
-  store_result(0, 4) =  2.755731922398589d-6
-  store_result(0, 5) =  2.505210838544172d-8
-  store_result(0, 6) =  1.605904383682161d-10
-  store_result(0, 7) =  7.647163731819816d-13
-  store_result(0, 8) =  2.811457254345521d-15
-  store_result(0, 9) =  8.220635246624329d-018
-  store_result(0,10) =  1.957294106339126d-020
-
-  store_result(1, 0) =  0.333333333333333d0     
-  store_result(1, 1) =  3.333333333333333d-002   
-  store_result(1, 2) =  1.190476190476191d-003
-  store_result(1, 3) =  2.204585537918871d-005   
-  store_result(1, 4) =  2.505210838544172d-007   
-  store_result(1, 5) =  1.927085260418594d-009   
-  store_result(1, 6) =  1.070602922454774d-011   
-  store_result(1, 7) =  4.498331606952833d-014   
-  store_result(1, 8) =  1.479714344392379d-016
-  store_result(1, 9) =  3.914588212678252d-019
-  store_result(1,10) =  8.509974375387505d-022
-
-  store_result(2, 0) =  6.666666666666667d-002
-  store_result(2, 1) =  4.761904761904762d-003
-  store_result(2, 2) =  1.322751322751323d-004
-  store_result(2, 3) =  2.004168670835337d-006
-  store_result(2, 4) =  1.927085260418594d-008
-  store_result(2, 5) =  1.284723506945729d-010
-  store_result(2, 6) =  6.297664249733966d-013
-  store_result(2, 7) =  2.367542951027807d-015
-  store_result(2, 8) =  7.046258782820854d-018
-  store_result(2, 9) =  1.701994875077501d-020
-  store_result(2,10) =  3.403989750155002d-023
-
-  if (n.LE.2) then
-    if (k.LE.10) then
-      coef_nk = store_result(n,k)
-      return 
-    endif
-  endif
-
-  if ((n+k).GE.60) then
-    coef_nk = 0.d0
-    return
-  endif
-
   gam=dble_fact(2*(n+k)+1)
-  coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
+
+!  coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
+  
+  coef_nk=1.d0/(2.d0**k*fact(k)*gam)
   
   return
 
@@ -1897,10 +1849,11 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
   integer n,k,m,q,l,kcp
   double precision gam,dble_fact,fact,pi,a,b
   double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
-  double precision dump
 
-  double precision, allocatable :: temp_array_a(:)
-  double precision, allocatable :: temp_array_b(:)
+  integer :: n_1,m_1,nlm
+  double precision :: term_A, term_B, term_rap,  expo
+  double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square
+  double precision :: int_prod_bessel_loc
 
   logical done
 
@@ -1928,41 +1881,58 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
     intold=-1.d0
     int=0.d0
     done=.false.
-    kcp=0
 
-    allocate( temp_array_a(0:100))
-    allocate( temp_array_b(0+m:200+m))
+    n_1 = 2*(n)+1
+    m_1 = 2*m+1
+    nlm = n+m+l
+    pi=dacos(-1.d0)
+    a_over_b_square = (a/b)**2
 
+    ! Calcul first term of the sequence
+
+    term_a =dble_fact(nlm-1) / (dble_fact(n_1)*dble_fact(m_1))
+    expo=0.5d0*dfloat(nlm+1)
+    term_rap = term_a / (2.d0*gam)**expo
+
+    s_0_0=term_rap*a**(n)*b**(m)
+    if(mod(nlm,2).eq.0)s_0_0=s_0_0*dsqrt(pi/2.d0)
+
+    ! Initialise the first recurence terme for the q loop
+    s_q_0 = s_0_0
+
+    ! Loop over q for the convergence of the sequence
     do while (.not.done)  
 
-      kcp=kcp+1
-      sum=0.d0
-
-      temp_array_a(q) = coef_nk(n,q)*(a**(n+2*q))
-      temp_array_b(q) = coef_nk(m,q)*b**(m+2*q)
+      ! Init
+      sum=0
+      s_q_k=s_q_0
 
+      ! Iteration of k
       do k=0,q
-        sum=sum+temp_array_a(k)*temp_array_b(q-2*k)
+        sum=sum+s_q_k
+        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)/dble(2*(k+n)+3) ) * ( dble(q-k)/dble(k+1)) * s_q_k
       enddo
+    
+      int=int+sum
 
-     int=int+sum*crochet(2*q+n+m+l,gam)
-     
-     if(dabs(int-intold).lt.1d-15)then
-       done=.true.
+      if(dabs(int-intold).lt.1d-15)then
+        done=.true.
       else
+
+        !Compute the s_q+1_0
+        s_q_0=s_q_0*(2.d0*q+nlm+1)*b**2/((2.d0*(m+q)+3)*4.d0*(q+1)*gam)
         
+        if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi/2.d0)
+        ! Increment q
         q=q+1
         intold=int
       endif
 
     enddo
-    
-    deallocate(temp_array_a)
-    deallocate(temp_array_b)
 
     int_prod_bessel=int*freal
-    
-    if(kcp.gt.100)then
+
+    if(q.gt.100)then
       print*,'**WARNING** bad convergence in int_prod_bessel'
     endif
 
@@ -1971,61 +1941,15 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
   if(a.eq.0.d0.and.b.ne.0.d0)then
        
-    if(n.ne.0)then
-      int_prod_bessel=0.d0
-      return
-    endif
-
-    q=0
-    intold=-1.d0
-    int=0.d0
-    done=.false.
-    kcp=0
-
-    do while (.not.done)
-        
-      kcp=kcp+1
-      int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
-      
-      if(dabs(int-intold).lt.1d-15)then
-        done=.true.
-      else
-        q=q+1
-        intold=int
-      endif
-    
-    enddo
-
+    int = int_prod_bessel_loc(l,gam,m,b)
     int_prod_bessel=int*freal
-    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
     return
   endif
 
   if(a.ne.0.d0.and.b.eq.0.d0)then
-    if(m.ne.0)then
-      int_prod_bessel=0.d0
-      return
-    endif
 
-    q=0
-    intold=-1.d0
-    int=0.d0
-    done=.false.
-    kcp=0
-
-    do while (.not.done)
-      kcp=kcp+1
-      int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
-      if(dabs(int-intold).lt.1d-15)then
-         done=.true.
-      else
-         q=q+1
-         intold=int
-      endif
-    enddo
-    
+    int = int_prod_bessel_loc(l,gam,n,a)
     int_prod_bessel=int*freal
-    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
     return
 
   endif
@@ -2101,30 +2025,49 @@ end
 !!
 !!    M_n(x) modified spherical bessel function 
 !!
-      double precision function int_prod_bessel_loc(l,gam,n,a)
-      implicit none
-      integer n,k,l,kcp
-      double precision gam,a
-      double precision int,intold,coef_nk,crochet
-      logical done
-      k=0
-      intold=-1.d0
-      int=0.d0
-      done=.false.
-      kcp=0
-      do while (.not.done)
-       kcp=kcp+1
-       int=int+coef_nk(n,k)*a**(n+2*k)*crochet(2*k+n+l,gam)
-       if(dabs(int-intold).lt.1d-15)then
-        done=.true.
-       else
-        k=k+1
-        intold=int
-       endif
-      enddo
-      int_prod_bessel_loc=int
-      if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
-      end
+double precision function int_prod_bessel_loc(l,gam,n,a)
+  implicit none
+  integer n,k,l,kcp
+  double precision gam,a
+  double precision int,intold,coef_nk,crochet,dble_fact, fact, pi, expo
+  double precision :: f_0, f_k
+  logical done
+  
+  pi=dacos(-1.d0)
+  intold=-1.d0
+  done=.false.
+  int=0
+  
+  ! Int f_0
+  coef_nk=1.d0/dble_fact(2*n+1)
+  expo=0.5d0*dfloat(n+l+1)
+  crochet=dble_fact(n+l-1)/(2.d0*gam)**expo
+  if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(pi/2.d0)
+
+  f_0 = coef_nk*a**n*crochet
+
+  k=0
+
+  f_k=f_0
+  do while (.not.done)
+
+    int=int+f_k
+
+    f_k = f_k*(a**2*(2*(k+1)+n+l-1)) / (2*(k+1)*(2*(n+k+1)+1)*2*gam)
+
+    if(dabs(int-intold).lt.1d-15)then
+      done=.true.
+    else
+      k=k+1
+      intold=int
+    endif
+
+  enddo
+  
+  int_prod_bessel_loc=int
+  if(k.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
+
+end
 
       double precision function int_prod_bessel_num(l,gam,n,m,a,b)
       implicit none

From 8147c04a21f21e61c4f0e8a48230712e04b22c8b Mon Sep 17 00:00:00 2001
From: Thomas Applencourt <t.applencourt@gmail.com>
Date: Thu, 7 May 2015 15:41:59 +0200
Subject: [PATCH 18/18] Remove warning print for convergance

---
 src/Pseudo_integrals/int.f90 | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 7db0dcd8..3c247918 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -1932,10 +1932,6 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
 
     int_prod_bessel=int*freal
 
-    if(q.gt.100)then
-      print*,'**WARNING** bad convergence in int_prod_bessel'
-    endif
-
     return
   endif
 
@@ -2065,8 +2061,6 @@ double precision function int_prod_bessel_loc(l,gam,n,a)
   enddo
   
   int_prod_bessel_loc=int
-  if(k.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
-
 end
 
       double precision function int_prod_bessel_num(l,gam,n,m,a,b)