diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f
index 7fcb980a..1b1595a3 100644
--- a/src/ao_basis/aos_in_r.irp.f
+++ b/src/ao_basis/aos_in_r.irp.f
@@ -65,46 +65,60 @@ double precision function primitive_value(i,j,r)
 
 end
 
+! ---
 
-subroutine give_all_aos_at_r(r,aos_array)
- implicit none
- BEGIN_dOC
-! input  : r == r(1) = x and so on
-!
-! output : aos_array(i) = aos(i) evaluated in $\textbf{r}$
- END_DOC
- double precision, intent(in) :: r(3)
- double precision, intent(out):: aos_array(ao_num)
+subroutine give_all_aos_at_r(r, tmp_array)
 
- integer :: power_ao(3)
- integer :: i,j,k,l,m
- double precision :: dx,dy,dz,r2
- double precision ::      dx2,dy2,dz2
- double precision :: center_ao(3)
- double precision :: beta
- do i = 1, nucl_num
-  center_ao(1:3) = nucl_coord(i,1:3)
-  dx = (r(1) - center_ao(1))
-  dy = (r(2) - center_ao(2))
-  dz = (r(3) - center_ao(3))
-  r2 = dx*dx + dy*dy + dz*dz
-  do j = 1,Nucl_N_Aos(i)
-   k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
-   aos_array(k) = 0.d0
-   power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
-   dx2 = dx**power_ao(1)
-   dy2 = dy**power_ao(2)
-   dz2 = dz**power_ao(3)
-   do l = 1,ao_prim_num(k)
-    beta = ao_expo_ordered_transp_per_nucl(l,j,i)
-    if(dabs(beta*r2).gt.40.d0)cycle
-    aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
-   enddo
-   aos_array(k) = aos_array(k) * dx2 * dy2 * dz2
+  BEGIN_dOC
+  !
+  ! input  : r == r(1) = x and so on
+  !
+  ! output : tmp_array(i) = aos(i) evaluated in $\textbf{r}$
+  !
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: tmp_array(ao_num)
+  integer                       :: p_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx2, dy2, dz2
+  double precision              :: c_ao(3)
+  double precision              :: beta
+
+  do i = 1, nucl_num
+
+    c_ao(1:3) = nucl_coord(i,1:3)
+    dx = r(1) - c_ao(1)
+    dy = r(2) - c_ao(2)
+    dz = r(3) - c_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+
+    do j = 1, Nucl_N_Aos(i)
+
+      k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
+      p_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i)
+      dx2 = dx**p_ao(1)
+      dy2 = dy**p_ao(2)
+      dz2 = dz**p_ao(3)
+
+      tmp_array(k) = 0.d0
+      do l = 1,ao_prim_num(k)
+        beta = ao_expo_ordered_transp_per_nucl(l,j,i)
+        if(dabs(beta*r2).gt.40.d0) cycle
+
+        tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+      enddo
+
+      tmp_array(k) = tmp_array(k) * dx2 * dy2 * dz2
+    enddo
   enddo
- enddo
+
+  return
 end
 
+! ---
 
 subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
  implicit none
diff --git a/src/ao_basis/aos_transp.irp.f b/src/ao_basis/aos_transp.irp.f
index ae6193bf..4e44a9f6 100644
--- a/src/ao_basis/aos_transp.irp.f
+++ b/src/ao_basis/aos_transp.irp.f
@@ -1,20 +1,28 @@
- BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
- implicit none
- BEGIN_DOC
- ! List of AOs attached on each atom
- END_DOC
- integer :: i
- integer, allocatable :: nucl_tmp(:)
- allocate(nucl_tmp(nucl_num))
- nucl_tmp = 0
- Nucl_Aos = 0
- do i = 1, ao_num
-  nucl_tmp(ao_nucl(i))+=1
-  Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
- enddo
- deallocate(nucl_tmp)
+
+! ---
+
+BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
+
+  BEGIN_DOC
+  ! List of AOs attached on each atom
+  END_DOC
+ 
+  implicit none
+  integer              :: i
+  integer, allocatable :: nucl_tmp(:)
+
+  allocate(nucl_tmp(nucl_num))
+  nucl_tmp = 0
+  do i = 1, ao_num
+    nucl_tmp(ao_nucl(i)) += 1
+    Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
+  enddo
+  deallocate(nucl_tmp)
+
 END_PROVIDER
 
+! ---
+
 BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ]
  implicit none
  integer :: i,j,k,l
diff --git a/src/bi_ort_ints/one_e_bi_ort.irp.f b/src/bi_ort_ints/one_e_bi_ort.irp.f
index 5f2795f1..49181182 100644
--- a/src/bi_ort_ints/one_e_bi_ort.irp.f
+++ b/src/bi_ort_ints/one_e_bi_ort.irp.f
@@ -6,7 +6,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
   implicit none
   integer :: i, j
 
-  ao_one_e_integrals_tc_tot = ao_one_e_integrals      
+  ao_one_e_integrals_tc_tot = ao_one_e_integrals
 
   !provide j1b_type
 
diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f
index 721ea0c8..c1bacbd0 100644
--- a/src/bi_ort_ints/total_twoe_pot.irp.f
+++ b/src/bi_ort_ints/total_twoe_pot.irp.f
@@ -80,6 +80,8 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
       enddo
     enddo
 
+    FREE ao_tc_int_chemist
+
   endif
 
 END_PROVIDER 
@@ -128,69 +130,99 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
 
   implicit none
   integer                       :: i, j, k, l, m, n, p, q
-  double precision, allocatable :: mo_tmp_1(:,:,:,:), mo_tmp_2(:,:,:,:)
+  double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:)
 
-  allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num))
-  mo_tmp_1 = 0.d0
+  PROVIDE mo_r_coef mo_l_coef
 
-  do m = 1, ao_num
-    do p = 1, ao_num
-      do n = 1, ao_num
-        do q = 1, ao_num
-          do k = 1, mo_num
-            !       (k n|p m)    = sum_q c_qk * (q n|p m)
-            mo_tmp_1(k,n,p,m) += mo_l_coef_transp(k,q) * ao_two_e_tc_tot(q,n,p,m)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
+  allocate(a2(ao_num,ao_num,ao_num,mo_num))
 
-  allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num))
-  mo_tmp_2 = 0.d0
+  call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0     &
+            , ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
+            , 0.d0 , a2(1,1,1,1), ao_num*ao_num*ao_num)
 
-  do m = 1, ao_num
-    do p = 1, ao_num
-      do n = 1, ao_num
-        do i = 1, mo_num
-          do k = 1, mo_num
-            !       (k i|p m) = sum_n c_ni * (k n|p m)
-            mo_tmp_2(k,i,p,m) += mo_r_coef_transp(i,n) * mo_tmp_1(k,n,p,m)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  deallocate(mo_tmp_1)
+  allocate(a1(ao_num,ao_num,mo_num,mo_num))
 
-  allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num))
-  mo_tmp_1 = 0.d0
-  do m = 1, ao_num
-    do p = 1, ao_num
-      do l = 1, mo_num
-        do i = 1, mo_num
-          do k = 1, mo_num
-            mo_tmp_1(k,i,l,m) += mo_l_coef_transp(l,p) * mo_tmp_2(k,i,p,m)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  deallocate(mo_tmp_2)
+  call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 &
+            , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num          &
+            , 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num)
 
-  mo_bi_ortho_tc_two_e_chemist = 0.d0 
-  do m = 1, ao_num
-    do j = 1, mo_num
-      do l = 1, mo_num
-        do i = 1, mo_num
-          do k = 1, mo_num
-            mo_bi_ortho_tc_two_e_chemist(k,i,l,j) += mo_r_coef_transp(j,m) * mo_tmp_1(k,i,l,m)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  deallocate(mo_tmp_1)
+  deallocate(a2)
+  allocate(a2(ao_num,mo_num,mo_num,mo_num))
+
+  call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
+            , a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num          &
+            , 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num)
+
+  deallocate(a1)
+
+  call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
+            , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num          &
+            , 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num)
+
+  deallocate(a2)
+
+
+  !allocate(a1(mo_num,ao_num,ao_num,ao_num))
+  !a1 = 0.d0
+
+  !do m = 1, ao_num
+  !  do p = 1, ao_num
+  !    do n = 1, ao_num
+  !      do q = 1, ao_num
+  !        do k = 1, mo_num
+  !          !       (k n|p m)    = sum_q c_qk * (q n|p m)
+  !          a1(k,n,p,m) += mo_l_coef_transp(k,q) * ao_two_e_tc_tot(q,n,p,m)
+  !        enddo
+  !      enddo
+  !    enddo
+  !  enddo
+  !enddo
+
+  !allocate(a2(mo_num,mo_num,ao_num,ao_num))
+  !a2 = 0.d0
+
+  !do m = 1, ao_num
+  !  do p = 1, ao_num
+  !    do n = 1, ao_num
+  !      do i = 1, mo_num
+  !        do k = 1, mo_num
+  !          !       (k i|p m) = sum_n c_ni * (k n|p m)
+  !          a2(k,i,p,m) += mo_r_coef_transp(i,n) * a1(k,n,p,m)
+  !        enddo
+  !      enddo
+  !    enddo
+  !  enddo
+  !enddo
+  !deallocate(a1)
+
+  !allocate(a1(mo_num,mo_num,mo_num,ao_num))
+  !a1 = 0.d0
+  !do m = 1, ao_num
+  !  do p = 1, ao_num
+  !    do l = 1, mo_num
+  !      do i = 1, mo_num
+  !        do k = 1, mo_num
+  !          a1(k,i,l,m) += mo_l_coef_transp(l,p) * a2(k,i,p,m)
+  !        enddo
+  !      enddo
+  !    enddo
+  !  enddo
+  !enddo
+  !deallocate(a2)
+
+  !mo_bi_ortho_tc_two_e_chemist = 0.d0 
+  !do m = 1, ao_num
+  !  do j = 1, mo_num
+  !    do l = 1, mo_num
+  !      do i = 1, mo_num
+  !        do k = 1, mo_num
+  !          mo_bi_ortho_tc_two_e_chemist(k,i,l,j) += mo_r_coef_transp(j,m) * a1(k,i,l,m)
+  !        enddo
+  !      enddo
+  !    enddo
+  !  enddo
+  !enddo
+  !deallocate(a1)
 
 END_PROVIDER 
 
@@ -209,6 +241,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e, (mo_num, mo_num, mo_num,
   implicit none
   integer :: i, j, k, l
 
+  PROVIDE mo_bi_ortho_tc_two_e_chemist
+
   do j = 1, mo_num
     do i = 1, mo_num
       do l = 1, mo_num
@@ -220,6 +254,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e, (mo_num, mo_num, mo_num,
     enddo
   enddo
 
+  FREE mo_bi_ortho_tc_two_e_chemist
+
 END_PROVIDER 
 
 ! ---
diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f
index c69309d1..13eedfb7 100644
--- a/src/bi_ortho_mos/mos_rl.irp.f
+++ b/src/bi_ortho_mos/mos_rl.irp.f
@@ -136,6 +136,7 @@ BEGIN_PROVIDER [ double precision, mo_r_coef, (ao_num, mo_num) ]
         mo_r_coef(j,i) = mo_coef(j,i)
       enddo
     enddo
+    call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
   endif
 
 END_PROVIDER
@@ -191,6 +192,7 @@ BEGIN_PROVIDER [ double precision, mo_l_coef, (ao_num, mo_num) ]
         mo_l_coef(j,i) = mo_coef(j,i)
       enddo
     enddo
+    call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
   endif
 
 END_PROVIDER
diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index b8beea76..16414f39 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -1,53 +1,64 @@
 
- BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
- implicit none
- BEGIN_DOC
- ! aos_in_r_array(i,j)        = value of the ith ao on the jth grid point
- END_DOC
- integer :: i,j
- double precision :: aos_array(ao_num), r(3)
- !$OMP PARALLEL DO &
- !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (i,r,aos_array,j) & 
- !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
- do i = 1, n_points_final_grid
-  r(1) = final_grid_points(1,i)
-  r(2) = final_grid_points(2,i)
-  r(3) = final_grid_points(3,i)
-  call give_all_aos_at_r(r,aos_array)
-  do j = 1, ao_num
-   aos_in_r_array(j,i) = aos_array(j)
+! ---
+
+BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
+
+  BEGIN_DOC
+  ! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: tmp_array(ao_num), r(3)
+
+  !$OMP PARALLEL DO               &
+  !$OMP DEFAULT (NONE)            &
+  !$OMP PRIVATE (i,r,tmp_array,j) & 
+  !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
+  do i = 1, n_points_final_grid
+    r(1) = final_grid_points(1,i)
+    r(2) = final_grid_points(2,i)
+    r(3) = final_grid_points(3,i)
+    call give_all_aos_at_r(r, tmp_array)
+    do j = 1, ao_num
+      aos_in_r_array(j,i) = tmp_array(j)
+    enddo
   enddo
- enddo
- !$OMP END PARALLEL DO
+  !$OMP END PARALLEL DO
 
- END_PROVIDER
+END_PROVIDER
 
+! ---
 
- BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
- implicit none
- BEGIN_DOC
- ! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
- END_DOC
- integer :: i,j
- double precision :: aos_array(ao_num), r(3)
- do i = 1, n_points_final_grid
-  do j = 1, ao_num
-   aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) 
+BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
+
+  BEGIN_DOC
+  ! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
+  END_DOC
+
+  implicit none
+  integer          :: i, j
+  double precision :: aos_array(ao_num), r(3)
+
+  do i = 1, n_points_final_grid
+    do j = 1, ao_num
+      aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) 
+    enddo
   enddo
- enddo
 
- END_PROVIDER
+END_PROVIDER
 
+! ---
 
+BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
 
- BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
- implicit none
  BEGIN_DOC
  ! aos_grad_in_r_array(i,j,k)          = value of the kth component of the gradient of ith ao on the jth grid point
  !
  ! k = 1 : x, k= 2, y, k  3, z
  END_DOC
+
+ implicit none
  integer :: i,j,m
  double precision :: aos_array(ao_num), r(3)
  double precision :: aos_grad_array(3,ao_num)
diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f
index afa10305..158ee2fb 100644
--- a/src/non_h_ints_mu/total_tc_int.irp.f
+++ b/src/non_h_ints_mu/total_tc_int.irp.f
@@ -1,4 +1,7 @@
 
+! TODO
+! remove ao_two_e_coul and use map directly
+
 ! ---
 
 BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
@@ -116,6 +119,7 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, a
 
   call wall_time(wall1)
   print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
+
 END_PROVIDER 
 
 ! ---
diff --git a/src/tc_bi_ortho/print_tc_dump.irp.f b/src/tc_bi_ortho/print_tc_dump.irp.f
new file mode 100644
index 00000000..327e0f02
--- /dev/null
+++ b/src/tc_bi_ortho/print_tc_dump.irp.f
@@ -0,0 +1,44 @@
+program tc_bi_ortho
+
+  BEGIN_DOC
+  ! TODO
+  END_DOC
+  implicit none
+
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  call KMat_tilde_dump()
+end
+
+! ---
+
+subroutine KMat_tilde_dump()
+
+  implicit none
+  integer                       :: i, j, k, l
+
+  PROVIDE mo_bi_ortho_tc_two_e_chemist
+
+  print *, ' Kmat_tilde in chem notation'
+
+  open(33, file='Kmat_tilde.dat', action='write')  
+    do l = 1, mo_num
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do i = 1, mo_num
+            write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(i,j,k,l)
+            ! TCHint convention
+            !write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(j,i,l,k)
+          enddo
+        enddo
+      enddo
+    enddo
+  close(33)
+
+  return
+end subroutine KMat_tilde_dump
+
+! ---
diff --git a/src/tc_bi_ortho/tc_utils.irp.f b/src/tc_bi_ortho/tc_utils.irp.f
index 9023e2f0..53fe5884 100644
--- a/src/tc_bi_ortho/tc_utils.irp.f
+++ b/src/tc_bi_ortho/tc_utils.irp.f
@@ -5,27 +5,34 @@ subroutine write_tc_energy()
   integer          :: i, j, k
   double precision :: hmono, htwoe, hthree, htot
   double precision :: E_TC, O_TC
+  double precision :: E_1e, E_2e, E_3e
 
   do k = 1, n_states
 
     E_TC = 0.d0
+    E_1e = 0.d0
+    E_2e = 0.d0
+    E_3e = 0.d0
     do i = 1, N_det
       do j = 1, N_det
-        !htot = htilde_matrix_elmt_bi_ortho(i,j)
         call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
         E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
-        !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
+        E_1e = E_1e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hmono
+        E_2e = E_2e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htwoe
+        E_3e = E_3e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hthree
       enddo
     enddo
 
     O_TC = 0.d0
     do i = 1, N_det
-      !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
       O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
     enddo
 
     print *, ' state :', k
     print *, " E_TC = ", E_TC / O_TC
+    print *, " E_1e = ", E_1e / O_TC
+    print *, " E_2e = ", E_2e / O_TC
+    print *, " E_3e = ", E_3e / O_TC
     print *, " O_TC = ", O_TC
 
   enddo
diff --git a/src/tc_scf/fock_three_hermit.irp.f b/src/tc_scf/fock_three_hermit.irp.f
index fe8fbfd7..89e6f620 100644
--- a/src/tc_scf/fock_three_hermit.irp.f
+++ b/src/tc_scf/fock_three_hermit.irp.f
@@ -1,30 +1,36 @@
+
+! ---
+
 BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)] 
- implicit none
+
+  implicit none
   integer :: i,j
   double precision :: contrib
+
   fock_3_mat = 0.d0
-  if(.not.bi_ortho.and.three_body_h_tc)then
-   call give_fock_ia_three_e_total(1,1,contrib)
-!!  !$OMP PARALLEL                  &
-!!  !$OMP DEFAULT (NONE)            &
-!!  !$OMP PRIVATE (i,j,m,integral) & 
-!!  !$OMP SHARED (mo_num,three_body_3_index)
-!!  !$OMP DO SCHEDULE (guided) COLLAPSE(3)
-   do i = 1, mo_num
-    do j = 1, mo_num
-     call give_fock_ia_three_e_total(j,i,contrib)
-     fock_3_mat(j,i) = -contrib
-    enddo
-   enddo
-  else if(bi_ortho.and.three_body_h_tc)then
-!!  !$OMP END DO
-!!  !$OMP END PARALLEL
-!!  do i = 1, mo_num
-!!   do j = 1, i-1
-!!    mat_three(j,i) = mat_three(i,j)
-!!   enddo
-!!  enddo
- endif
+  if(.not.bi_ortho .and. three_body_h_tc) then
+
+    call give_fock_ia_three_e_total(1, 1, contrib)
+    !!  !$OMP PARALLEL                  &
+    !!  !$OMP DEFAULT (NONE)            &
+    !!  !$OMP PRIVATE (i,j,m,integral) & 
+    !!  !$OMP SHARED (mo_num,three_body_3_index)
+    !!  !$OMP DO SCHEDULE (guided) COLLAPSE(3)
+    do i = 1, mo_num
+      do j = 1, mo_num
+        call give_fock_ia_three_e_total(j,i,contrib)
+        fock_3_mat(j,i) = -contrib
+       enddo
+     enddo
+    !else if(bi_ortho.and.three_body_h_tc) then
+    !!  !$OMP END DO
+    !!  !$OMP END PARALLEL
+    !!  do i = 1, mo_num
+    !!   do j = 1, i-1
+    !!    mat_three(j,i) = mat_three(i,j)
+    !!   enddo
+    !!  enddo
+  endif
 
 END_PROVIDER 
 
diff --git a/src/tc_scf/molden_lr_mos.irp.f b/src/tc_scf/molden_lr_mos.irp.f
index 735349ba..e12fcd1c 100644
--- a/src/tc_scf/molden_lr_mos.irp.f
+++ b/src/tc_scf/molden_lr_mos.irp.f
@@ -1,6 +1,9 @@
-program molden
+! ---
+
+program molden_lr_mos
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! TODO : Put the documentation of the program here
   END_DOC
 
   implicit none
@@ -14,13 +17,21 @@ program molden
 !  my_n_pt_a_grid = 26 ! small grid for quick debug
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
-  call molden_lr
+  !call molden_lr
+  call molden_l()
+  call molden_r()
+
 end
+
+! ---
+
 subroutine molden_lr
-  implicit none
+
   BEGIN_DOC
   ! Produces a Molden file
   END_DOC
+
+  implicit none
   character*(128)                :: output
   integer                        :: i_unit_output,getUnitAndOpen
   integer                        :: i,j,k,l
@@ -37,7 +48,7 @@ subroutine molden_lr
 
   write(i_unit_output,'(A)') '[Atoms] Angs'
   do i = 1, nucl_num
-    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')                   &
+    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')             &
         trim(element_name(int(nucl_charge(i)))),                     &
         i,                                                           &
         int(nucl_charge(i)),                                         &
@@ -174,3 +185,314 @@ subroutine molden_lr
   close(i_unit_output)
 end
 
+! ---
+
+subroutine molden_l()
+
+  BEGIN_DOC
+  ! Produces a Molden file
+  END_DOC
+
+  implicit none
+  character*(128)             :: output
+  integer                     :: i_unit_output, getUnitAndOpen
+  integer                     :: i, j, k, l
+  double precision, parameter :: a0 = 0.529177249d0
+
+  PROVIDE ezfio_filename
+  PROVIDE mo_l_coef
+
+  output=trim(ezfio_filename)//'_left.mol'
+  print*,'output = ',trim(output)
+
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  write(i_unit_output,'(A)') '[Molden Format]'
+
+  write(i_unit_output,'(A)') '[Atoms] Angs'
+  do i = 1, nucl_num
+    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')             &
+        trim(element_name(int(nucl_charge(i)))),                     &
+        i,                                                           &
+        int(nucl_charge(i)),                                         &
+        nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
+  enddo
+
+  write(i_unit_output,'(A)') '[GTO]'
+
+  character*(1)                  :: character_shell
+  integer                        :: i_shell,i_prim,i_ao
+  integer                        :: iorder(ao_num)
+  integer                        :: nsort(ao_num)
+
+  i_shell = 0
+  i_prim = 0
+  do i=1,nucl_num
+    write(i_unit_output,*) i, 0
+    do j=1,nucl_num_shell_aos(i)
+      i_shell +=1
+      i_ao = nucl_list_shell_aos(i,j)
+      character_shell = trim(ao_l_char(i_ao))
+      write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
+      do k = 1, ao_prim_num(i_ao)
+        i_prim +=1
+        write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
+      enddo
+      l = i_ao
+      do while ( ao_l(l) == ao_l(i_ao) )
+        nsort(l) = i*10000 + j*100
+        l += 1
+        if (l > ao_num) exit
+      enddo
+    enddo
+    write(i_unit_output,*)''
+  enddo
+
+
+  do i=1,ao_num
+    iorder(i) = i
+    ! p
+    if      ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 3
+    ! d
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    ! f
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 10
+    ! g
+    else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 10
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 11
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 12
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 13
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 14
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 15
+    endif
+  enddo
+
+  call isort(nsort,iorder,ao_num)
+  write(i_unit_output,'(A)') '[MO]'
+  do i=1,mo_num
+    write (i_unit_output,*) 'Sym= 1'
+    write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
+    write (i_unit_output,*) 'Spin= Alpha'
+    write (i_unit_output,*) 'Occup=', mo_occ(i)
+    do j=1,ao_num
+      write(i_unit_output, '(I6,2X,E20.10)') j, mo_l_coef(iorder(j),i)
+    enddo
+  enddo
+  close(i_unit_output)
+
+end
+
+! ---
+
+subroutine molden_r()
+
+  BEGIN_DOC
+  ! Produces a Molden file
+  END_DOC
+
+  implicit none
+  character*(128)             :: output
+  integer                     :: i_unit_output, getUnitAndOpen
+  integer                     :: i, j, k, l
+  double precision, parameter :: a0 = 0.529177249d0
+
+  PROVIDE ezfio_filename
+
+  output=trim(ezfio_filename)//'_right.mol'
+  print*,'output = ',trim(output)
+
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  write(i_unit_output,'(A)') '[Molden Format]'
+
+  write(i_unit_output,'(A)') '[Atoms] Angs'
+  do i = 1, nucl_num
+    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')             &
+        trim(element_name(int(nucl_charge(i)))),                     &
+        i,                                                           &
+        int(nucl_charge(i)),                                         &
+        nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
+  enddo
+
+  write(i_unit_output,'(A)') '[GTO]'
+
+  character*(1)                  :: character_shell
+  integer                        :: i_shell,i_prim,i_ao
+  integer                        :: iorder(ao_num)
+  integer                        :: nsort(ao_num)
+
+  i_shell = 0
+  i_prim = 0
+  do i=1,nucl_num
+    write(i_unit_output,*) i, 0
+    do j=1,nucl_num_shell_aos(i)
+      i_shell +=1
+      i_ao = nucl_list_shell_aos(i,j)
+      character_shell = trim(ao_l_char(i_ao))
+      write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
+      do k = 1, ao_prim_num(i_ao)
+        i_prim +=1
+        write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
+      enddo
+      l = i_ao
+      do while ( ao_l(l) == ao_l(i_ao) )
+        nsort(l) = i*10000 + j*100
+        l += 1
+        if (l > ao_num) exit
+      enddo
+    enddo
+    write(i_unit_output,*)''
+  enddo
+
+
+  do i=1,ao_num
+    iorder(i) = i
+    ! p
+    if      ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 3
+    ! d
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    ! f
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 10
+    ! g
+    else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 10
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 11
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 12
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 13
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 14
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 15
+    endif
+  enddo
+
+  call isort(nsort, iorder, ao_num)
+  write(i_unit_output,'(A)') '[MO]'
+  do i=1,mo_num
+    write (i_unit_output,*) 'Sym= 1'
+    write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
+    write (i_unit_output,*) 'Spin= Alpha'
+    write (i_unit_output,*) 'Occup=', mo_occ(i)
+    do j=1,ao_num
+      write(i_unit_output, '(I6,2X,E20.10)') j, mo_r_coef(iorder(j),i)
+    enddo
+  enddo
+  close(i_unit_output)
+
+end
+