diff --git a/TODO b/TODO
index 734ed083..046510ed 100644
--- a/TODO
+++ b/TODO
@@ -75,4 +75,6 @@ Davidson Diagonalization
 ------------------------
 
  Not enough memory: aborting in davidson_diag_hjj_sjj
->>>>>>> 94bacff2d093aa9b32c653ab59bcdb79d13f3264
+
+
+qp man does not find the programs in external plugins
diff --git a/ocaml/.gitignore b/ocaml/.gitignore
index 74e83c02..228bb7f6 100644
--- a/ocaml/.gitignore
+++ b/ocaml/.gitignore
@@ -9,12 +9,13 @@ Input_ao_two_e_erf_ints.ml
 Input_ao_two_e_ints.ml
 Input_auto_generated.ml
 Input_becke_numerical_grid.ml
-Input_champ.ml
 Input_davidson.ml
 Input_density_for_dft.ml
 Input_determinants.ml
 Input_dft_keywords.ml
 Input_dressing.ml
+Input_firth_order_der.ml
+Input_ijkl_ints_in_r3.ml
 Input_mo_one_e_ints.ml
 Input_mo_two_e_erf_ints.ml
 Input_mo_two_e_ints.ml
@@ -22,7 +23,6 @@ Input_nuclei.ml
 Input_perturbation.ml
 Input_pseudo.ml
 Input_scf_utils.ml
-Input_variance.ml
 qp_create_ezfio
 qp_create_ezfio.native
 qp_edit
diff --git a/src/dft_utils_in_r/dm_in_r.irp.f b/src/dft_utils_in_r/dm_in_r.irp.f
index 8d277453..774061cf 100644
--- a/src/dft_utils_in_r/dm_in_r.irp.f
+++ b/src/dft_utils_in_r/dm_in_r.irp.f
@@ -93,7 +93,6 @@ end
    grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
    grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
    grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
-   grad_dm_a *= 2.d0
    ! aos_grad_array_bis = \rho_ao * aos_grad_array
 
    ! beta density
@@ -104,9 +103,10 @@ end
    grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
    grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
    grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
-   grad_dm_b *= 2.d0
    ! aos_grad_array_bis = \rho_ao * aos_grad_array
  enddo
+   grad_dm_a *= 2.d0
+   grad_dm_b *= 2.d0
  end
 
  BEGIN_PROVIDER [double precision, one_e_dm_alpha_in_r, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
diff --git a/src/tools/NEED b/src/tools/NEED
index 98afd7a1..c07c9109 100644
--- a/src/tools/NEED
+++ b/src/tools/NEED
@@ -1,3 +1,4 @@
 fci
 mo_two_e_erf_ints
 aux_quantities
+hartree_fock
diff --git a/src/tools/sort_by_fock_energies.irp.f b/src/tools/sort_by_fock_energies.irp.f
new file mode 100644
index 00000000..89b30780
--- /dev/null
+++ b/src/tools/sort_by_fock_energies.irp.f
@@ -0,0 +1,32 @@
+program sort_by_fock_energies
+ BEGIN_DOC
+! programs that save the current mos ordered by Diagonal element of the Fock operator. 
+!
+! Warning : the Fock operator, and therefore its matrix elements, depends on the occupancy. 
+ END_DOC
+ implicit none
+ integer :: i,j,k
+ integer, allocatable :: iorder(:)
+ double precision, allocatable :: fock_energies_tmp(:), new_mo_coef(:,:)
+ allocate(iorder(mo_num), fock_energies_tmp(mo_num),new_mo_coef(ao_num,mo_num))
+ do i = 1, mo_num
+   fock_energies_tmp(i) = Fock_matrix_diag_mo(i)
+   print*,'fock_energies_tmp(i) = ',fock_energies_tmp(i)
+   iorder(i) = i
+ enddo
+ print*,''
+ print*,'Sorting by Fock energies'
+ print*,''
+ call dsort(fock_energies_tmp,iorder,mo_num)
+ do i = 1, mo_num
+  k = iorder(i) 
+  print*,'fock_energies_new(i) = ',fock_energies_tmp(i)
+  do j = 1, ao_num
+   new_mo_coef(j,i) = mo_coef(j,k)
+  enddo
+ enddo
+ mo_coef = new_mo_coef
+ touch mo_coef
+ call save_mos
+
+end
diff --git a/src/tools/swap_mos.irp.f b/src/tools/swap_mos.irp.f
new file mode 100644
index 00000000..80a29965
--- /dev/null
+++ b/src/tools/swap_mos.irp.f
@@ -0,0 +1,14 @@
+program swap_mos
+  implicit none
+  integer :: i,j, i1, i2
+  double precision :: x
+  print *,  'MOs to swap?'
+  read(*,*) i1, i2
+  do i=1,ao_num
+    x = mo_coef(i,i1)
+    mo_coef(i,i1) = mo_coef(i,i2)
+    mo_coef(i,i2) = x
+  enddo
+  call save_mos
+
+end