diff --git a/scripts/module/create_Makefile_depend.sh b/scripts/module/create_Makefile_depend.sh
index 56e1c5e0..2485bb97 100755
--- a/scripts/module/create_Makefile_depend.sh
+++ b/scripts/module/create_Makefile_depend.sh
@@ -23,12 +23,12 @@ DEPS="$NEEDED_MODULES"
 for M in ${DEPS}
 do
   # X is the list of external source files
-  X=$(grep 'SRC=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2) 
+  X=$(grep '^SRC=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2) 
   for f in ${X}
   do
     SRC+=" ${M}/${f}"
   done
-  X=$(grep 'OBJ=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2) 
+  X=$(grep '^OBJ=' "${QPACKAGE_ROOT}/src/${M}/Makefile" 2>/dev/null |cut -d '=' -f 2) 
   for f in ${X}
   do
     OBJ+=" IRPF90_temp/${M}/${f/IRPF90_temp//}"
diff --git a/src/AOs/dimensions_integrals.irp.f b/src/AOs/dimensions_integrals.irp.f
index 635e1a0e..d1beea61 100644
--- a/src/AOs/dimensions_integrals.irp.f
+++ b/src/AOs/dimensions_integrals.irp.f
@@ -3,7 +3,7 @@
  implicit none
  integer :: n_pt_sup
  integer :: prim_power_l_max
- include 'constants.F'
+ include 'include/constants.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
diff --git a/src/Determinants/Makefile b/src/Determinants/Makefile
index 092d879d..51375558 100644
--- a/src/Determinants/Makefile
+++ b/src/Determinants/Makefile
@@ -1,6 +1,7 @@
 # Define here all new external source files and objects.Don't forget to prefix the
 # object files with IRPF90_temp/
-SRC=H_apply_template.f
+#SRC=H_apply_template.f
+SRC=
 OBJ=
 
 include $(QPACKAGE_ROOT)/src/Makefile.common
diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f
index a9fc3eb1..e213c2b9 100644
--- a/src/Integrals_Bielec/ao_bi_integrals.irp.f
+++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f
@@ -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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 'constants.F'
+  include 'include/constants.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 1e84d3d4..300bf99f 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 'constants.F'
+include 'include/constants.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 'constants.F'
+ include 'include/constants.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 'constants.F'
+ include 'include/constants.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 'constants.F'
+include 'include/constants.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 'constants.F'
+include 'include/constants.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 'constants.F'
+include 'include/constants.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 'constants.F'
+include 'include/constants.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/Makefile.common b/src/Makefile.common
index 606e6f7c..d6c9dd90 100644
--- a/src/Makefile.common
+++ b/src/Makefile.common
@@ -1,5 +1,5 @@
 .PHONY: default silent
-
+export
 ifneq ($(IN_MAKE),1) 
 
 default:
diff --git a/src/Perturbation/Makefile b/src/Perturbation/Makefile
index b668f388..fe655ac1 100644
--- a/src/Perturbation/Makefile
+++ b/src/Perturbation/Makefile
@@ -1,6 +1,7 @@
 # Define here all new external source files and objects.Don't forget to prefix the
 # object files with IRPF90_temp/
-SRC=perturbation_template.f
+SRC=
+#SRC=perturbation_template.f
 OBJ=
 
 include $(QPACKAGE_ROOT)/src/Makefile.common
diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f
index 89dc559d..ddeee2a9 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 'include/constants.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 'include/constants.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 'include/constants.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
diff --git a/src/Utils/need.irp.f b/src/Utils/need.irp.f
index 22cb6a48..f34dbad6 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 '../include/constants.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 d356454b..591e07a5 100644
--- a/src/Utils/one_e_integration.irp.f
+++ b/src/Utils/one_e_integration.irp.f
@@ -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 'include/constants.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