diff --git a/.travis.yml b/.travis.yml
index 8d3938a3..a4e27b3f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -12,7 +12,7 @@ before_script:
 script: 
   - ./setup_environment.sh --robot 
   - source ./quantum_package.rc
-  - ./scripts/compilation/create_ninja_build.py --production ./config/gfortran_example.cfg
-  - ./ninja/ninja
+  - qp_create_ninja.py --production ./config/gfortran_example.cfg
+  - ninja
   - cd ocaml ; make ; cd -
   - cd testing_no_regression ; ./unit_test.py
diff --git a/scripts/compilation/create_ninja_build.py b/scripts/compilation/qp_create_ninja.py
similarity index 94%
rename from scripts/compilation/create_ninja_build.py
rename to scripts/compilation/qp_create_ninja.py
index 47de8120..450193b4 100755
--- a/scripts/compilation/create_ninja_build.py
+++ b/scripts/compilation/qp_create_ninja.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 """
-Usage: create_ninja_build.py (--development | --production) CONFIG_FILE
+Usage: qp_create_ninja.py (--development | --production) CONFIG_FILE
 """
 
 import os
@@ -12,7 +12,7 @@ from collections import namedtuple
 from collections import defaultdict
 
 try:
-    from module_handler import get_dict_module_boss, get_dict_genealogy_desc
+    from module_handler import ModuleHandler
     from read_compilation_cfg import get_compilation_option
     from docopt import docopt
 except ImportError:
@@ -307,26 +307,21 @@ def get_l_file_for_module(path_module_abs):
     l_template = []
 
     for f in os.listdir(path_module_abs):
-        if f.endswith(".template.f"):
+        if f.lower().endswith(tuple([".template.f", ".include.f"])):
             l_template.append(join(path_module_abs, f))
         elif f.endswith(".irp.f"):
             l_irp.append(join(path_module_abs, f))
-        elif any(f.endswith("{0}".format(i)) for i in [".f",
-                                                       ".f90",
-                                                       ".c",
-                                                       ".cpp",
-                                                       ".cxx"]):
+        elif f.lower().endswith(tuple([".f", ".f90", ".c", ".cpp", ".cxx"])):
             l_src.append(join(path_module_abs, f))
             obj = '{0}.o'.format(os.path.splitext(f)[0])
             l_obj.append(join(path_module_abs, obj))
         elif f == "EZFIO.cfg":
             l_irp.append(join(path_module_abs, "ezfio_interface.irp.f"))
 
-    d = dict()
-    d["l_irp"] = l_irp
-    d["l_src"] = l_src
-    d["l_obj"] = l_obj
-    d["l_template"] = l_template
+    d = {"l_irp": l_irp,
+         "l_src": l_src,
+         "l_obj": l_obj,
+         "l_template": l_template}
 
     return d
 
@@ -459,7 +454,7 @@ def ninja_readme_build(path_module):
 # |_) o ._   _. ._
 # |_) | | | (_| | \/
 #                 /
-def get_program(path_module):
+def get_binaries(path_module):
     """
     Return the list of binaries (Path= namedtuple('Path', ['abs', 'rel']) for a module
     """
@@ -484,7 +479,7 @@ def get_program(path_module):
             return []
 
 
-def get_dict_binaries(mode="production"):
+def get_dict_binaries(l_module, mode="production"):
     """
     Return a dict [module] = list_binaries
     If a the production mode is enable only header module
@@ -497,20 +492,20 @@ def get_dict_binaries(mode="production"):
 
     # Create d_binaries
     # Ake module => binaries generated
-    for module in d_genealogy_path.keys():
-        l_binaries = get_program(module)
+    for module in l_module:
+        l_binaries = get_binaries(module)
 
         if l_binaries:
             d_binaries[module] += l_binaries
 
     if mode == "production":
 
-        dict_boss = get_dict_module_boss()
-        dict_boss_path = dict_module_genelogy_path(dict_boss)
+        dict_root = ModuleHandler.dict_root
+        dict_root_path = dict_module_genelogy_path(dict_root)
 
         d_binaries_condensed = defaultdict(list)
         for module in d_binaries:
-            d_binaries_condensed[dict_boss_path[module]] += d_binaries[module]
+            d_binaries_condensed[dict_root_path[module]] += d_binaries[module]
 
         d_binaries = d_binaries_condensed
 
@@ -649,11 +644,14 @@ if __name__ == "__main__":
     # G e n e a l o g y _ d i c t #
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
 
-    d_genealogy = get_dict_genealogy_desc()
+    d_genealogy = ModuleHandler.dict_descendant
     d_genealogy_path = dict_module_genelogy_path(d_genealogy)
     d_irp = get_file_dependency(d_genealogy_path)
 
-    d_binaries_production = get_dict_binaries(mode="production")
+    l_module = d_genealogy_path.keys()
+
+    d_binaries_production = get_dict_binaries(l_module,
+                                              mode="production")
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # M o d u l e _ t o _ i r p #
diff --git a/scripts/ezfio_interface/prepare_ezfio.sh b/scripts/ezfio_interface/prepare_ezfio.sh
deleted file mode 100755
index 7c4b6c2f..00000000
--- a/scripts/ezfio_interface/prepare_ezfio.sh
+++ /dev/null
@@ -1,18 +0,0 @@
-#!/bin/bash
-# Prepares the EZFIO library before running make
-# Wed Mar 25 21:25:11 CET 2015
-
-echo "Building EZFIO library"
-echo "=-=-=-=-=-=-=-=-=-=-=-"
-
-# Interpret all the EZFIO.cfg files
-for dir in ${QPACKAGE_ROOT}/src/*/ 
-do 
-  cd $dir || exit -1
-  ${QPACKAGE_ROOT}/scripts/ezfio_interface/ei_handler.py --ezfio_config --ezfio_default
-done
-
-# For old-style directories. Will be removed some day...
-cp ${QPACKAGE_ROOT}/src{/*,}/*.ezfio_config ${QPACKAGE_ROOT}/EZFIO/config 
-
-
diff --git a/scripts/module/build_modules.sh b/scripts/module/build_modules.sh
deleted file mode 100755
index 049fb2aa..00000000
--- a/scripts/module/build_modules.sh
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/bin/bash 
-#
-# This script should run from the directory $QPACKAGE_ROOT/src
-
-# Check is we are in `QPACKAGE_ROOT_SRC` and all the dependency are correct
-${QPACKAGE_ROOT}/scripts/module/check_src.sh
-
-NPROC=$(cat /proc/cpuinfo | grep MHz | wc -l)
-
-export IN_MAKE=1
-
-for MODULE in $@
-do
-  if [[ ! -d ${MODULE} ]]
-  then
-    error "Module ${MODULE} doesn't exist"
-  fi
-  cd ${MODULE}
-  echo ${MODULE}  
-
-  if [[ $# -eq 1 ]] 
-  then
-    env make -j ${NPROC} all 
-  else
-    env make -j ${NPROC} all &> make.log 
-    if [[ $? -ne 0 ]]
-    then
-      cat make.log
-      error "
-Build failed for module $MODULE
-"
-    fi
-  fi
-
-  cd ${OLDPWD}
-done
-${QPACKAGE_ROOT}/scripts/module/create_executables_list.sh
diff --git a/scripts/module/clean_modules.sh b/scripts/module/clean_modules.sh
index a11deb70..19db6f46 100755
--- a/scripts/module/clean_modules.sh
+++ b/scripts/module/clean_modules.sh
@@ -14,7 +14,7 @@ function do_clean()
 {
   rm -rf -- \
     IRPF90_temp IRPF90_man Makefile.depend \
-    $(module_handler.py print_genealogy) include \
+    $(module_handler.py print_descendant) include \
     ezfio_interface.irp.f irpf90.make irpf90_entities tags $(ls_exe) *.mod
 
   touch -c EZFIO.cfg *.ezfio_config
diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py
index 833083a6..f4ff006e 100755
--- a/scripts/module/module_handler.py
+++ b/scripts/module/module_handler.py
@@ -6,44 +6,31 @@ Create the NEEDED_MODULE
 of a NEEDED_CHILDREN_MODULES file
 
 Usage:
-    module_handler.py print_genealogy       [<NEEDED_CHILDREN_MODULES>]
+    module_handler.py print_descendant      [<NEEDED_CHILDREN_MODULES>]
     module_handler.py create_png            [<NEEDED_CHILDREN_MODULES>]
     module_handler.py head_module
 
 Options:
-    print_genealogy         Print the genealogy of the NEEDED_CHILDREN_MODULES
+    print_descendant         Print the genealogy of the NEEDED_CHILDREN_MODULES
                                  aka (children, subchildren, etc)
     create_png              Create a png of the file
     NEEDED_CHILDREN_MODULES The path of NEEDED_CHILDREN_MODULES
                                 by default try to open the file in the current path
 """
-
-from docopt import docopt
-
 import os
 import sys
 import os.path
-from cache import cache
 
-from collections import namedtuple
-Dependency = namedtuple('Dependency', ['src', 'obj'])
-Module_info = namedtuple('Module_info', ['l_children', 'l_dependency'])
-
-
-def get_list_from_makefile(data, sep):
-    # Split for sep
-    dump = [l.split(sep)[1] for l in data if l.startswith(sep)]
-
-    # Delete the empy one
-    l_unique = [k for k in map(str.strip, dump) if k]
-
-    # Return the flat one (if multi in l_unique)
-    l_flat = [j for i in l_unique for j in i.split()]
-    return l_flat
+try:
+    from docopt import docopt
+    from decorator import classproperty
+except ImportError:
+    print "source .quantum_package.rc"
+    raise
 
 
 # Canot cache for namedtuple are not hashable
-def get_dict_genealogy():
+def get_dict_child():
     """Loop over MODULE in  QPACKAGE_ROOT/src, open all the NEEDED_CHILDREN_MODULES
     and create a dict[MODULE] = [sub module needed, ...]
     """
@@ -59,23 +46,14 @@ def get_dict_genealogy():
             with open(path_file, "r") as f:
                 l_children = f.read().split()
         except IOError:
-            continue
-
-        try:
-            path_file = os.path.join(dir_, o, "Makefile")
-            with open(os.path.join(dir_, o, "Makefile"), "r") as f:
-                data = f.readlines()
-                l_depend = Dependency(get_list_from_makefile(data, "SRC="),
-                                      get_list_from_makefile(data, "OBJ="))
-        except IOError:
-            l_depend = []
-
-        d_ref[o] = Module_info(l_children, l_depend)
+            pass
+        else:
+            d_ref[o] = l_children
 
     return d_ref
 
 
-def him_and_all_children(d_ref, l_module):
+def l_module_generalogy_rec(d_chidlren, l_module):
     """
     From a list of module return the module and all of the genealogy
     """
@@ -85,7 +63,7 @@ def him_and_all_children(d_ref, l_module):
         if module not in l:
             l.append(module)
             try:
-                l.extend(him_and_all_children(d_ref, d_ref[module].l_children))
+                l.extend(l_module_generalogy_rec(d_chidlren, d_chidlren[module]))
             except KeyError:
                 print >> sys.stderr, "`{0}` not submodule".format(module)
                 print >> sys.stderr, "Check the corresponding NEEDED_CHILDREN_MODULES"
@@ -94,116 +72,116 @@ def him_and_all_children(d_ref, l_module):
     return list(set(l))
 
 
-def get_dict_genealogy_desc():
-    """
-    Get a dic of all the genealogy desc (children and all_children)
-    """
-    d_ref = get_dict_genealogy()
+class ModuleHandler:
 
-    d = {}
+    dict_child = get_dict_child()
 
-    for module_name in d_ref:
-        d[module_name] = him_and_all_children(d_ref,
-                                              d_ref[module_name].l_children)
+    @classproperty
+    def l_module(self):
+        return self.dict_child.keys()
 
-    return d
+    @classproperty
+    def dict_parent(self):
+        """
+        Get a dic of the first parent
+        """
+        d_child = self.dict_child
 
+        d = {}
 
-def get_dict_parent():
-    """
-    Get a dic of the first parent
-    """
-    d_ref = get_dict_genealogy()
+        for module_name in d_child:
+            d[module_name] = [i for i in d_child.keys() if module_name in d_child[i]]
 
-    d = {}
+        return d
 
-    for module_name in d_ref:
-        d[module_name] = [i for i in d_ref.keys()
-                          if module_name in d_ref[i].l_children]
+    @classproperty
+    def dict_descendant(self):
+        """
+        Get a dic of all the genealogy desc (children and all_children)
+        """
+        d = {}
 
-    return d
+        d_child = self.dict_child
 
+        for module_name in d_child:
+            d[module_name] = l_module_generalogy_rec(d_child,
+                                                     d_child[module_name])
 
-def get_dict_module_boss():
-    """
-    Return a dict(module_name) = module_boss
-    Module boss is a module who have not parent (a edge) and have module_name
-    in is genealogy
-    """
-    d_ref_asc = get_dict_parent()
-    d_ref_desc = get_dict_genealogy_desc()
+        return d
 
-    l_all_module = d_ref_asc.keys()
+    @classproperty
+    def dict_root(self):
+        """
+        Return a dict(module_name) = module_boss
+        The top node in a tree.
+        """
+        d_asc = self.dict_parent
+        d_desc = self.dict_descendant
 
-    d_module_boss = {}
+        l_all_module = self.l_module
 
-    for module in l_all_module:
-        d_module_boss[module] = [
-            p for p in l_all_module
-            if module in [p] + d_ref_desc[p] and not d_ref_asc[p]
-        ][0]
+        dict_root = {}
 
-    return d_module_boss
+        for module in l_all_module:
+            dict_root[module] = [ p for p in l_all_module if module in [p] + d_desc[p] and not d_asc[p]][0]
 
+        return dict_root
 
-def module_genealogy(module_name):
-    """
-    Take a name of a NEEDED_CHILDREN_MODULES
-    and return a list of all the {sub, subsub, ...}children
-    """
+    @classmethod
+    def l_descendant_unique(cls, l_module):
+        d_desc = cls.dict_descendant
 
-    d_ref = get_dict_genealogy()
-    return him_and_all_children(d_ref, d_ref[module_name].l_children)
+        d = {}
+        for module in l_module:
+            for e in d_desc[module]:
+                d[e] = 1
 
+        return d.keys()
 
-def file_dependency(module_name):
+    @classmethod
+    def l_reduce_tree(cls, l_module):
 
-    d_ref = get_dict_genealogy()
-    l_src, l_obj = d_ref[module_name].l_dependency
+        l_d_u = cls.l_descendant_unique(l_module)
+        l_module_reduce = []
+        for module in l_module:
+            if module not in l_d_u:
+                l_module_reduce.append(module)
 
-    l_children_module = him_and_all_children(d_ref,
-                                             d_ref[module_name].l_children)
-    for module in l_children_module:
-        l_src_dump, l_obj_dump = d_ref[module].l_dependency
-        l_src.extend("{0}/{1}".format(module, i) for i in l_src_dump)
-        l_obj.extend("IRPF90_temp/{0}/{1}".format(module, os.path.basename(i))
-                     for i in l_obj_dump)
+        return l_module_reduce
 
-    return Dependency(l_src, l_obj)
+    @classmethod
+    def create_png(cls, l_module):
+        """Create the png of the dependency tree for a l_module"""
 
+        # Init
+        import pydot
+        all_ready_done = []
 
-def create_png(l_module):
-    """Create the png of the dependency tree for a l_module"""
+        def draw_module_edge(module, l_children):
+            "Draw all the module recursifly"
 
-    # Init
-    import pydot
-    all_ready_done = []
+            if module not in all_ready_done:
+                for children in l_children:
+                    # Add Edge
+                    edge = pydot.Edge(module, children)
+                    graph.add_edge(edge)
+                    # Recurs
+                    draw_module_edge(children, d_ref[children])
+                all_ready_done.append(module)
 
-    def draw_module_edge(module, l_children):
-        "Draw all the module recursifly"
+        # Init
+        graph = pydot.Dot(graph_type='digraph')
+        d_ref = cls.dict_child
 
-        if module not in all_ready_done:
-            for children in l_children:
-                # Add Edge
-                edge = pydot.Edge(module, children)
-                graph.add_edge(edge)
-                # Recurs
-                draw_module_edge(children, d_ref[children].l_children)
-            all_ready_done.append(module)
+        # Create all the edge
+        for module in l_module:
+            node_a = pydot.Node(module, fontcolor="red")
+            graph.add_node(node_a)
+            draw_module_edge(module, d_ref[module])
 
-    # Init
-    graph = pydot.Dot(graph_type='digraph')
-    d_ref = get_dict_genealogy()
-
-    # Create all the edge
-    for module in l_module:
-        node_a = pydot.Node(module, fontcolor="red")
-        graph.add_node(node_a)
-        draw_module_edge(module, d_ref[module].l_children)
-
-    # Save
-    path = '{0}.png'.format("tree_dependency")
-    graph.write_png(path)
+        # Save
+        path = '{0}.png'.format("tree_dependency")
+        graph.write_png(path)
 
 
 if __name__ == '__main__':
@@ -220,13 +198,8 @@ if __name__ == '__main__':
 
     path_file = os.path.basename(dir_)
 
-    if arguments['print_genealogy']:
-        l_all_needed_molule = module_genealogy(path_file)
-        print " ".join(sorted(l_all_needed_molule))
+    if arguments['print_descendant']:
+        print " ".join(sorted(ModuleHandler.l_module))
 
     if arguments["create_png"]:
-        create_png([path_file])
-
-    if arguments["head_module"]:
-        for module, boss in get_dict_module_boss().iteritems():
-            print module, boss
+        ModuleHandler.create_png([path_file])
diff --git a/scripts/module/qp_create_module.py b/scripts/module/qp_create_module.py
new file mode 100755
index 00000000..923d88cd
--- /dev/null
+++ b/scripts/module/qp_create_module.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+Usage: qp_create_module.py list
+       qp_create_module.py create -n <name> [<children_module>...]
+"""
+
+import sys
+import os
+
+try:
+    from docopt import docopt
+    from module_handler import ModuleHandler
+    from update_README import Doc_key, Needed_key
+except ImportError:
+    print "source .quantum_package.rc"
+    raise
+
+
+def save_new_module(path, l_child):
+
+    # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
+    # N E E D E D _ C H I L D R E N _ M O D U L E S #
+    # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
+
+    try:
+        os.makedirs(path)
+    except OSError:
+        print "The module ({0}) already exist...".format(path)
+        sys.exit(1)
+
+    with open(os.path.join(path, "NEEDED_CHILDREN_MODULES"), "w") as f:
+        f.write(" ".join(l_child))
+        f.write("\n")
+
+    # ~#~#~#~#~#~#~ #
+    # R E A D _ M E #
+    # ~#~#~#~#~#~#~ #
+
+    module_name = os.path.basename(path)
+
+    header = "{0}\n{1}\n{0}\n".format("=" * len(module_name), module_name)
+
+    with open(os.path.join(path, "README.rst"), "w") as f:
+        f.write(header + "\n")
+        f.write(Doc_key + "\n")
+        f.write(Needed_key + "\n")
+
+
+if __name__ == '__main__':
+    arguments = docopt(__doc__)
+
+    if arguments["list"]:
+        for module in ModuleHandler.l_module:
+            print module
+
+    elif arguments["create"]:
+        l_children = arguments["<children_module>"]
+
+        qpackage_root = os.environ['QPACKAGE_ROOT']
+        path = os.path.join(qpackage_root, "src", arguments["<name>"])
+
+        print "You will create the module:"
+        print path
+
+        for children in l_children:
+            if children not in ModuleHandler.dict_descendant:
+                print "This module ({0}) is not a valide module.".format(children)
+                print "Run `list` flag for the list of module avalaible"
+                print "Aborting..."
+                sys.exit(1)
+
+        print "You ask for this submodule:"
+        print l_children
+
+        print "You can use all the routine in this module"
+        print l_children + ModuleHandler.l_descendant_unique(l_children)
+
+        print "This can be reduce to:"
+        l_child_reduce = ModuleHandler.l_reduce_tree(l_children)
+        print l_child_reduce
+        save_new_module(path, l_child_reduce)
diff --git a/scripts/module/qp_create_module.sh b/scripts/module/qp_create_module.sh
deleted file mode 100755
index 449bd5a2..00000000
--- a/scripts/module/qp_create_module.sh
+++ /dev/null
@@ -1,140 +0,0 @@
-#!/bin/bash
-# 
-# usage:
-#  create_module.sh MOs AOs Electrons
-# Prepares all the files for the creation of a new module.
-# The first argument is the name of the module
-# All remaining aruments are dependencies.
-# Thu Apr  3 01:44:58 CEST 2014
-
-if [[ -z ${QPACKAGE_ROOT} ]]
-then
-  print "The QPACKAGE_ROOT environment variable is not set."
-  print "Please reload the quantum_package.rc file."
-  exit -1
-fi
-source ${QPACKAGE_ROOT}/scripts/qp_include.sh
-
-check_current_dir_is_src
-
-DEBUG=0
-
-# If DEBUG=1, the print debug info.
-function debug()
-{
-  if [[ $DEBUG -eq 1 ]]
-  then
-    function debug()
-    {
-      echo "$@"
-    }
-  else
-    function debug()
-    {
-      :
-    }
-  fi
-  debug $@
-}
-
-function fail()
-{
-  echo "Error:  $@"
-  cd "${QPACKAGE_ROOT}/src"
-  rm -rf -- "${MODULE}"
-  exit 1
-}
-
-
-
-MODULE=$1
-
-# Check command line
-if [[ -z $MODULE ]]
-then
-  echo "usage: $(basename $0) <NewModuleName>"
-  exit 1
-fi
-
-shift
-
-
-
-# Check if module already exists
-if [ -d ${MODULE} ]
-then
-  echo "Error: Module $MODULE already exists"
-  exit 1
-fi
-
-debug "Module does not already exist: OK"
-
-
-# Set up dependencies
-ALL_MODULES="$(cat NEEDED_MODULES)"
-echo "Select which modules you are sure you will need:  (press q to quit)"  
-NEEDED_MODULES=""
-select M in ${ALL_MODULES}
-do
-  if [[ -z $M ]]
-  then
-    break
-  fi
-  NEEDED_MODULES+=" $M"
-  echo "$NEEDED_MODULES"
-done
-
-
-
-# Create module directory and go into it
-mkdir "${QPACKAGE_ROOT}/src/${MODULE}"
-if [[ $? != 0 ]]
-then
-  echo  "Error: Unable to create module directory."
-  exit 1
-fi
-
-if [[ ! -d  "${QPACKAGE_ROOT}/src/${MODULE}" ]]
-then
-  echo  "Something strange happened: the"
-  echo  "${QPACKAGE_ROOT}/src/${MODULE}"
-  echo  "directory was not created."
-  exit 1
-fi
-
-cd "${QPACKAGE_ROOT}/src/${MODULE}"
-if [[ "${PWD}" != "${QPACKAGE_ROOT}/src/${MODULE}" ]]
-then
-  echo  "Something strange happened: we should be in"
-  echo  "${QPACKAGE_ROOT}/src/${MODULE}"
-  echo  "but we are in"
-  echo  "${PWD}"
-  exit 1
-fi
-
-debug "Module directory is created."
-
-
-
-# Create the Makefile
-"${QPACKAGE_ROOT}/scripts/module/create_Makefile.sh" || fail "Unable to create Makefile"
-if [[ ! -f Makefile ]]
-then
-  fail "Makefile was not created"
-fi
-debug "Makefile created"
-
-# Create the NEEDED_MODULES file
-echo "$NEEDED_MODULES" > NEEDED_CHILDREN_MODULES
-debug "NEEDED_MODULES created"
-
-# Create rst templates
-"${QPACKAGE_ROOT}/scripts/module/create_rst_templates.sh" || fail "Unable to create rst templates"
-
-
-# Update module list in main NEEDED_MODULES
-ALL_MODULES+=" ${MODULE}"
-echo "${ALL_MODULES}" > "${QPACKAGE_ROOT}/src/NEEDED_MODULES"
-debug "Updated NEEDED_MODULES"
-
-
diff --git a/scripts/module/update_README.py b/scripts/module/update_README.py
index 7f4624b3..1293f6f2 100755
--- a/scripts/module/update_README.py
+++ b/scripts/module/update_README.py
@@ -23,19 +23,15 @@ header = """
 
 def fetch_splitted_data():
     """Read the README.rst file and split it in strings:
-    * The description
-    * The assumptions
     * The documentation
     * The needed modules
     The result is given as a list of strings
     """
 
-    file = open(README, 'r')
-    data = file.read()
-    file.close()
+    with open(README, 'r') as f:
+        data = f.read()
 
     # Place sentinels
-    data = data.replace(Assum_key, Sentinel + Assum_key)
     data = data.replace(Doc_key, Sentinel + Doc_key)
     data = data.replace(Needed_key, Sentinel + Needed_key)
 
diff --git a/scripts/utility/cache.py b/scripts/utility/cache.py
deleted file mode 100644
index 9be4ea2c..00000000
--- a/scripts/utility/cache.py
+++ /dev/null
@@ -1,17 +0,0 @@
-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/AOs/dimensions_integrals.irp.f b/src/AOs/dimensions_integrals.irp.f
index 86db70f6..19241819 100644
--- a/src/AOs/dimensions_integrals.irp.f
+++ b/src/AOs/dimensions_integrals.irp.f
@@ -3,13 +3,13 @@
  implicit none
  integer :: n_pt_sup
  integer :: prim_power_l_max
- include 'Utils/constants.F'
+ include 'Utils/constants.include.F'
  prim_power_l_max = maxval(ao_power)
  n_pt_max_integrals = 24 * prim_power_l_max + 4
  n_pt_max_i_x = 8 * prim_power_l_max
  ASSERT (n_pt_max_i_x-1 <= max_dim)
  if (n_pt_max_i_x-1 > max_dim) then
-   print *, 'Increase max_dim in Utils/constants.F to ', n_pt_max_i_x-1
+   print *, 'Increase max_dim in Utils/constants.include.F to ', n_pt_max_i_x-1
    stop 1
  endif
 END_PROVIDER
diff --git a/src/DensityFit/overlap.irp.f b/src/DensityFit/overlap.irp.f
index 24875c4b..6c671198 100644
--- a/src/DensityFit/overlap.irp.f
+++ b/src/DensityFit/overlap.irp.f
@@ -8,7 +8,7 @@ double precision function aux_basis_four_overlap(i,j,k,l)
   double precision               :: I_center(3),J_center(3),K_center(3),L_center(3)
   integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
   double precision               :: overlap_x,overlap_y,overlap_z, overlap
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
   double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
   integer                        :: iorder_p(3), iorder_q(3)
diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f
index a6ce3dc3..ba3bbcc1 100644
--- a/src/Integrals_Bielec/ao_bi_integrals.irp.f
+++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f
@@ -9,7 +9,7 @@ double precision function ao_bielec_integral(i,j,k,l)
   double precision               :: I_center(3),J_center(3),K_center(3),L_center(3)
   integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
   double precision               :: integral
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
   double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
   integer                        :: iorder_p(3), iorder_q(3)
@@ -115,7 +115,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
   double precision               :: I_center(3),J_center(3),K_center(3),L_center(3)
   integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
   double precision               :: integral
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
   double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
   integer                        :: iorder_p(3), iorder_q(3)
@@ -519,7 +519,7 @@ double precision function general_primitive_integral(dim,            &
   ! Computes the integral <pq|rs> where p,q,r,s are Gaussian primitives
   END_DOC
   integer,intent(in)             :: dim
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision, intent(in)   :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv
   double precision, intent(in)   :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv
   integer, intent(in)            :: iorder_p(3)
@@ -665,7 +665,7 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,
   integer                        :: n_pt_sup
   double precision               :: p,q,denom,coeff
   double precision               :: I_f
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   if(iand(a_x+b_x+c_x+d_x,1).eq.1.or.iand(a_y+b_y+c_y+d_y,1).eq.1.or.iand(a_z+b_z+c_z+d_z,1).eq.1)then
     ERI = 0.d0
     return
@@ -859,7 +859,7 @@ subroutine give_polynom_mult_center_x(P_center,Q_center,a_x,d_x,p,q,n_pt_in,pq_i
   integer, intent(in)            :: a_x,d_x
   double precision, intent(in)   :: P_center, Q_center
   double precision, intent(in)   :: p,q,pq_inv,p10_1,p01_1,p10_2,p01_2,pq_inv_2
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision,intent(out)   :: d(0:max_dim)
   double precision               :: accu
   accu = 0.d0
@@ -916,7 +916,7 @@ subroutine I_x1_pol_mult(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
   ! recursive function involved in the bielectronic integral
   END_DOC
   integer , intent(in)           :: n_pt_in
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision,intent(inout) :: d(0:max_dim)
   integer,intent(inout)          :: nd
   integer, intent(in)            :: a,c
@@ -950,7 +950,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
   ! recursive function involved in the bielectronic integral
   END_DOC
   integer , intent(in)           :: n_pt_in
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision,intent(inout) :: d(0:max_dim)
   integer,intent(inout)          :: nd
   integer, intent(in)            :: a,c
@@ -1036,7 +1036,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
   ! recursive function involved in the bielectronic integral
   END_DOC
   integer , intent(in)           :: n_pt_in
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision,intent(inout) :: d(0:max_dim)
   integer,intent(inout)          :: nd
   integer, intent(in)            :: c
@@ -1090,7 +1090,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
   ! recursive function involved in the bielectronic integral
   END_DOC
   integer , intent(in)           :: n_pt_in
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision,intent(inout) :: d(0:max_dim)
   integer,intent(inout)          :: nd
   integer, intent(in)            :: c
@@ -1152,7 +1152,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
   ! recursive function involved in the bielectronic integral
   END_DOC
   integer , intent(in)           :: dim
-  include 'Utils/constants.F'
+  include 'Utils/constants.include.F'
   double precision               :: d(0:max_dim)
   integer,intent(inout)          :: nd
   integer, intent(in)            :: c
diff --git a/src/Integrals_Monoelec/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f
index ea8e95a3..09b67425 100644
--- a/src/Integrals_Monoelec/pot_ao_ints.irp.f
+++ b/src/Integrals_Monoelec/pot_ao_ints.irp.f
@@ -148,7 +148,7 @@ double precision :: P_center(3)
 double precision :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor
 double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
 double precision ::  V_e_n,const_factor,dist_integral,tmp
-include 'Utils/constants.F'
+include 'Utils/constants.include.F'
   if ( (A_center(1)/=B_center(1)).or. &
        (A_center(2)/=B_center(2)).or. &
        (A_center(3)/=B_center(3)).or. &
@@ -351,7 +351,7 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
  integer,intent(inout) :: nd
  integer, intent(in):: a,c
  double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2)
- include 'Utils/constants.F'
+ include 'Utils/constants.include.F'
   double precision :: X(0:max_dim)
   double precision :: Y(0:max_dim)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
@@ -417,7 +417,7 @@ end
 recursive subroutine I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,d,nd,dim) 
  implicit none
  integer , intent(in) :: dim
- include 'Utils/constants.F'
+ include 'Utils/constants.include.F'
  double precision :: d(0:max_dim)
  integer,intent(inout) :: nd
  integer, intent(in):: c
@@ -492,7 +492,7 @@ implicit none
 double precision :: alpha
 integer :: n
 double precision :: dble_fact
-include 'Utils/constants.F'
+include 'Utils/constants.include.F'
 
 !if(iand(n,1).eq.1)then
 !  int_gaus_pol= 0.d0
@@ -521,7 +521,7 @@ double precision function V_r(n,alpha)
 implicit none
 double precision ::  alpha, fact
 integer :: n
-include 'Utils/constants.F'
+include 'Utils/constants.include.F'
 if(iand(n,1).eq.1)then
  V_r = 0.5d0 * fact(ishft(n,-1)) / (alpha ** (ishft(n,-1) + 1))
 else
@@ -549,7 +549,7 @@ implicit none
 !!     integral on "theta" with boundaries ( 0 ; pi) of [ cos(theta) **n  sin(theta) **m ]
 integer :: n,m,i
 double precision ::  Wallis, prod
-include 'Utils/constants.F'
+include 'Utils/constants.include.F'
   V_theta = 0.d0
   prod = 1.d0
   do i = 0,ishft(n,-1)-1
@@ -565,7 +565,7 @@ double precision function Wallis(n)
 implicit none
 double precision :: fact
 integer :: n,p
-include 'Utils/constants.F'
+include 'Utils/constants.include.F'
 if(iand(n,1).eq.0)then
   Wallis = fact(ishft(n,-1))
   Wallis = pi * fact(n) / (dble(ibset(0_8,n)) * (Wallis+Wallis)*Wallis) 
diff --git a/src/Utils/constants.F b/src/Utils/constants.F
deleted file mode 100644
index 632dc50b..00000000
--- a/src/Utils/constants.F
+++ /dev/null
@@ -1,10 +0,0 @@
-integer, parameter :: max_dim = 255
-integer, parameter :: SIMD_vector = 32
-
-double precision, parameter :: pi =  dacos(-1.d0)
-double precision, parameter :: sqpi =  dsqrt(dacos(-1.d0))
-double precision, parameter :: pi_5_2 =  34.9868366552d0
-double precision, parameter :: dfour_pi =  4.d0*dacos(-1.d0)
-double precision, parameter :: dtwo_pi =  2.d0*dacos(-1.d0)
-double precision, parameter :: inv_sq_pi =  1.d0/dsqrt(dacos(-1.d0))
-double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f
index a98efb3d..e9c3f9ab 100644
--- a/src/Utils/integration.irp.f
+++ b/src/Utils/integration.irp.f
@@ -6,7 +6,7 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph
   !        fact_k  (x-x_P)^iorder(1)  (y-y_P)^iorder(2)  (z-z_P)^iorder(3) exp(-p(r-P)^2)
   END_DOC
   implicit none
-  include 'constants.F'
+  include 'constants.include.F'
   integer, intent(in)            :: dim
   integer, intent(in)            :: a,b               ! powers : (x-xa)**a_x = (x-A(1))**a(1)
   double precision, intent(in)   :: alpha, beta       ! exponents
@@ -53,7 +53,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   !               * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
   END_DOC
   implicit none
-  include 'constants.F'
+  include 'constants.include.F'
   integer, intent(in)            :: dim
   integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
   double precision, intent(in)   :: alpha, beta       ! exponents
@@ -131,7 +131,7 @@ subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder
   !               * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
   END_DOC
   implicit none
-  include 'constants.F'
+  include 'constants.include.F'
   integer, intent(in)            :: dim
   integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
   double precision, intent(in)   :: alpha, beta, gama ! exponents
@@ -415,7 +415,7 @@ double precision function F_integral(n,p)
   double precision               :: p
   integer                        :: i,j
   double precision               :: accu,sqrt_p,fact_ratio,tmp,fact
-  include 'constants.F'
+  include 'constants.include.F'
   if(n < 0)then
     F_integral = 0.d0
   endif
@@ -441,7 +441,7 @@ double precision function rint(n,rho)
 !  \int_0^1 dx \exp(-p x^2) x^n
 !
   END_DOC
-  include 'constants.F'
+  include 'constants.include.F'
   double precision               :: rho,u,rint1,v,val0,rint_large_n,u_inv
   integer                        :: n,k
   double precision               :: two_rho_inv
@@ -486,7 +486,7 @@ double precision function rint_sum(n_pt_out,rho,d1)
   BEGIN_DOC
   ! Needed for the calculation of two-electron integrals.
   END_DOC
-  include 'constants.F'
+  include 'constants.include.F'
   integer, intent(in)            :: n_pt_out
   double precision, intent(in)   :: rho,d1(0:n_pt_out)
   double precision               :: u,rint1,v,val0,rint_large_n,u_inv
diff --git a/src/Utils/need.irp.f b/src/Utils/need.irp.f
index 22cb6a48..35b61c2e 100644
--- a/src/Utils/need.irp.f
+++ b/src/Utils/need.irp.f
@@ -46,7 +46,7 @@
 
       double precision function rinteg(n,u)
       implicit double precision(a-h,o-z)
-      include 'constants.F'
+      include 'constants.include.F'
 !     pi=dacos(-1.d0)
       ichange=1
       factor=1.d0
diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f
index 6602118f..fd15054e 100644
--- a/src/Utils/one_e_integration.irp.f
+++ b/src/Utils/one_e_integration.irp.f
@@ -6,7 +6,7 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   ! \sum_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2) dx
   !
   END_DOC
-  include 'constants.F'
+  include 'constants.include.F'
   integer,intent(in)             :: dim ! dimension maximum for the arrays representing the polynomials
   double precision,intent(in)    :: A_center,B_center  ! center of the x1 functions
   integer,intent(in)             :: power_A, power_B ! power of the x1 functions
@@ -34,7 +34,7 @@ end
 
 subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap)
  implicit none
-  include 'constants.F'
+  include 'constants.include.F'
  integer, intent(in)            :: dim                                                             
  integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
  double precision, intent(in)   :: alpha, beta, gama ! exponents
@@ -120,7 +120,7 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
   !   S = S_x S_y S_z
   !
   END_DOC
-  include 'constants.F'
+  include 'constants.include.F'
   integer,intent(in)             :: dim ! dimension maximum for the arrays representing the polynomials
   double precision,intent(in)    :: A_center(3),B_center(3)  ! center of the x1 functions
   double precision, intent(in)   :: alpha,beta
diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f
index 1a579662..02cd87f9 100644
--- a/src/Utils/util.irp.f
+++ b/src/Utils/util.irp.f
@@ -51,7 +51,7 @@ integer function align_double(n)
   ! Compute 1st dimension such that it is aligned for vectorization.
   END_DOC
   integer                        :: n
-  include 'constants.F'
+  include 'constants.include.F'
   if (mod(n,SIMD_vector/4) /= 0) then
     align_double= n + SIMD_vector/4 - mod(n,SIMD_vector/4)
   else