diff --git a/docs/source/modules/cipsi.rst b/docs/source/modules/cipsi.rst
index cb5fb19d..5415a4e2 100644
--- a/docs/source/modules/cipsi.rst
+++ b/docs/source/modules/cipsi.rst
@@ -1213,6 +1213,7 @@ Subroutines / functions
        :columns: 3
 
        * :c:data:`generators_bitmask`
+       * :c:data:`selection_weight`
        * :c:data:`pt2_stoch_istate`
        * :c:data:`psi_selectors`
        * :c:data:`psi_det`
@@ -1235,6 +1236,7 @@ Subroutines / functions
        * :c:data:`n_det_generators`
        * :c:data:`psi_det_generators`
        * :c:data:`n_int`
+       * :c:data:`pseudo_sym`
        * :c:data:`psi_det_generators`
 
     Called by:
@@ -1590,6 +1592,7 @@ Subroutines / functions
        * :c:data:`h_apply_buffer_allocated`
        * :c:data:`n_det`
        * :c:data:`s2_eig`
+       * :c:data:`variance_max`
        * :c:data:`do_pt2`
        * :c:data:`psi_energy`
        * :c:data:`pt2_relative_error`
@@ -1810,16 +1813,18 @@ Subroutines / functions
        * :c:data:`psi_bilinear_matrix_values`
        * :c:data:`psi_bilinear_matrix_transp_values`
        * :c:data:`mo_num`
+       * :c:data:`psi_det_sorted`
        * :c:data:`psi_bilinear_matrix_transp_rows_loc`
        * :c:data:`n_states`
        * :c:data:`pt2_f`
-       * :c:data:`psi_bilinear_matrix_values`
        * :c:data:`elec_alpha_num`
+       * :c:data:`psi_bilinear_matrix_values`
+       * :c:data:`weight_selection`
        * :c:data:`psi_bilinear_matrix_transp_values`
        * :c:data:`n_int`
        * :c:data:`psi_det_alpha_unique`
        * :c:data:`psi_det_sorted`
-       * :c:data:`psi_det_sorted`
+       * :c:data:`pseudo_sym`
        * :c:data:`psi_bilinear_matrix_columns_loc`
        * :c:data:`psi_selectors_coef_transp`
 
@@ -1893,6 +1898,7 @@ Subroutines / functions
        * :c:data:`pt2_e0_denominator`
        * :c:data:`pt2_stoch_istate`
        * :c:data:`read_wf`
+       * :c:data:`selection_weight`
        * :c:data:`state_average_weight`
        * :c:data:`threshold_generators`
 
@@ -1913,6 +1919,7 @@ Subroutines / functions
     .. hlist::
        :columns: 3
 
+       * :c:data:`selection_weight`
        * :c:data:`pt2_stoch_istate`
        * :c:data:`psi_det`
        * :c:data:`zmq_state`
@@ -1969,6 +1976,7 @@ Subroutines / functions
 
        * :c:data:`pt2_e0_denominator`
        * :c:data:`pt2_stoch_istate`
+       * :c:data:`selection_weight`
        * :c:data:`state_average_weight`
        * :c:data:`threshold_generators`
 
@@ -2413,6 +2421,7 @@ Subroutines / functions
        :columns: 3
 
        * :c:func:`zmq_pt2`
+       * :c:func:`zmq_selection`
 
     Touches:
 
@@ -2440,6 +2449,7 @@ Subroutines / functions
        :columns: 3
 
        * :c:data:`psi_det_sorted`
+       * :c:data:`selection_weight`
        * :c:data:`pt2_stoch_istate`
        * :c:data:`psi_selectors`
        * :c:data:`psi_bilinear_matrix_values`
@@ -2450,6 +2460,7 @@ Subroutines / functions
        * :c:data:`state_average_weight`
        * :c:data:`n_det`
        * :c:data:`s2_eig`
+       * :c:data:`psi_det_sorted`
        * :c:data:`pt2_j`
        * :c:data:`mo_two_e_integrals_in_map`
        * :c:data:`psi_bilinear_matrix_transp_values`
@@ -2467,7 +2478,7 @@ Subroutines / functions
        * :c:data:`threshold_generators`
        * :c:data:`psi_det_beta_unique`
        * :c:data:`det_to_occ_pattern`
-       * :c:data:`psi_det_sorted`
+       * :c:data:`global_selection_buffer`
        * :c:data:`psi_bilinear_matrix_transp_rows_loc`
        * :c:data:`n_states`
        * :c:data:`pt2_f`
@@ -2476,7 +2487,7 @@ Subroutines / functions
        * :c:data:`n_int`
        * :c:data:`psi_det_hii`
        * :c:data:`pt2_j`
-       * :c:data:`global_selection_buffer`
+       * :c:data:`pseudo_sym`
        * :c:data:`pt2_w`
        * :c:data:`pt2_u`
 
@@ -2545,12 +2556,12 @@ Subroutines / functions
        :columns: 3
 
        * :c:data:`psi_det_sorted`
+       * :c:data:`selection_weight`
        * :c:data:`psi_selectors`
        * :c:data:`psi_bilinear_matrix_values`
        * :c:data:`n_det`
        * :c:data:`psi_bilinear_matrix_columns_loc`
        * :c:data:`n_det_selectors`
-       * :c:data:`psi_bilinear_matrix_transp_values`
        * :c:data:`psi_det_alpha_unique`
        * :c:data:`psi_bilinear_matrix_transp_values`
        * :c:data:`state_average_weight`
@@ -2566,8 +2577,9 @@ Subroutines / functions
        * :c:data:`n_states`
        * :c:data:`pt2_f`
        * :c:data:`n_det_generators`
-       * :c:data:`pt2_match_weight`
+       * :c:data:`psi_bilinear_matrix_transp_values`
        * :c:data:`n_int`
+       * :c:data:`pseudo_sym`
 
     Called by:
 
@@ -2592,6 +2604,7 @@ Subroutines / functions
        * :c:func:`save_wavefunction`
        * :c:func:`selection_collector`
        * :c:func:`selection_slave_inproc`
+       * :c:func:`update_pt2_and_variance_weights`
        * :c:func:`write_double`
 
     Touches:
@@ -2607,4 +2620,5 @@ Subroutines / functions
        * :c:data:`psi_det_size`
        * :c:data:`psi_det_sorted_bit`
        * :c:data:`pt2_match_weight`
+       * :c:data:`variance_match_weight`
 
diff --git a/docs/source/modules/iterations.rst b/docs/source/modules/iterations.rst
index ab1d5065..81adc6a9 100644
--- a/docs/source/modules/iterations.rst
+++ b/docs/source/modules/iterations.rst
@@ -160,6 +160,13 @@ Subroutines / functions
        * :c:func:`run_cipsi`
        * :c:func:`run_stochastic_cipsi`
 
+    Calls:
+
+    .. hlist::
+       :columns: 3
+
+       * :c:func:`print_energy_components`
+
  
 .. c:function:: save_iterations:
 
diff --git a/docs/source/programs/fci.rst b/docs/source/programs/fci.rst
index e0441d91..5979c28d 100644
--- a/docs/source/programs/fci.rst
+++ b/docs/source/programs/fci.rst
@@ -89,6 +89,7 @@ fci
     * :c:data:`pt2_match_weight` 
     * :c:data:`pt2_stoch_istate` 
     * :c:data:`read_wf` 
+    * :c:data:`selection_weight` 
     * :c:data:`state_average_weight` 
     * :c:data:`threshold_generators` 
     * :c:data:`variance_match_weight` 
diff --git a/docs/source/programs/pt2.rst b/docs/source/programs/pt2.rst
index 45726526..4878bb4d 100644
--- a/docs/source/programs/pt2.rst
+++ b/docs/source/programs/pt2.rst
@@ -55,5 +55,6 @@ pt2
     * :c:data:`pt2_e0_denominator` 
     * :c:data:`pt2_stoch_istate` 
     * :c:data:`read_wf` 
+    * :c:data:`selection_weight` 
     * :c:data:`state_average_weight` 
     * :c:data:`threshold_generators` 
diff --git a/docs/source/programs/rotate_mos.rst b/docs/source/programs/rotate_mos.rst
new file mode 100644
index 00000000..e5b4353c
--- /dev/null
+++ b/docs/source/programs/rotate_mos.rst
@@ -0,0 +1,37 @@
+.. _rotate_mos: 
+ 
+.. program:: rotate_mos 
+ 
+========== 
+rotate_mos 
+========== 
+ 
+ 
+ 
+ 
+ Rotates molecular orbitals i and j by combining them as 
+ $1/\sqrt{2} ( \phi_i + \phi_j )$ and 
+ $1/\sqrt{2} ( \phi_i - \phi_j )$. 
+ 
+ Needs: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:data:`ao_num` 
+    * :c:data:`mo_num` 
+    * :c:data:`mo_coef` 
+ 
+ Calls: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:func:`save_mos` 
+ 
+ Touches: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:data:`mo_coef` 
diff --git a/docs/source/programs/sort_by_fock_energies.rst b/docs/source/programs/sort_by_fock_energies.rst
new file mode 100644
index 00000000..5d0e0d9f
--- /dev/null
+++ b/docs/source/programs/sort_by_fock_energies.rst
@@ -0,0 +1,39 @@
+.. _sort_by_fock_energies: 
+ 
+.. program:: sort_by_fock_energies 
+ 
+===================== 
+sort_by_fock_energies 
+===================== 
+ 
+ 
+ 
+ 
+ Program that saves the current |MOs| ordered by diagonal element of the Fock operator. 
+  
+ Warning : the Fock operator, and therefore its matrix elements, depends on the occupancy. 
+ 
+ Needs: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:data:`fock_matrix_mo` 
+    * :c:data:`ao_num` 
+    * :c:data:`mo_num` 
+    * :c:data:`mo_coef` 
+ 
+ Calls: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:func:`dsort` 
+    * :c:func:`save_mos` 
+ 
+ Touches: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:data:`mo_coef` 
diff --git a/docs/source/programs/swap_mos.rst b/docs/source/programs/swap_mos.rst
new file mode 100644
index 00000000..c2d3440e
--- /dev/null
+++ b/docs/source/programs/swap_mos.rst
@@ -0,0 +1,27 @@
+.. _swap_mos: 
+ 
+.. program:: swap_mos 
+ 
+======== 
+swap_mos 
+======== 
+ 
+ 
+ 
+ 
+ Swaps the indices of two molecular orbitals 
+ 
+ Needs: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:data:`ao_num` 
+    * :c:data:`mo_coef` 
+ 
+ Calls: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:func:`save_mos` 
diff --git a/docs/source/programs/test.rst b/docs/source/programs/test.rst
new file mode 100644
index 00000000..6595f24d
--- /dev/null
+++ b/docs/source/programs/test.rst
@@ -0,0 +1,19 @@
+.. _test: 
+ 
+.. program:: test 
+ 
+==== 
+test 
+==== 
+ 
+ 
+ 
+ 
+ 
+ Calls: 
+ 
+ .. hlist:: 
+    :columns: 3 
+ 
+    * :c:func:`two_e_integrals_index` 
+    * :c:func:`two_e_integrals_index_reverse` 
diff --git a/man/fci.1 b/man/fci.1
index 9f565169..e4b4dd50 100644
--- a/man/fci.1
+++ b/man/fci.1
@@ -125,10 +125,10 @@ Touches:
 \fBpsi_occ_pattern\fP
 .IP \(bu 2
 \fBc0_weight\fP
-.UNINDENT
-.INDENT 2.0
 .IP \(bu 2
 \fBdistributed_davidson\fP
+.UNINDENT
+.INDENT 2.0
 .IP \(bu 2
 \fBpsi_coef\fP
 .IP \(bu 2
@@ -143,10 +143,10 @@ Touches:
 \fBpsi_energy\fP
 .IP \(bu 2
 \fBpsi_occ_pattern\fP
-.UNINDENT
-.INDENT 2.0
 .IP \(bu 2
 \fBpsi_energy\fP
+.UNINDENT
+.INDENT 2.0
 .IP \(bu 2
 \fBpt2_e0_denominator\fP
 .IP \(bu 2
@@ -156,6 +156,8 @@ Touches:
 .IP \(bu 2
 \fBread_wf\fP
 .IP \(bu 2
+\fBselection_weight\fP
+.IP \(bu 2
 \fBstate_average_weight\fP
 .IP \(bu 2
 \fBthreshold_generators\fP
diff --git a/man/pt2.1 b/man/pt2.1
index 2c70c66e..5902f81a 100644
--- a/man/pt2.1
+++ b/man/pt2.1
@@ -104,6 +104,8 @@ Touches:
 .IP \(bu 2
 \fBread_wf\fP
 .IP \(bu 2
+\fBselection_weight\fP
+.IP \(bu 2
 \fBstate_average_weight\fP
 .IP \(bu 2
 \fBthreshold_generators\fP
diff --git a/man/rotate_mos.1 b/man/rotate_mos.1
new file mode 100644
index 00000000..54acbc0d
--- /dev/null
+++ b/man/rotate_mos.1
@@ -0,0 +1,85 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "ROTATE_MOS" "1" "Jun 15, 2019" "2.0" "Quantum Package"
+.SH NAME
+rotate_mos \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.INDENT 0.0
+.INDENT 3.5
+Rotates molecular orbitals i and j by combining them as
+$1/sqrt{2} ( phi_i + phi_j )$ and
+$1/sqrt{2} ( phi_i \- phi_j )$.
+.sp
+Needs:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBao_num\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_num\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_coef\fP
+.UNINDENT
+.UNINDENT
+.sp
+Calls:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBsave_mos()\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.sp
+Touches:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_coef\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2019, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/sort_by_fock_energies.1 b/man/sort_by_fock_energies.1
new file mode 100644
index 00000000..ac6dd64c
--- /dev/null
+++ b/man/sort_by_fock_energies.1
@@ -0,0 +1,89 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "SORT_BY_FOCK_ENERGIES" "1" "Jun 15, 2019" "2.0" "Quantum Package"
+.SH NAME
+sort_by_fock_energies \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.INDENT 0.0
+.INDENT 3.5
+Program that saves the current MOs ordered by diagonal element of the Fock operator.
+.sp
+Warning : the Fock operator, and therefore its matrix elements, depends on the occupancy.
+.sp
+Needs:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBfock_matrix_mo\fP
+.IP \(bu 2
+\fBao_num\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_num\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_coef\fP
+.UNINDENT
+.UNINDENT
+.sp
+Calls:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBdsort()\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBsave_mos()\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.sp
+Touches:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_coef\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2019, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/man/swap_mos.1 b/man/swap_mos.1
new file mode 100644
index 00000000..73bfc4b0
--- /dev/null
+++ b/man/swap_mos.1
@@ -0,0 +1,69 @@
+.\" Man page generated from reStructuredText.
+.
+.TH "SWAP_MOS" "1" "Jun 15, 2019" "2.0" "Quantum Package"
+.SH NAME
+swap_mos \-  | Quantum Package >
+.
+.nr rst2man-indent-level 0
+.
+.de1 rstReportMargin
+\\$1 \\n[an-margin]
+level \\n[rst2man-indent-level]
+level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
+-
+\\n[rst2man-indent0]
+\\n[rst2man-indent1]
+\\n[rst2man-indent2]
+..
+.de1 INDENT
+.\" .rstReportMargin pre:
+. RS \\$1
+. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin]
+. nr rst2man-indent-level +1
+.\" .rstReportMargin post:
+..
+.de UNINDENT
+. RE
+.\" indent \\n[an-margin]
+.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.nr rst2man-indent-level -1
+.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]]
+.in \\n[rst2man-indent\\n[rst2man-indent-level]]u
+..
+.INDENT 0.0
+.INDENT 3.5
+Swaps the indices of two molecular orbitals
+.sp
+Needs:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBao_num\fP
+.UNINDENT
+.INDENT 2.0
+.IP \(bu 2
+\fBmo_coef\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.sp
+Calls:
+.INDENT 0.0
+.INDENT 2.0
+.IP \(bu 2
+\fBsave_mos()\fP
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.INDENT 2.0
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.UNINDENT
+.SH AUTHOR
+A. Scemama, E. Giner
+.SH COPYRIGHT
+2019, A. Scemama, E. Giner
+.\" Generated by docutils manpage writer.
+.
diff --git a/src/tools/rotate_mos.irp.f b/src/tools/rotate_mos.irp.f
index 3e79c4f7..30185474 100644
--- a/src/tools/rotate_mos.irp.f
+++ b/src/tools/rotate_mos.irp.f
@@ -1,18 +1,27 @@
 program rotate_mos
- implicit none
- integer :: iorb,jorb
- read(5,*)iorb,jorb
- double precision, allocatable :: mo_coef_tmp(:,:)
- allocate(mo_coef_tmp(ao_num,mo_num))
- mo_coef_tmp = mo_coef
- integer :: i,j
- double precision :: dsqrt2_inv
- dsqrt2_inv = 1.d0/dsqrt(2.d0)
- do i = 1, ao_num
-  mo_coef(i,iorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) + mo_coef_tmp(i,jorb) )
-  mo_coef(i,jorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) - mo_coef_tmp(i,jorb) )
- enddo
- touch mo_coef 
- call save_mos
+  implicit none
+  BEGIN_DOC
+  ! Rotates molecular orbitals i and j by combining them as
+  ! $1/\sqrt{2} ( \phi_i + \phi_j )$ and
+  ! $1/\sqrt{2} ( \phi_i - \phi_j )$.
+  END_DOC
+  integer                        :: iorb,jorb
+  integer                        :: i,j
+  double precision               :: dsqrt2_inv
+  double precision, allocatable  :: mo_coef_tmp(:,:)
 
+  read(5,*)iorb,jorb
+
+  allocate(mo_coef_tmp(ao_num,mo_num))
+  mo_coef_tmp = mo_coef
+
+  dsqrt2_inv = 1.d0/dsqrt(2.d0)
+  do i = 1, ao_num
+    mo_coef(i,iorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) + mo_coef_tmp(i,jorb) )
+    mo_coef(i,jorb) = dsqrt2_inv * ( mo_coef_tmp(i,iorb) - mo_coef_tmp(i,jorb) )
+  enddo
+
+  touch mo_coef
+  call save_mos
+  
 end
diff --git a/src/tools/sort_by_fock_energies.irp.f b/src/tools/sort_by_fock_energies.irp.f
index 89b30780..4f903c7b 100644
--- a/src/tools/sort_by_fock_energies.irp.f
+++ b/src/tools/sort_by_fock_energies.irp.f
@@ -1,32 +1,38 @@
 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
+  BEGIN_DOC
+  ! Program that saves 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
index 80a29965..bba9cb34 100644
--- a/src/tools/swap_mos.irp.f
+++ b/src/tools/swap_mos.irp.f
@@ -1,7 +1,10 @@
 program swap_mos
   implicit none
-  integer :: i,j, i1, i2
-  double precision :: x
+  BEGIN_DOC
+  ! Swaps the indices of two molecular orbitals
+  END_DOC
+  integer                        :: i,j, i1, i2
+  double precision               :: x
   print *,  'MOs to swap?'
   read(*,*) i1, i2
   do i=1,ao_num
@@ -10,5 +13,5 @@ program swap_mos
     mo_coef(i,i2) = x
   enddo
   call save_mos
-
+  
 end