diff --git a/MOBasis/Localisation.ml b/MOBasis/Localisation.ml
index 153ed9b..9ddc004 100644
--- a/MOBasis/Localisation.ml
+++ b/MOBasis/Localisation.ml
@@ -1241,24 +1241,24 @@ let localize mo_basis =
       else
       (let multipoles = AOBasis.multipole ao_basis in
       let n_mo = Mat.dim2 mat in
-      let phi_x_phi =
+      let x =
           Multipole.matrix_x multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
       in
-      let phi_y_phi =
+      let y =
           Multipole.matrix_y multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
       in
-      let phi_z_phi =
+      let z =
           Multipole.matrix_z multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
       in
-      let phi_r2_phi =
+      let r2 =
           Multipole.matrix_r2 multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
       in
       Vec.init n_mo (fun i ->
-      phi_r2_phi.{i,i} -. phi_x_phi.{i,i}**2. -. phi_y_phi.{i,i}**2. -. phi_z_phi.{i,i}**2.))
+      r2.{i,i} -. x.{i,i}*.x.{i,i} -. y.{i,i}*.y.{i,i} -. z.{i,i}*.z.{i,i}))
   in
 
   (* Sorting function *)
@@ -1300,15 +1300,15 @@ let localize mo_basis =
 
   (* Function to compute the variance *)
   let f_variance mat =
-          if (Vec.dim (distrib mat)) = 0
+          if (Vec.dim (distrib mat)) <= 1
           then 0.
           else  Vec.sum(Vec.init (Vec.dim (distrib mat)) (fun i ->
-     ((distrib mat).{i}-. (f_average mat))**2.)) /. float_of_int(Vec.dim (distrib mat))
+     ((distrib mat).{i}-. (f_average mat))**2.)) /. float_of_int(Vec.dim (distrib mat) - 1)
   in
 
   (* Function to compute the standard deviation *)
   let f_stand_dev mat =
-          if (Vec.dim (distrib mat)) = 0
+          if (Vec.dim (distrib mat)) <= 1
           then 0.
           else  sqrt(abs_float(f_variance mat))
   in