diff --git a/REPLACE b/REPLACE
index d0b6adce..1f9fff4b 100755
--- a/REPLACE
+++ b/REPLACE
@@ -463,7 +463,7 @@ qp_name ao_nucl_elec_integrals_per_atom -r ao_integrals_n_e_per_atom
 qp_name bi_elec_ref_bitmask_energy -r ref_bitmask_two_e_energy
 qp_name mono_elec_ref_bitmask_energy -r ref_bitmask_one_e_energy
 qp_name kinetic_ref_bitmask_energy -r ref_bitmask_kinetic_energy
-qp_name nucl_elec_ref_bitmask_energy -r ref_bitmask_e_n_energy
+qp_name nucl_elec_ref_bitmask_energy -r ref_bitmask_n_e_energy
 qp_name disk_access_ao_integrals_erf
 qp_name mo_bielec_integral_jj
 qp_name mo_bielec_integral_jj -r mo_two_e_integrals_jj
@@ -835,3 +835,7 @@ qp_name potential_sr_c_beta_ao_pbe --rename=potential_c_beta_ao_sr_pbe
 qp_name potential_sr_xc_alpha_ao_pbe --rename=potential_xc_alpha_ao_sr_pbe
 qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe
 qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion
+qp_name nucl_elec_ref_bitmask_energy -r ref_bitmask_n_e_energy
+qp_name ref_bitmask_e_n_energy -r ref_bitmask_n_e_energy
+qp_name read_ao_integrals_e_n -r read_ao_integrals_n_e
+qp_name write_ao_integrals_e_n -r write_ao_integrals_n_e
diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml
index 12654aad..631b6bbe 100644
--- a/ocaml/Input_mo_basis.ml
+++ b/ocaml/Input_mo_basis.ml
@@ -10,6 +10,7 @@ module Mo_basis : sig
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array;
         ao_md5          : MD5.t;
       } [@@deriving sexp]
   val read : unit -> t option
@@ -24,6 +25,7 @@ end = struct
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array;
         ao_md5          : MD5.t;
       } [@@deriving sexp]
   let get_default = Qpackage.get_ezfio_default "mo_basis"
@@ -37,11 +39,17 @@ end = struct
 
 
   let reorder b ordering =
-    { b with mo_coef =
+    { b with
+      mo_coef =
       Array.map (fun mo ->
         Array.init (Array.length mo)
           (fun i -> mo.(ordering.(i)))
-      ) b.mo_coef 
+      ) b.mo_coef  ;
+      mo_coef_imag =
+      Array.map (fun mo ->
+        Array.init (Array.length mo)
+          (fun i -> mo.(ordering.(i)))
+      ) b.mo_coef_imag 
     }
 
   let read_ao_md5 () =
@@ -121,6 +129,17 @@ end = struct
         Array.sub a (j*ao_num) (ao_num) 
       )
 
+  let read_mo_coef_imag () =
+    let a = Ezfio.get_mo_basis_mo_coef_imag ()
+            |> Ezfio.flattened_ezfio
+            |> Array.map MO_coef.of_float
+    in
+    let mo_num = read_mo_num () |> MO_number.to_int in
+    let ao_num = (Array.length a)/mo_num in
+    Array.init mo_num (fun j ->
+        Array.sub a (j*ao_num) (ao_num) 
+      )
+
 
   let read () =
     if (Ezfio.has_mo_basis_mo_num ()) then
@@ -130,6 +149,7 @@ end = struct
           mo_class        = read_mo_class ();
           mo_occ          = read_mo_occ ();
           mo_coef         = read_mo_coef ();
+          mo_coef_imag    = read_mo_coef_imag ();
           ao_md5          = read_ao_md5 ();
         }
     else
@@ -137,6 +157,7 @@ end = struct
 
 
   let mo_coef_to_string mo_coef =
+  (*TODO : add imaginary part here *)
     let ao_num = Array.length mo_coef.(0)
     and mo_num = Array.length mo_coef in
     let rec print_five imin imax =
@@ -222,6 +243,7 @@ MO coefficients ::
 
 
   let to_string b =
+  (*TODO : add imaginary part here *)
     Printf.sprintf "
 mo_label        = \"%s\"
 mo_num          = %s
@@ -281,7 +303,19 @@ mo_coef         = %s
     |> Array.to_list
     |> List.concat
     in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data
-    |> Ezfio.set_mo_basis_mo_coef
+    |> Ezfio.set_mo_basis_mo_coef;
+  ;;
+
+  let write_mo_coef_imag a =
+    let mo_num = Array.length a in
+    let ao_num = Array.length a.(0) in
+    let data =
+      Array.map (fun mo -> Array.map MO_coef.to_float mo
+      |> Array.to_list) a
+    |> Array.to_list
+    |> List.concat
+    in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data
+    |> Ezfio.set_mo_basis_mo_coef_imag;
   ;;
 
   let write 
@@ -290,6 +324,7 @@ mo_coef         = %s
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array;
         ao_md5          : MD5.t;
       } =
       write_mo_num mo_num;
@@ -297,6 +332,7 @@ mo_coef         = %s
       write_mo_class mo_class;
       write_mo_occ mo_occ;
       write_mo_coef mo_coef;
+      write_mo_coef_imag mo_coef_imag;
       write_md5 ao_md5
   ;;
 
diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg
index 930f7a77..9ef019fa 100644
--- a/src/ao_one_e_ints/EZFIO.cfg
+++ b/src/ao_one_e_ints/EZFIO.cfg
@@ -1,10 +1,16 @@
-[ao_integrals_e_n]
+[ao_integrals_n_e]
 type: double precision
 doc: Nucleus-electron integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
-[io_ao_integrals_e_n]
+[ao_integrals_n_e_imag]
+type: double precision
+doc: Imaginary part of the nucleus-electron integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
+[io_ao_integrals_n_e]
 type: Disk_access
 doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
@@ -17,6 +23,12 @@ doc: Kinetic energy integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_kinetic_imag]
+type: double precision
+doc: Imaginary part of the kinetic energy integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_kinetic]
 type: Disk_access
 doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ]
@@ -30,6 +42,12 @@ doc: Pseudopotential integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_pseudo_imag]
+type: double precision
+doc: Imaginary part of the pseudopotential integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_pseudo]
 type: Disk_access
 doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ]
@@ -43,6 +61,12 @@ doc: Overlap integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_overlap_imag]
+type: double precision
+doc: Imaginary part of the overlap integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_overlap]
 type: Disk_access
 doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ]
@@ -56,6 +80,12 @@ doc: Combined integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num) 
 interface: ezfio
 
+[ao_one_e_integrals_imag]
+type: double precision
+doc: Imaginary part of the combined integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num) 
+interface: ezfio
+
 [io_ao_one_e_integrals]
 type: Disk_access
 doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ]
diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f
index 1df7fd37..c3084ae2 100644
--- a/src/ao_one_e_ints/ao_one_e_ints.irp.f
+++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f
@@ -27,3 +27,24 @@
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)]
+  implicit none
+  integer :: i,j,n,l
+  BEGIN_DOC
+ ! One-electron Hamiltonian in the |AO| basis.
+  END_DOC
+
+  IF (read_ao_one_e_integrals) THEN
+     call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag)
+  ELSE
+     print *,  irp_here, ': Not yet implemented'
+     stop -1
+  ENDIF
+
+  IF (write_ao_one_e_integrals) THEN
+       call ezfio_set_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag)
+       print *,  'AO one-e integrals written to disk'
+  ENDIF
+
+END_PROVIDER
+
diff --git a/src/ao_one_e_ints/kin_ao_ints.irp.f b/src/ao_one_e_ints/kin_ao_ints.irp.f
index b87e2a22..4f117deb 100644
--- a/src/ao_one_e_ints/kin_ao_ints.irp.f
+++ b/src/ao_one_e_ints/kin_ao_ints.irp.f
@@ -149,3 +149,25 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)]
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)]
+  implicit none
+  BEGIN_DOC
+  ! Kinetic energy integrals in the |AO| basis.
+  !
+  ! $\langle \chi_i |\hat{T}| \chi_j \rangle$
+  !
+  END_DOC
+  integer                        :: i,j,k,l
+
+  if (read_ao_integrals_kinetic) then
+    call ezfio_get_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
+    print *,  'AO kinetic integrals read from disk'
+  else
+    print *,  irp_here, ': Not yet implemented'
+  endif
+  if (write_ao_integrals_kinetic) then
+    call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
+    print *,  'AO kinetic integrals written to disk'
+  endif
+END_PROVIDER
+
diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f
index bb23f3b1..1d4cab7d 100644
--- a/src/ao_one_e_ints/pot_ao_ints.irp.f
+++ b/src/ao_one_e_ints/pot_ao_ints.irp.f
@@ -12,8 +12,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
   integer                        :: i,j,k,l,n_pt_in,m
   double precision               :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
 
-  if (read_ao_integrals_e_n) then
-    call ezfio_get_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e)
+  if (read_ao_integrals_n_e) then
+    call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
     print *,  'AO N-e integrals read from disk'
   else
 
@@ -76,13 +76,36 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
     !$OMP END DO
     !$OMP END PARALLEL
   endif
-  if (write_ao_integrals_e_n) then
-    call ezfio_set_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e)
+  if (write_ao_integrals_n_e) then
+    call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
     print *,  'AO N-e integrals written to disk'
   endif
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)]
+  BEGIN_DOC
+  !  Nucleus-electron interaction, in the |AO| basis set.
+  !
+  !  :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
+  END_DOC
+  implicit none
+  double precision               :: alpha, beta, gama, delta
+  integer                        :: num_A,num_B
+  double precision               :: A_center(3),B_center(3),C_center(3)
+  integer                        :: power_A(3),power_B(3)
+  integer                        :: i,j,k,l,n_pt_in,m
+  double precision               :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
+
+  if (read_ao_integrals_n_e) then
+    call ezfio_get_ao_one_e_ints_ao_integrals_n_e_imag(ao_integrals_n_e_imag)
+    print *,  'AO N-e integrals read from disk'
+  else
+   print *,  irp_here, ': Not yet implemented'
+  endif
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)]
   BEGIN_DOC
 ! Nucleus-electron interaction in the |AO| basis set, per atom A.
@@ -166,7 +189,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
   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
+  double precision               :: V_n_e,const_factor,dist_integral,tmp
   double precision               :: accu,epsilo,rint
   integer                        :: n_pt_out,lmax
   include 'utils/constants.include.F'
@@ -178,7 +201,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
         (A_center(3)/=C_center(3))) then
     continue
   else
-    NAI_pol_mult = V_e_n(power_A(1),power_A(2),power_A(3),           &
+    NAI_pol_mult = V_n_e(power_A(1),power_A(2),power_A(3),           &
         power_B(1),power_B(2),power_B(3),alpha,beta)
     return
   endif
@@ -476,7 +499,7 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
   endif
 end
 
-double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
+double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
   implicit none
   BEGIN_DOC
 ! Primitve nuclear attraction between the two primitves centered on the same atom.
@@ -489,9 +512,9 @@ double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
   double precision               :: alpha,beta
   double precision               :: V_r, V_phi, V_theta
   if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then
-    V_e_n = 0.d0
+    V_n_e = 0.d0
   else
-    V_e_n =   V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta)              &
+    V_n_e =   V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta)              &
         * V_phi(a_x+b_x,a_y+b_y)                                     &
         * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1)
   endif
diff --git a/src/ao_one_e_ints_periodic/EZFIO.cfg b/src/ao_one_e_ints_periodic/EZFIO.cfg
deleted file mode 100644
index 5dac91f7..00000000
--- a/src/ao_one_e_ints_periodic/EZFIO.cfg
+++ /dev/null
@@ -1,72 +0,0 @@
-[ao_integrals_e_n]
-type: double precision
-doc: Nucleus-electron integrals in |AO| basis set
-size: (ao_basis.ao_num,ao_basis.ao_num)
-interface: ezfio
-
-[io_ao_integrals_e_n]
-type: Disk_access
-doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-
-
-[ao_integrals_kinetic]
-type: double precision
-doc: Kinetic energy integrals in |AO| basis set
-size: (ao_basis.ao_num,ao_basis.ao_num)
-interface: ezfio
-
-[io_ao_integrals_kinetic]
-type: Disk_access
-doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-
-
-[ao_integrals_pseudo]
-type: double precision
-doc: Pseudopotential integrals in |AO| basis set
-size: (ao_basis.ao_num,ao_basis.ao_num)
-interface: ezfio
-
-[io_ao_integrals_pseudo]
-type: Disk_access
-doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-
-
-[ao_integrals_overlap]
-type: double precision
-doc: Overlap integrals in |AO| basis set
-size: (ao_basis.ao_num,ao_basis.ao_num)
-interface: ezfio
-
-[io_ao_integrals_overlap]
-type: Disk_access
-doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-
-
-[ao_one_e_integrals]
-type: double precision
-doc: Combined integrals in |AO| basis set
-size: (ao_basis.ao_num,ao_basis.ao_num) 
-interface: ezfio
-
-[io_ao_one_e_integrals]
-type: Disk_access
-doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: Read
-
-
-[num_kpts]
-type: integer
-doc: Number of k-points
-interface: ezfio,provider,ocaml
-default: 1
-
-
diff --git a/src/ao_one_e_ints_periodic/NEED b/src/ao_one_e_ints_periodic/NEED
deleted file mode 100644
index 61d23b1e..00000000
--- a/src/ao_one_e_ints_periodic/NEED
+++ /dev/null
@@ -1,2 +0,0 @@
-ao_basis
-pseudo
diff --git a/src/ao_one_e_ints_periodic/README.rst b/src/ao_one_e_ints_periodic/README.rst
deleted file mode 100644
index 48ec8d01..00000000
--- a/src/ao_one_e_ints_periodic/README.rst
+++ /dev/null
@@ -1,15 +0,0 @@
-==================
-ao_one_e_integrals
-==================
-
-All the one-electron integrals in the periodic |AO| basis are here.
-
-Warning: this is incompatible with non-periodic |AOs|.
-
-The most important providers for usual quantum-chemistry calculation are:
-
-* `ao_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis 
-* `ao_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis
-* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis
-
-
diff --git a/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f b/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f
deleted file mode 100644
index 27e17754..00000000
--- a/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f
+++ /dev/null
@@ -1,29 +0,0 @@
- BEGIN_PROVIDER [ complex*16, ao_one_e_integrals,(ao_num,ao_num)]
-&BEGIN_PROVIDER [ double precision, ao_one_e_integrals_diag,(ao_num)]
-  implicit none
-  integer :: i,j,n,l
-  BEGIN_DOC
- ! One-electron Hamiltonian in the |AO| basis.
-  END_DOC
-
-  IF (read_ao_one_e_integrals) THEN
-     call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals)
-  ELSE
-        ao_one_e_integrals = ao_integrals_n_e + ao_kinetic_integrals
-
-        IF (DO_PSEUDO) THEN
-              ao_one_e_integrals += ao_pseudo_integrals
-        ENDIF
-  ENDIF
-
-  DO j = 1, ao_num
-    ao_one_e_integrals_diag(j) = real(ao_one_e_integrals(j,j))
-  ENDDO
-
-  IF (write_ao_one_e_integrals) THEN
-       call ezfio_set_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals)
-       print *,  'AO one-e integrals written to disk'
-  ENDIF
-
-END_PROVIDER
-
diff --git a/src/ao_one_e_ints_periodic/ao_overlap.irp.f b/src/ao_one_e_ints_periodic/ao_overlap.irp.f
deleted file mode 100644
index 9c74091f..00000000
--- a/src/ao_one_e_ints_periodic/ao_overlap.irp.f
+++ /dev/null
@@ -1,136 +0,0 @@
-BEGIN_PROVIDER [ complex*16, ao_overlap,(ao_num,ao_num) ]
-  implicit none
-  BEGIN_DOC  
-! Overlap between atomic basis functions:
-! :math:`\int \chi_i(r) \chi_j(r) dr)`
-  END_DOC
-  if (read_ao_integrals_overlap) then
-    call read_one_e_integrals_complex('ao_overlap', ao_overlap,&
-        size(ao_overlap,1), size(ao_overlap,2))
-    print *,  'AO overlap integrals read from disk'
-  else
-    print *, 'complex AO overlap integrals must be provided'
-  endif
-END_PROVIDER
-
-
-BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
-  implicit none
-  BEGIN_DOC  
-! Overlap between absolute value of atomic basis functions:
-! :math:`\int |\chi_i(r)| |\chi_j(r)| dr)`
-  END_DOC
-  integer :: i,j
-  !$OMP PARALLEL DO SCHEDULE(GUIDED) &
-  !$OMP DEFAULT(NONE) &
-  !$OMP PRIVATE(i,j) &
-  !$OMP SHARED(ao_overlap_abs,ao_overlap,ao_num)
-  do j=1,ao_num
-   do i= 1,ao_num
-    ao_overlap_abs(i,j)= cdabs(ao_overlap(i,j))
-   enddo
-  enddo
-  !$OMP END PARALLEL DO
-END_PROVIDER
-
-BEGIN_PROVIDER [ complex*16, S_inv,(ao_num,ao_num) ]
- implicit none
- BEGIN_DOC
-! Inverse of the overlap matrix
- END_DOC
- call get_pseudo_inverse_complex(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1))
-END_PROVIDER
-
-BEGIN_PROVIDER [ complex*16, S_half_inv, (AO_num,AO_num) ]
-
-  BEGIN_DOC
-!   :math:`X = S^{-1/2}` obtained by SVD
-  END_DOC
-
-  implicit none
-
-  integer                         :: num_linear_dependencies
-  integer                         :: LDA, LDC
-  double precision, allocatable   :: D(:)
-  complex*16, allocatable   :: U(:,:),Vt(:,:)
-  integer                         :: info, i, j, k
-  double precision, parameter     :: threshold_overlap_AO_eigenvalues = 1.d-6
-
-  LDA = size(AO_overlap,1)
-  LDC = size(S_half_inv,1)
-
-  allocate(         &
-    U(LDC,AO_num),  &
-    Vt(LDA,AO_num), &
-    D(AO_num))
-
-  call svd_complex(    &
-       AO_overlap,LDA, &
-       U,LDC,          &
-       D,              &
-       Vt,LDA,         &
-       AO_num,AO_num)
-
-  num_linear_dependencies = 0
-  do i=1,AO_num
-    print*,D(i)
-    if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
-      D(i) = 0.d0
-      num_linear_dependencies += 1
-    else
-      ASSERT (D(i) > 0.d0)
-      D(i) = 1.d0/sqrt(D(i))
-    endif
-    do j=1,AO_num
-      S_half_inv(j,i) = (0.d0,0.d0)
-    enddo
-  enddo
-  write(*,*) 'linear dependencies',num_linear_dependencies
-
-  do k=1,AO_num
-    if(D(k) /= 0.d0) then
-      do j=1,AO_num
-        do i=1,AO_num
-          S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j)
-        enddo
-      enddo
-    endif
-  enddo
-
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ complex*16, S_half, (ao_num,ao_num)  ]
- implicit none
- BEGIN_DOC
- ! :math:`S^{1/2}`
- END_DOC
-
-  integer :: i,j,k
-  complex*16, allocatable  :: U(:,:)
-  complex*16, allocatable  :: Vt(:,:)
-  double precision, allocatable  :: D(:)
-
-  allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num))
-
-  call svd_complex(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num)
-
-  do i=1,ao_num
-    D(i) = dsqrt(D(i))
-    do j=1,ao_num
-      S_half(j,i) = (0.d0,0.d0)
-    enddo
-  enddo
-
-  do k=1,ao_num
-      do j=1,ao_num
-        do i=1,ao_num
-          S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
-        enddo
-      enddo
-  enddo
-
-  deallocate(U,Vt,D)
-
-END_PROVIDER
-
diff --git a/src/ao_one_e_ints_periodic/aos_kpts.irp.f b/src/ao_one_e_ints_periodic/aos_kpts.irp.f
deleted file mode 100644
index e44d3a23..00000000
--- a/src/ao_one_e_ints_periodic/aos_kpts.irp.f
+++ /dev/null
@@ -1,5 +0,0 @@
-
-BEGIN_PROVIDER [integer, ao_num_per_kpt ]
-  implicit none
-  ao_num_per_kpt = ao_num / num_kpts
-END_PROVIDER
diff --git a/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f b/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f
deleted file mode 100644
index ea84bd4f..00000000
--- a/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f
+++ /dev/null
@@ -1,18 +0,0 @@
-BEGIN_PROVIDER [complex*16, ao_kinetic_integrals, (ao_num,ao_num)]
-  implicit none
-  BEGIN_DOC
-  ! array of the priminitve basis kinetic integrals
-  !  \langle \chi_i |\hat{T}| \chi_j \rangle
-  END_DOC
-  
-  if (read_ao_integrals_kinetic) then
-    call read_one_e_integrals_complex('ao_kinetic_integrals', ao_kinetic_integrals,&
-        size(ao_kinetic_integrals,1), size(ao_kinetic_integrals,2))
-    print *,  'AO kinetic integrals read from disk'
-  else
-    print *, 'complex AO kinetic integrals must be provided'
-    stop
-  endif
-END_PROVIDER
-
-
diff --git a/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f b/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f
deleted file mode 100644
index 36b7f7ae..00000000
--- a/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f
+++ /dev/null
@@ -1,28 +0,0 @@
-BEGIN_PROVIDER [ complex*16, ao_integrals_n_e, (ao_num,ao_num)]
-  BEGIN_DOC
-  !  Nucleus-electron interaction, in the |AO| basis set.
-  !
-  !  :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
-  END_DOC
-
-   if (read_ao_integrals_e_n) then
-    call read_one_e_integrals_complex('ao_ne_integral', ao_integrals_n_e,      &
-            size(ao_integrals_n_e,1), size(ao_integrals_n_e,2))
-     print *,  'AO N-e integrals read from disk'
-   else
-     print *, 'complex AO N-e integrals must be provided'
-     stop
-   endif
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)]
-  BEGIN_DOC
-! Nucleus-electron interaction in the |AO| basis set, per atom A.
-!
-! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle`
-  END_DOC
-  print *, 'ao_nucl_elec_integral_per_atom not implemented for k-points'
-  stop
-
-END_PROVIDER
-
diff --git a/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f
deleted file mode 100644
index 243efdcb..00000000
--- a/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f
+++ /dev/null
@@ -1,25 +0,0 @@
-BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
-  implicit none
-  BEGIN_DOC
-  ! Pseudo-potential integrals in the |AO| basis set.
-  END_DOC
-
-  if (read_ao_integrals_pseudo) then
-    call ezfio_get_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals)
-    print *,  'AO pseudopotential integrals read from disk'
-  else
-
-    if (do_pseudo) then
-      print *, irp_here, 'Not yet implemented'
-      stop -1
-    endif
-  endif
-
-  if (write_ao_integrals_pseudo) then
-    call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals)
-    print *,  'AO pseudopotential integrals written to disk'
-  endif
-
-END_PROVIDER
-
-
diff --git a/src/ao_one_e_ints_periodic/test.irp.f b/src/ao_one_e_ints_periodic/test.irp.f
deleted file mode 100644
index 3519c17f..00000000
--- a/src/ao_one_e_ints_periodic/test.irp.f
+++ /dev/null
@@ -1,3 +0,0 @@
-program test
- print *, ' hello'
-end
diff --git a/src/mo_basis/track_orb.irp.f b/src/bitmask/track_orb.irp.f
similarity index 97%
rename from src/mo_basis/track_orb.irp.f
rename to src/bitmask/track_orb.irp.f
index c49c545b..1cdde9cb 100644
--- a/src/mo_basis/track_orb.irp.f
+++ b/src/bitmask/track_orb.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, mo_coef_begin_iteration, (ao_num,mo_num) ]
    BEGIN_DOC
    ! Void provider to store the coefficients of the |MO| basis at the beginning of the SCF iteration
    !
-   ! Usefull to track some orbitals
+   ! Useful to track some orbitals
    END_DOC
 END_PROVIDER
 
diff --git a/src/determinants/ref_bitmask.irp.f b/src/determinants/ref_bitmask.irp.f
index ef6af40e..4e029ceb 100644
--- a/src/determinants/ref_bitmask.irp.f
+++ b/src/determinants/ref_bitmask.irp.f
@@ -1,7 +1,7 @@
  BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]
 &BEGIN_PROVIDER [ double precision, ref_bitmask_one_e_energy ]
 &BEGIN_PROVIDER [ double precision, ref_bitmask_kinetic_energy ]
-&BEGIN_PROVIDER [ double precision, ref_bitmask_e_n_energy ]
+&BEGIN_PROVIDER [ double precision, ref_bitmask_n_e_energy ]
 &BEGIN_PROVIDER [ double precision, ref_bitmask_two_e_energy ]
 &BEGIN_PROVIDER [ double precision, ref_bitmask_energy_ab ]
 &BEGIN_PROVIDER [ double precision, ref_bitmask_energy_bb ]
@@ -23,19 +23,19 @@
   ref_bitmask_energy = 0.d0
   ref_bitmask_one_e_energy = 0.d0
   ref_bitmask_kinetic_energy   = 0.d0
-  ref_bitmask_e_n_energy = 0.d0
+  ref_bitmask_n_e_energy = 0.d0
   ref_bitmask_two_e_energy = 0.d0
 
   do i = 1, elec_beta_num
     ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1)) + mo_one_e_integrals(occ(i,2),occ(i,2))
     ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1)) + mo_kinetic_integrals(occ(i,2),occ(i,2))
-    ref_bitmask_e_n_energy += mo_integrals_n_e(occ(i,1),occ(i,1)) + mo_integrals_n_e(occ(i,2),occ(i,2))
+    ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1)) + mo_integrals_n_e(occ(i,2),occ(i,2))
   enddo
 
   do i = elec_beta_num+1,elec_alpha_num
     ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1))
     ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1))
-    ref_bitmask_e_n_energy += mo_integrals_n_e(occ(i,1),occ(i,1))
+    ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1))
   enddo
 
   do j= 1, elec_alpha_num
@@ -55,7 +55,7 @@
       ref_bitmask_energy += mo_two_e_integrals_jj(occ(i,1),occ(j,2))
     enddo
   enddo
-  ref_bitmask_one_e_energy = ref_bitmask_kinetic_energy +   ref_bitmask_e_n_energy
+  ref_bitmask_one_e_energy = ref_bitmask_kinetic_energy +   ref_bitmask_n_e_energy
 
  ref_bitmask_energy_ab = 0.d0
  do i = 1, elec_alpha_num
diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg
index a055aad3..874af46a 100644
--- a/src/mo_basis/EZFIO.cfg
+++ b/src/mo_basis/EZFIO.cfg
@@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO|
 interface: ezfio
 size: (ao_basis.ao_num,mo_basis.mo_num)
 
+[mo_coef_imag]
+type: double precision
+doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)
+
 [mo_label]
 type: character*(64)
 doc: Label characterizing the MOS (Local, Canonical, Natural, *etc*)
diff --git a/src/mo_basis/NEED b/src/mo_basis/NEED
index 1a6a943c..5a504c20 100644
--- a/src/mo_basis/NEED
+++ b/src/mo_basis/NEED
@@ -1,3 +1,3 @@
 ao_basis
-ao_one_e_ints
 electrons
+ao_one_e_ints
diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f
index aa04fb01..73d33901 100644
--- a/src/mo_basis/mos.irp.f
+++ b/src/mo_basis/mos.irp.f
@@ -93,6 +93,59 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ]
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, mo_coef_imag, (ao_num,mo_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Molecular orbital coefficients on |AO| basis set
+  !
+  ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO|
+  !
+  ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc)
+  END_DOC
+  integer                        :: i, j
+  double precision, allocatable  :: buffer(:,:)
+  logical                        :: exists
+  PROVIDE ezfio_filename
+
+
+  if (mpi_master) then
+    ! Coefs
+    call ezfio_has_mo_basis_mo_coef_imag(exists)
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_coef_imag with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag)
+      write(*,*) 'Read  mo_coef_imag'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_coef_imag with MPI'
+      endif
+    IRP_ENDIF
+  else
+    ! Orthonormalized AO basis
+    do i=1,mo_num
+      do j=1,ao_num
+        mo_coef_imag(j,i) = 0.d0
+      enddo
+    enddo
+  endif
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ]
  implicit none
  BEGIN_DOC
diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg
index ebdcfead..0f9b63eb 100644
--- a/src/nuclei/EZFIO.cfg
+++ b/src/nuclei/EZFIO.cfg
@@ -31,3 +31,9 @@ default: None
 doc: Nuclear repulsion (Computed automaticaly or Read in the |EZFIO|)
 type:double precision
 interface: ezfio
+
+[periodic]
+type: logical
+doc: If true, the calculation uses periodic boundary conditions
+interface: ezfio, provider, ocaml
+default: false