diff --git a/configure b/configure
index c32420a0..9b59b209 100755
--- a/configure
+++ b/configure
@@ -480,24 +480,24 @@ def create_ninja_and_rc(l_installed):
         'export QP_PYTHON={0}'.format(":".join(l_python)), "",
         'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
         'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
-        'qp_append_export () {',
-        '   #Append path $2:${!1}. Add the semicolon only if ${!1} is defined',
-        '   echo ${2}${!1:+:${!1}}',
+        'function qp_append_export () {',
+        '  #Append path $2:${!1}. Add the semicolon only if ${!1} is defined',
+        '  eval "value_1=\"\${$1}\""',
+        '  echo ${2}${value_1:+:${value_1}}',
         '}',
         'export PYTHONPATH=$(qp_append_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")',
         'export PATH=$(qp_append_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)',
         'export LD_LIBRARY_PATH=$(qp_append_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
         'export LIBRARY_PATH=$(qp_append_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)',
-        'export C_INCLUDE_PATH=$(qp_append_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)'
+        'export C_INCLUDE_PATH=$(qp_append_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)',
         '',
-        'source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh',
+        'if [[ $SHELL == "bash" ]] ; then',
+        '  source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh',
+        'fi',
         '',
         '# Choose the correct network interface',
         '# export QP_NIC=ib0',
         '# export QP_NIC=eth0',
-        '',
-        '# Choose how to start MPI processes',
-        '# export QP_MPIRUN="mpirun"',
         ''
     ]
 
diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml
index ab265202..a46be191 100644
--- a/ocaml/Gto.ml
+++ b/ocaml/Gto.ml
@@ -36,7 +36,10 @@ let of_prim_coef_list pc =
 
 let read_one in_channel =
   (* Fetch number of lines to read on first line *)
-  let buffer = input_line in_channel in
+  let buffer = 
+    try input_line in_channel with
+    | End_of_file -> raise End_Of_Basis
+  in
   if ( (String_ext.strip buffer) = "" ) then
      raise End_Of_Basis;
   let sym_str = String.sub buffer 0 2 in
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index 48887ca0..90174e18 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@@ -11,6 +11,7 @@ module Determinants_by_hand : sig
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
+      state_average_weight   : Positive_float.t array;
     } [@@deriving sexp]
   val read  : unit -> t
   val read_maybe  : unit -> t option
@@ -30,6 +31,7 @@ end = struct
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
+      state_average_weight   : Positive_float.t array;
     } [@@deriving sexp]
   ;;
 
@@ -82,7 +84,6 @@ end = struct
     |> Ezfio.set_determinants_n_det
   ;;
 
-
   let read_n_states () =
     if not (Ezfio.has_determinants_n_states ()) then
       Ezfio.set_determinants_n_states 1
@@ -96,6 +97,36 @@ end = struct
     |> Ezfio.set_determinants_n_states
   ;;
 
+  let write_state_average_weight data =
+      let n_states =
+        read_n_states ()
+        |> States_number.to_int
+      in
+      let data =
+        Array.map ~f:Positive_float.to_float data
+        |> Array.to_list
+      in
+      Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_states |] ~data
+      |> Ezfio.set_determinants_state_average_weight
+  ;;
+
+  let read_state_average_weight () =
+    if not (Ezfio.has_determinants_state_average_weight ()) then
+      begin
+        let n_states =
+          read_n_states ()
+          |> States_number.to_int
+        in
+        let data = 
+          Array.init n_states (fun _ -> 1./.(float_of_int n_states))
+          |> Array.map ~f:Positive_float.of_float
+        in
+        write_state_average_weight data;
+      end;
+    Ezfio.get_determinants_state_average_weight ()
+    |> Ezfio.flattened_ezfio
+    |> Array.map ~f:Positive_float.of_float
+  ;;
 
   let read_expected_s2 () =
     if not (Ezfio.has_determinants_expected_s2 ()) then
@@ -205,6 +236,7 @@ end = struct
           psi_coef               = read_psi_coef ()             ;
           psi_det                = read_psi_det ()              ;
           n_states               = read_n_states ()             ;
+          state_average_weight   = read_state_average_weight () ;
         }
     else
       failwith "No molecular orbitals, so no determinants"
@@ -228,6 +260,7 @@ end = struct
               psi_coef             ;
               psi_det              ;
               n_states             ;
+              state_average_weight ;
             } =
      write_n_int n_int ;
      write_bit_kind bit_kind;
@@ -236,6 +269,7 @@ end = struct
      write_expected_s2 expected_s2;
      write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
      write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
+     write_state_average_weight state_average_weight;
   ;;
 
 
@@ -288,12 +322,17 @@ Number of determinants ::
 
   n_det = %s
 
+State average weights ::
+
+  state_average_weight = (%s)
+
 Determinants ::
 
 %s
 "
      (b.expected_s2   |> Positive_float.to_string)
      (b.n_det         |> Det_number.to_string)
+     (b.state_average_weight |> Array.to_list |> List.map ~f:Positive_float.to_string |> String.concat ~sep:"\t")
      det_text
      |> Rst_string.of_string
   ;;
@@ -307,6 +346,7 @@ bit_kind               = %s
 n_det                  = %s
 n_states               = %s
 expected_s2            = %s
+state_average_weight   = %s
 psi_coef               = %s
 psi_det                = %s
 "
@@ -315,6 +355,7 @@ psi_det                = %s
      (b.n_det         |> Det_number.to_string)
      (b.n_states      |> States_number.to_string)
      (b.expected_s2   |> Positive_float.to_string)
+     (b.state_average_weight |> Array.to_list |> List.map ~f:Positive_float.to_string |> String.concat ~sep:",")
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")
      (b.psi_det   |> Array.to_list |> List.map ~f:(Determinant.to_string
diff --git a/ocaml/Zmatrix.ml b/ocaml/Zmatrix.ml
index 0aae3441..b963479f 100644
--- a/ocaml/Zmatrix.ml
+++ b/ocaml/Zmatrix.ml
@@ -196,17 +196,6 @@ let rotation_matrix axis angle =
       (2. *. (b *. d -. a *. c),
        2. *. (c *. d +. a *. b),
        a *. a +. d *. d -. b *. b -. c *. c)]
-(*
-     [(a *. a +. b *. b -. c *. c -. d *. d,
-       2. *. (b *. c +. a *. d),
-       2. *. (b *. d -. a *. c));
-      (2. *. (b *. c -. a *. d),
-       a *. a +. c *. c -.b *. b -. d *. d,
-       2. *. (c *. d +. a *. b));
-      (2. *. (b *. d +. a *. c),
-       2. *. (c *. d -. a *. b),
-       a *. a +. d *. d -. b *. b -. c *. c)]
-*)
       
 
 
diff --git a/plugins/CIS/super_ci.irp.f b/plugins/CIS/super_ci.irp.f
index 630b9599..979a48a4 100644
--- a/plugins/CIS/super_ci.irp.f
+++ b/plugins/CIS/super_ci.irp.f
@@ -12,12 +12,12 @@ subroutine super_CI
   integer                        :: k
   character                      :: save_char
 
-  call write_time(output_hartree_fock)
-  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )')           &
+  call write_time(6)
+  write(6,'(A4,X,A16, X, A16, X, A16 )')           &
       '====','================','================','================'
-  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )')           &
+  write(6,'(A4,X,A16, X, A16, X, A16 )')           &
       '  N ', 'Energy  ', 'Energy diff  ', 'Save  '
-  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )')           &
+  write(6,'(A4,X,A16, X, A16, X, A16 )')           &
       '====','================','================','================'
   
   E = HF_energy + 1.d0
@@ -39,7 +39,7 @@ subroutine super_CI
       save_char = ' '
     endif
     E_min = min(E,E_min)
-    write(output_hartree_fock,'(I4,X,F16.10, X, F16.10, X, A8 )') &
+    write(6,'(I4,X,F16.10, X, F16.10, X, A8 )') &
         k, E, delta_E, save_char
     if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then
       exit
@@ -55,8 +55,8 @@ subroutine super_CI
     TOUCH mo_coef
   enddo
   
-  write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )')           &
+  write(6,'(A4,X,A16, X, A16, X, A16 )')           &
       '====','================','================','================'
-  call write_time(output_hartree_fock)
+  call write_time(6)
 end
     
diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
index 5848cec0..9d1c50d4 100644
--- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
+++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
@@ -25,8 +25,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
   
   double precision               :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
   double precision, external     :: omp_get_wtime
+  double precision               :: state_average_weight_save(N_states), w(N_states)
   double precision               :: time
-  double precision               :: w(N_states)
   integer(ZMQ_PTR), external     :: new_zmq_to_qp_run_socket
   
   if (N_det < max(10,N_states)) then
@@ -35,18 +35,19 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
     error(:) = 0.d0
   else
     
+    state_average_weight_save(:) = state_average_weight(:)
     do pt2_stoch_istate=1,N_states
       SOFT_TOUCH pt2_stoch_istate
-      w(:) = 0.d0
-      w(pt2_stoch_istate) = 1.d0
-      call update_psi_average_norm_contrib(w)
+      state_average_weight(:) = 0.d0
+      state_average_weight(pt2_stoch_istate) = 1.d0
+      TOUCH state_average_weight 
       
       allocate(pt2_detail(N_states,N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
       sumabove = 0d0
       sum2above = 0d0
       Nabove = 0d0
       
-      provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors
+      provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors 
       
       computed = .false.
       
@@ -141,7 +142,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
       
       deallocate(pt2_detail, comb, computed, tbc)
     enddo
-    FREE psi_average_norm_contrib pt2_stoch_istate
+    FREE pt2_stoch_istate 
+    state_average_weight(:) = state_average_weight_save(:)
+    TOUCH state_average_weight
   endif
   do k=N_det+1,N_states
     pt2(k) = 0.d0
diff --git a/plugins/Generators_CAS/generators.irp.f b/plugins/Generators_CAS/generators.irp.f
index 4e2fcd58..259af99d 100644
--- a/plugins/Generators_CAS/generators.irp.f
+++ b/plugins/Generators_CAS/generators.irp.f
@@ -7,7 +7,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
   END_DOC
   integer                        :: i,k,l
   logical                        :: good
-  call write_time(output_determinants)
+  call write_time(6)
   N_det_generators = 0
   do i=1,N_det
     do l=1,n_cas_bitmask
@@ -28,7 +28,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
     endif
   enddo
   N_det_generators = max(N_det_generators,1)
-  call write_int(output_determinants,N_det_generators,'Number of generators')
+  call write_int(6,N_det_generators,'Number of generators')
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
diff --git a/plugins/Generators_full/generators.irp.f b/plugins/Generators_full/generators.irp.f
index eea5821b..4f2c715e 100644
--- a/plugins/Generators_full/generators.irp.f
+++ b/plugins/Generators_full/generators.irp.f
@@ -8,7 +8,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  END_DOC
  integer :: i
  double precision :: norm
- call write_time(output_determinants)
+ call write_time(6)
  norm = 0.d0
  N_det_generators = N_det
  do i=1,N_det
@@ -19,7 +19,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
    endif
  enddo
  N_det_generators = max(N_det_generators,1)
- call write_int(output_determinants,N_det_generators,'Number of generators')
+ call write_int(6,N_det_generators,'Number of generators')
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
@@ -30,15 +30,8 @@ END_PROVIDER
  ! Hartree-Fock determinant
  END_DOC
  integer                        :: i, k
- psi_coef_generators = 0.d0
- psi_det_generators  = 0_bit_kind
- do i=1,N_det_generators
-   do k=1,N_int
-     psi_det_generators(k,1,i) = psi_det_sorted(k,1,i)
-     psi_det_generators(k,2,i) = psi_det_sorted(k,2,i)
-   enddo
-   psi_coef_generators(i,:) = psi_coef_sorted(i,:)
- enddo
+ psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
+ psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
 
 END_PROVIDER
 
diff --git a/plugins/Generators_restart/generators.irp.f b/plugins/Generators_restart/generators.irp.f
index 17854330..bcd8d0d2 100644
--- a/plugins/Generators_restart/generators.irp.f
+++ b/plugins/Generators_restart/generators.irp.f
@@ -14,7 +14,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  else
   print*,'PB in generators restart !!!'
  endif
- call write_int(output_determinants,N_det_generators,'Number of generators')
+ call write_int(6,N_det_generators,'Number of generators')
 END_PROVIDER
 
 
diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
index 241721ae..c66c8985 100644
--- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
+++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
@@ -14,19 +14,21 @@ END_DOC
  
   integer                        :: i,j
   double precision, allocatable :: mo_coef_save(:,:)
+
+  PROVIDE ao_md5 mo_occ level_shift 
  
   allocate(mo_coef_save(ao_num,mo_tot_num),                          &
       Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS),                 &
       error_matrix_DIIS(ao_num,ao_num,max_dim_DIIS)                  &
       )
 
-  call write_time(output_hartree_fock)
+  call write_time(6)
 
-  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+  write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
     '====','================','================','================'
-  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+  write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
     '  N ', 'Energy  ', 'Energy diff  ', 'DIIS error  '
-  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+  write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
     '====','================','================','================'
 
 ! Initialize energies and density matrices
@@ -115,7 +117,7 @@ END_DOC
 
 !   Print results at the end of each iteration
 
-    write(output_hartree_fock,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)')  &
+    write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)')  &
       iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, dim_DIIS
 
     if (Delta_energy_SCF < 0.d0) then
@@ -128,18 +130,18 @@ END_DOC
 ! End of Main SCF loop
 !
 
-  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)') &
+  write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)') &
     '====','================','================','================'
-  write(output_hartree_fock,*)
+  write(6,*)
   
   if(.not.no_oa_or_av_opt)then
    call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
-  call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy')
+  call write_double(6, Energy_SCF, 'Hartree-Fock energy')
   call ezfio_set_hartree_fock_energy(Energy_SCF)
 
-  call write_time(output_hartree_fock)
+  call write_time(6)
 
 end
 
diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f
index 3d71d826..75a07145 100644
--- a/plugins/Hartree_Fock/SCF.irp.f
+++ b/plugins/Hartree_Fock/SCF.irp.f
@@ -23,7 +23,7 @@ subroutine create_guess
       mo_coef = ao_ortho_lowdin_coef
       TOUCH mo_coef
       mo_label = 'Guess'
-      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
+      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,1,.false.)
       SOFT_TOUCH mo_coef mo_label
     else if (mo_guess_type == "Huckel") then
       call huckel_guess
diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f
index 20a8abd7..f97fbf82 100644
--- a/plugins/Hartree_Fock/damping_SCF.irp.f
+++ b/plugins/Hartree_Fock/damping_SCF.irp.f
@@ -28,13 +28,13 @@ subroutine damping_SCF
   enddo
   
   
-  call write_time(output_hartree_fock)
+  call write_time(6)
 
-  write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
+  write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
     '====','================','================','================', '===='
-  write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
+  write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
     '  N ', 'Energy  ', 'Energy diff  ', 'Density diff  ', 'Save'
-  write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
+  write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  &
     '====','================','================','================', '===='
 
   E = HF_energy + 1.d0
@@ -58,7 +58,7 @@ subroutine damping_SCF
       save_char = ' '
     endif
 
-    write(output_hartree_fock,'(I4,1X,F16.10, 1X, F16.10, 1X, F16.10, 3X, A )')  &
+    write(6,'(I4,1X,F16.10, 1X, F16.10, 1X, F16.10, 3X, A )')  &
       k, E, delta_E, delta_D, save_char
     
     D_alpha = HF_density_matrix_ao_alpha
@@ -115,17 +115,17 @@ subroutine damping_SCF
     TOUCH mo_coef
 
   enddo
-  write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  '====','================','================','================', '===='
-  write(output_hartree_fock,*)
+  write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )')  '====','================','================','================', '===='
+  write(6,*)
   
   if(.not.no_oa_or_av_opt)then
    call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
   endif
 
-  call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
+  call write_double(6, E_min, 'Hartree-Fock energy')
   call ezfio_set_hartree_fock_energy(E_min)
 
-  call write_time(output_hartree_fock)
+  call write_time(6)
 
   deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta)
 end
diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index 46b08de2..6609790b 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -211,7 +211,7 @@ END_PROVIDER
        call davidson_diag_mrcc_HS2(psi_det,eigenvectors,             &
            size(eigenvectors,1),                                     &
            eigenvalues,N_det,N_states,N_states_diag,N_int,           &
-           output_determinants,mrcc_state)
+           6,mrcc_state)
        CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
        CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
      enddo
@@ -316,12 +316,12 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
   
   integer                        :: j
   character*(8)                  :: st
-  call write_time(output_determinants)
+  call write_time(6)
   do j=1,min(N_det,N_states)
     write(st,'(I4)') j
     CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
-    call write_double(output_determinants,CI_energy_dressed(j),'Energy of state '//trim(st))
-    call write_double(output_determinants,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
+    call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st))
+    call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
   enddo
 
 END_PROVIDER
diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f
index 34d26127..e186116d 100644
--- a/plugins/MRPT_Utils/mrpt_utils.irp.f
+++ b/plugins/MRPT_Utils/mrpt_utils.irp.f
@@ -354,12 +354,12 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states) ]
   
   integer                        :: j
   character*(8)                  :: st
-  call write_time(output_determinants)
+  call write_time(6)
   do j=1,N_states
     CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
     write(st,'(I4)') j
-    call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))
-    call write_double(output_determinants,CI_dressed_pt2_new_eigenvectors_s2(j),'S^2 of state '//trim(st))
+    call write_double(6,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))
+    call write_double(6,CI_dressed_pt2_new_eigenvectors_s2(j),'S^2 of state '//trim(st))
   enddo
 
 END_PROVIDER
diff --git a/plugins/Perturbation/selection.irp.f b/plugins/Perturbation/selection.irp.f
index a9a5b3bb..6ec5062b 100644
--- a/plugins/Perturbation/selection.irp.f
+++ b/plugins/Perturbation/selection.irp.f
@@ -125,7 +125,7 @@ subroutine remove_small_contributions
   if (N_removed > 0) then
     N_det = N_det - N_removed
     SOFT_TOUCH N_det psi_det psi_coef 
-    call write_int(output_determinants,N_removed, 'Removed determinants')
+    call write_int(6,N_removed, 'Removed determinants')
   endif
 end
 
diff --git a/plugins/Selectors_Utils/zmq.irp.f b/plugins/Selectors_Utils/zmq.irp.f
index 5f40cd4f..b32436aa 100644
--- a/plugins/Selectors_Utils/zmq.irp.f
+++ b/plugins/Selectors_Utils/zmq.irp.f
@@ -45,6 +45,7 @@ integer function zmq_get_$X(zmq_to_qp_run_socket, worker_id)
   integer                        :: rc
   character*(256)                :: msg
 
+  PROVIDE zmq_state
   zmq_get_$X = 0
   if (mpi_master) then
 
diff --git a/plugins/Selectors_full/selectors.irp.f b/plugins/Selectors_full/selectors.irp.f
index 42e3c87b..3d58bdcc 100644
--- a/plugins/Selectors_full/selectors.irp.f
+++ b/plugins/Selectors_full/selectors.irp.f
@@ -8,7 +8,7 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
   END_DOC
   integer                        :: i
   double precision               :: norm, norm_max
-  call write_time(output_determinants)
+  call write_time(6)
   N_det_selectors = N_det
   if (threshold_generators < 1.d0) then
     norm = 0.d0
@@ -21,7 +21,7 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
     enddo
     N_det_selectors = max(N_det_selectors,N_det_generators)
   endif
-  call write_int(output_determinants,N_det_selectors,'Number of selectors')
+  call write_int(6,N_det_selectors,'Number of selectors')
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
diff --git a/plugins/Selectors_no_sorted/selectors.irp.f b/plugins/Selectors_no_sorted/selectors.irp.f
index 3ac8218d..e81aa795 100644
--- a/plugins/Selectors_no_sorted/selectors.irp.f
+++ b/plugins/Selectors_no_sorted/selectors.irp.f
@@ -8,11 +8,11 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
  END_DOC
  integer :: i
  double precision :: norm
- call write_time(output_determinants)
+ call write_time(6)
  norm = 0.d0
  N_det_selectors = N_det
  N_det_selectors = max(N_det_selectors,1)
- call write_int(output_determinants,N_det_selectors,'Number of selectors')
+ call write_int(6,N_det_selectors,'Number of selectors')
 END_PROVIDER
 
 
diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f
index 9540295f..ffce8082 100644
--- a/plugins/Symmetry/Symmetry.main.irp.f
+++ b/plugins/Symmetry/Symmetry.main.irp.f
@@ -3,9 +3,17 @@ program Symmetry
   BEGIN_DOC
 ! TODO
   END_DOC
-  integer :: i, j
+  integer :: i, j, k
   character*8 :: sym
 
+do k=1,n_irrep
+  print *,  sym_operation(k)
+  do i=1,mo_tot_num
+    print '(1000(F8.4,X))', mo_symm(i,:,k), sum(mo_symm(i,:,k))
+  enddo
+  print *,  ''
+enddo
+
   print *,  'Molecule is linear:  ', molecule_is_linear
   print *,  'Has center of inversion:  ', molecule_has_center_of_inversion
   print *,  'Has S2n improper rotation:  ', molecule_has_improper_rotation
@@ -17,5 +25,4 @@ program Symmetry
   do i=1,n_irrep
     print *,  i, real(character_table(i,:))
   enddo
-  PROVIDE mo_sym
 end
diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f
index 1ed567bc..ed746a40 100644
--- a/plugins/Symmetry/aos.irp.f
+++ b/plugins/Symmetry/aos.irp.f
@@ -22,15 +22,19 @@ subroutine generate_sym_coord(n_sym_points,result)
   BEGIN_DOC
   ! xyz coordinates of points to check the symmetry, drawn uniformly in the molecular box.
   END_DOC
-  integer                        :: i, xyz
+  integer                        :: i, iop
   
-  do i=1,n_sym_points
-    call random_number(result(1,i))
-    call random_number(result(2,i))
-    call random_number(result(3,i))
-  enddo
-  do xyz=1,3
-    result(xyz,1:n_sym_points) = sym_box(xyz,1) + result(xyz,:) * (sym_box(xyz,2)-sym_box(xyz,1))
+  double precision, external :: halton_ranf
+  do i=1,n_sym_points,n_irrep
+    result(1,i) = sym_box(1,1) + halton_ranf(1) * (sym_box(1,2)-sym_box(1,1))
+    result(2,i) = sym_box(1,1) + halton_ranf(2) * (sym_box(2,2)-sym_box(2,1))
+    result(3,i) = sym_box(1,1) + halton_ranf(3) * (sym_box(3,2)-sym_box(3,1))
+    do iop=2,n_irrep
+      if (iop-1+i > n_sym_points) exit
+      call dgemm('N','N',3,1,3,1.d0,sym_transformation_matrices(1,1,iop), &
+        size(sym_transformation_matrices,1),&
+        result(1,i),size(result,1),0.d0,result(1,i+iop-1),size(result,1))
+    enddo
   enddo
   
 end
diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f
index 38c0a6b7..817638b4 100644
--- a/plugins/Symmetry/find_sym.irp.f
+++ b/plugins/Symmetry/find_sym.irp.f
@@ -47,6 +47,7 @@ BEGIN_PROVIDER [ integer, sym_rotation_axis, (3) ]
   logical                        :: found
   double precision, external     :: u_dot_u
   integer                        :: iorder, iaxis
+
   do iaxis=1,3
     do iorder=12,2,-1
       sym_rotation_axis(iaxis) = iorder
@@ -300,14 +301,6 @@ BEGIN_PROVIDER [ character*16, point_group ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ character*8, mo_symmetry ]
- implicit none
- BEGIN_DOC
- ! Symmetry of the MOs
- END_DOC
- integer :: i,j
-END_PROVIDER
-
 
 BEGIN_PROVIDER [ integer, n_irrep ]
  implicit none
@@ -369,105 +362,6 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ]
   double precision               :: sym_operations_on_mos(mo_tot_num)
   logical                        :: possible_irrep(n_irrep,mo_tot_num)
 
-  n_sym_points = 10
-  allocate(val(n_sym_points,mo_tot_num,2), sym_points(3,n_sym_points), ref_points(3,n_sym_points))
-
-  call generate_sym_coord(n_sym_points,ref_points)
-  call compute_sym_mo_values(ref_points,n_sym_points,val(1,1,2))
-
-  possible_irrep = .True.
-  do iop=1,n_irrep
-    if (sym_operation(iop) == 'E') then
-      cycle
-    endif
-
-    if (sym_operation(iop) == 'i') then
-      do ipoint=1,n_sym_points
-        call sym_apply_inversion(ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 'sh') then
-      do ipoint=1,n_sym_points
-        call sym_apply_reflexion(molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 's') then
-      do ipoint=1,n_sym_points
-        call sym_apply_reflexion(molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 'sv') then 
-      do ipoint=1,n_sym_points
-        call sym_apply_reflexion(molecule_ternary_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 'sd') then
-      angle = dble(maxval(sym_rotation_axis))
-      do ipoint=1,n_sym_points
-        call sym_apply_diagonal_reflexion(angle,molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 'C2''') then 
-      angle = 2.d0
-      do ipoint=1,n_sym_points
-        call sym_apply_rotation(angle,molecule_secondary_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else if (sym_operation(iop) == 'C2"') then 
-      angle = 2.d0
-      do ipoint=1,n_sym_points
-        call sym_apply_rotation(angle,molecule_ternary_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-      enddo
-    else
-      do l=2,len(sym_operation(iop))
-        if (sym_operation(iop)(l:l) == '^') exit
-      enddo
-      read(sym_operation(iop)(2:l-1), *) iangle
-      if (l == len(sym_operation(iop))+1) then
-        l=1
-      else
-        read(sym_operation(iop)(l+1:), *, err=10, end=10) l
-        10 continue
-      endif
-      angle = dble(iangle)/(dble(l))
-      if (sym_operation(iop)(1:1) == 'C') then 
-        do ipoint=1,n_sym_points
-          call sym_apply_rotation(angle,molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-        enddo
-      else if (sym_operation(iop)(1:1) == 'S') then 
-        do ipoint=1,n_sym_points
-          call sym_apply_improper_rotation(angle,molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint))
-        enddo
-      endif
-    endif
-
-    call compute_sym_mo_values(sym_points,n_sym_points,val(1,1,1))
-
-    print *,  sym_operation(iop)
-    double precision :: icount
-    do imo=1,mo_tot_num
-      sym_operations_on_mos(imo) = 0.d0
-      icount = 0
-      do ipoint=1,n_sym_points
-        double precision :: x
-        if (dabs(val(ipoint,imo,1)) < 1.d-5) cycle
-        icount += 1.d0 
-        x =  val(ipoint,imo,1)/val(ipoint,imo,2)
-        if (dabs(x) > 1.d0) then
-          x = 1.d0/x
-        endif
-        sym_operations_on_mos(imo) += x
-      enddo
-      sym_operations_on_mos(imo) *= 1.d0/icount
-      if (dabs(sym_operations_on_mos(imo) - 1.d0) < 1.d-2) then
-        sym_operations_on_mos(imo) = 1.d0
-      else if (dabs(sym_operations_on_mos(imo) + 1.d0) < 1.d-2) then
-        sym_operations_on_mos(imo) = -1.d0
-      else if (dabs(sym_operations_on_mos(imo)) < 1.d-2) then
-        sym_operations_on_mos(imo) = 0.d0
-      endif
-      print *, imo,  sym_operations_on_mos(imo)
-      do i=1,n_irrep
-        if (dabs(character_table(i,iop) - sym_operations_on_mos(imo)) > 1.d-2) then
-          possible_irrep(i,imo) = .False.
-        endif
-      enddo
-    enddo
-  enddo
   do imo=1,mo_tot_num
     print *,  'MO ', imo
     do i=1,n_irrep
diff --git a/plugins/Symmetry/nuclei.irp.f b/plugins/Symmetry/nuclei.irp.f
index d680393d..405b529a 100644
--- a/plugins/Symmetry/nuclei.irp.f
+++ b/plugins/Symmetry/nuclei.irp.f
@@ -69,26 +69,26 @@ BEGIN_PROVIDER [ double precision, nucl_coord_sym,  (nucl_num,3) ]
      character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
      double precision, parameter    :: a0= 0.529177249d0
      
-     call write_time(output_Nuclei)
-     write(output_Nuclei,'(A)') ''
-     write(output_Nuclei,'(A)') 'Nuclear Coordinates in standard orientation (Angstroms)'
-     write(output_Nuclei,'(A)') '======================================================='
-     write(output_Nuclei,'(A)') ''
-     write(output_Nuclei,ft)                                         &
+     call write_time(6)
+     write(6,'(A)') ''
+     write(6,'(A)') 'Nuclear Coordinates in standard orientation (Angstroms)'
+     write(6,'(A)') '======================================================='
+     write(6,'(A)') ''
+     write(6,ft)                                         &
          '================','============','============','============','============'
-     write(output_Nuclei,*)                                          &
+     write(6,*)                                          &
          '     Atom          Charge          X            Y            Z '
-     write(output_Nuclei,ft)                                         &
+     write(6,ft)                                         &
          '================','============','============','============','============'
      do i=1,nucl_num
-       write(output_Nuclei,f) nucl_label(i), nucl_charge(i),         &
+       write(6,f) nucl_label(i), nucl_charge(i),         &
            nucl_coord_sym(i,1)*a0,                                       &
            nucl_coord_sym(i,2)*a0,                                       &
            nucl_coord_sym(i,3)*a0
      enddo
-     write(output_Nuclei,ft)                                         &
+     write(6,ft)                                         &
          '================','============','============','============','============'
-     write(output_Nuclei,'(A)') ''
+     write(6,'(A)') ''
      
    endif
    
diff --git a/plugins/Symmetry/sym_operation.irp.f b/plugins/Symmetry/sym_operation.irp.f
index cfc86621..ccf72ec3 100644
--- a/plugins/Symmetry/sym_operation.irp.f
+++ b/plugins/Symmetry/sym_operation.irp.f
@@ -23,7 +23,7 @@ subroutine sym_apply_diagonal_reflexion(angle,iaxis,point_in,point_out)
   double precision :: point_tmp1(3), point_tmp2(3)
   integer :: iaxis2
   iaxis2 = mod(iaxis,3)+1
-  iaxis2 = mod(iaxis2,3)+1
+!  iaxis2 = mod(iaxis2,3)+1
   call sym_apply_rotation(-angle,iaxis,point_in,point_tmp1)
   call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2)
   call sym_apply_rotation(angle,iaxis,point_tmp2,point_out)
diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f
index 5dfa8556..727bdba7 100644
--- a/plugins/mrcepa0/dressing.irp.f
+++ b/plugins/mrcepa0/dressing.irp.f
@@ -963,6 +963,19 @@ END_PROVIDER
         enddo
       end do
     end do
+!  else if(mrmode == 10) then
+!    do i = 1, N_det_ref
+!      do i_state = 1, N_states
+!        delta_ii(i_state,i)= delta_ii_mrsc2(i_state,i)
+!        delta_ii_s2(i_state,i)= delta_ii_s2_mrsc2(i_state,i)
+!      enddo
+!      do j = 1, N_det_non_ref
+!        do i_state = 1, N_states
+!          delta_ij(i_state,j,i) = delta_ij_mrsc2(i_state,j,i)
+!          delta_ij_s2(i_state,j,i) = delta_ij_s2_mrsc2(i_state,j,i)
+!        enddo
+!      end do
+!    end do
   else if(mrmode == 5) then
     do i = 1, N_det_ref
       do i_state = 1, N_states
diff --git a/plugins/mrcepa0/mrcc_stoch_routines.irp.f b/plugins/mrcepa0/mrcc_stoch_routines.irp.f
index 78940d5e..d68b5137 100644
--- a/plugins/mrcepa0/mrcc_stoch_routines.irp.f
+++ b/plugins/mrcepa0/mrcc_stoch_routines.irp.f
@@ -32,17 +32,15 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error)
   
   double precision, external     :: omp_get_wtime
   double precision               :: time
-  double precision               :: w(N_states)
 
 
+  state_average_weight(:) = 0.d0
+  state_average_weight(mrcc_stoch_istate) = 1.d0
+  TOUCH state_average_weight 
+  
 
   provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral mrcc_weight psi_selectors
 
-  w(:) = 0.d0
-  w(mrcc_stoch_istate) = 1.d0
-  call update_psi_average_norm_contrib(w)
-
-  
   
   
   print *, '========== ================= ================= ================='
diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py
index ee44a1e1..8d154fc2 100755
--- a/scripts/ezfio_interface/ei_handler.py
+++ b/scripts/ezfio_interface/ei_handler.py
@@ -329,7 +329,8 @@ def create_ezfio_provider(dict_ezfio_cfg):
             ez_p.set_doc(dict_info['doc'])
             ez_p.set_ezfio_dir(dict_info['ezfio_dir'])
             ez_p.set_ezfio_name(dict_info['ezfio_name'])
-            ez_p.set_output("output_%s" % dict_info['module'].lower)
+            ez_p.set_output("6")
+#            ez_p.set_output("output_%s" % dict_info['module'].lower)
 
             # (nuclei.nucl_num,pseudo.klocmax) => (nucl_num,klocmax)
             ez_p.set_size(re.sub(r'\w+\.', "", dict_info['size']))
diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py
index da1aac5a..3b5d96c2 100755
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@@ -142,7 +142,7 @@ class H_apply(object):
   endif
   SOFT_TOUCH psi_det psi_coef N_det
 """
-    s["printout_now"]   = """write(output_determinants,*)  &
+    s["printout_now"]   = """write(6,*)  &
        100.*float(i_generator)/float(N_det_generators), '% in ', wall_1-wall_0, 's'"""
     self.data = s
 
@@ -370,9 +370,9 @@ class H_apply(object):
     delta_pt2(k) = 0.d0
     pt2_old(k) = 0.d0
   enddo
-        write(output_determinants,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
+        write(6,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
                  'N_generators', 'Norm', 'Delta PT2', 'PT2', 'Est. PT2', 'secs'
-        write(output_determinants,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
+        write(6,'(A12, 1X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
                  '============', '========', '=========', '=========', '=========', &
                  '========='
       """ 
@@ -385,7 +385,7 @@ class H_apply(object):
       """
       self.data["printout_now"] = """
       do k=1,N_st
-        write(output_determinants,'(I10, 4(2X, F9.6), 2X, F8.1)') &
+        write(6,'(I10, 4(2X, F9.6), 2X, F8.1)') &
                  i_generator, norm_psi(k), delta_pt2(k), pt2(k), &
                  pt2(k)/(norm_psi(k)*norm_psi(k)), &
                  wall_1-wall_0
@@ -416,7 +416,7 @@ class H_apply(object):
       SOFT_TOUCH psi_det psi_coef N_det
       selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0
       selection_criterion = selection_criterion_min
-      call write_double(output_determinants,selection_criterion,'Selection criterion')
+      call write_double(6,selection_criterion,'Selection criterion')
       """
       self.data["keys_work"] = """
       e_2_pert_buffer = 0.d0
diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f
index 5f3fc7f5..4eac3e0a 100644
--- a/src/Bitmask/bitmasks.irp.f
+++ b/src/Bitmask/bitmasks.irp.f
@@ -100,7 +100,7 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ]
  ! Number of bitmasks for generators
  END_DOC
  logical                        :: exists
- PROVIDE ezfio_filename
+ PROVIDE ezfio_filename N_int
  
  if (mpi_master) then
   call ezfio_has_bitmasks_N_mask_gen(exists)
diff --git a/src/Davidson/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f
index f01cfb28..c65d9763 100644
--- a/src/Davidson/diagonalize_CI.irp.f
+++ b/src/Davidson/diagonalize_CI.irp.f
@@ -7,14 +7,14 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
   
   integer                        :: j
   character*(8)                  :: st
-  call write_time(output_determinants)
+  call write_time(6)
   do j=1,min(N_det,N_states_diag)
     CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
   enddo
   do j=1,min(N_det,N_states)
     write(st,'(I4)') j
-    call write_double(output_determinants,CI_energy(j),'Energy of state '//trim(st))
-    call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
+    call write_double(6,CI_energy(j),'Energy of state '//trim(st))
+    call write_double(6,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
   enddo
 
 END_PROVIDER
@@ -58,14 +58,14 @@ END_PROVIDER
      
 !     call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy,  &
 !         size(CI_eigenvectors,1),               &
-!         N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
+!         N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,6)
 !
 !     call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
 !         min(N_det,N_states_diag),size(CI_eigenvectors,1))
 
      call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, &
          size(CI_eigenvectors,1),CI_electronic_energy,               &
-         N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
+         N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,6)
      
    else if (diag_algorithm == "Lapack") then
      
diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f
index ef396f9c..e5197a21 100644
--- a/src/Determinants/H_apply.irp.f
+++ b/src/Determinants/H_apply.irp.f
@@ -265,7 +265,7 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
       endif
     enddo
     N_det = k
-    call write_bool(output_determinants,found_duplicates,'Found duplicate determinants')
+    call write_bool(6,found_duplicates,'Found duplicate determinants')
     SOFT_TOUCH N_det psi_det psi_coef
   endif
   deallocate (duplicate,bit_tmp)
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index be28183b..bd5f0741 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -368,13 +368,13 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
  END_DOC
  logical :: exists
 
- state_average_weight = 1.d0
+ state_average_weight(:) = 1.d0
  call ezfio_has_determinants_state_average_weight(exists)
  if (exists) then
   call ezfio_get_determinants_state_average_weight(state_average_weight)
  endif
- state_average_weight = state_average_weight+1.d-31
- state_average_weight = state_average_weight/(sum(state_average_weight))
+ state_average_weight(:) = state_average_weight(:)+1.d-31
+ state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:)))
 END_PROVIDER
 
 
diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f
index 40c4524a..8530fa64 100644
--- a/src/Determinants/determinants.irp.f
+++ b/src/Determinants/determinants.irp.f
@@ -43,7 +43,7 @@ BEGIN_PROVIDER [ integer, N_det ]
     else
       N_det = 1
     endif
-    call write_int(output_determinants,N_det,'Number of determinants')
+    call write_int(6,N_det,'Number of determinants')
   endif
   IRP_IF MPI
     include 'mpif.h'
@@ -77,7 +77,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
   BEGIN_DOC
   ! Size of the psi_det/psi_coef arrays
   END_DOC
-  PROVIDE ezfio_filename output_determinants
+  PROVIDE ezfio_filename 
   logical                        :: exists
   if (mpi_master) then
     call ezfio_has_determinants_n_det(exists)
@@ -87,7 +87,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
       psi_det_size = 1
     endif
     psi_det_size = max(psi_det_size,100000)
-    call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays')
+    call write_int(6,psi_det_size,'Dimension of the psi arrays')
   endif
   IRP_IF MPI
     include 'mpif.h'
@@ -225,34 +225,6 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
   
 END_PROVIDER
 
-subroutine update_psi_average_norm_contrib(w)
-  implicit none
-  BEGIN_DOC
-  ! Compute psi_average_norm_contrib for different state average weights w(:)
-  END_DOC
-  double precision, intent(in)   :: w(N_states)
-  double precision               :: w0(N_states), f
-  w0(:) = w(:)/sum(w(:))
-  
-  integer                        :: i,j,k
-  do i=1,N_det
-    psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*w(1)
-  enddo
-  do k=2,N_states
-    do i=1,N_det
-      psi_average_norm_contrib(i) = psi_average_norm_contrib(i) +    &
-          psi_coef(i,k)*psi_coef(i,k)*w(k)
-    enddo
-  enddo
-  f = 1.d0/sum(psi_average_norm_contrib(1:N_det))
-  do i=1,N_det
-    psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f
-  enddo
-  SOFT_TOUCH psi_average_norm_contrib
-  
-end subroutine
-
-
 BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
   implicit none
   BEGIN_DOC
@@ -260,14 +232,12 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
   END_DOC
   integer                        :: i,j,k
   double precision               :: f
-  f = 1.d0/dble(N_states)
-  do i=1,N_det
-    psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f
-  enddo
-  do k=2,N_states
+
+  psi_average_norm_contrib(:) = 0.d0
+  do k=1,N_states
     do i=1,N_det
       psi_average_norm_contrib(i) = psi_average_norm_contrib(i) +    &
-          psi_coef(i,k)*psi_coef(i,k)*f
+          psi_coef(i,k)*psi_coef(i,k)*state_average_weight(k)
     enddo
   enddo
   f = 1.d0/sum(psi_average_norm_contrib(1:N_det))
@@ -562,7 +532,7 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
     
     call ezfio_set_determinants_psi_coef(psi_coef_save)
     deallocate (psi_coef_save)
-    call write_int(output_determinants,ndet,'Saved determinants')
+    call write_int(6,ndet,'Saved determinants')
   endif
 end
 
@@ -634,7 +604,7 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde
   enddo
   
   call ezfio_set_determinants_psi_coef(psi_coef_save)
-  call write_int(output_determinants,ndet,'Saved determinants')
+  call write_int(6,ndet,'Saved determinants')
   deallocate (psi_coef_save)
 end
 
diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f
index 968ced57..591843f7 100644
--- a/src/Determinants/psi_cas.irp.f
+++ b/src/Determinants/psi_cas.irp.f
@@ -41,7 +41,7 @@ use bitmasks
       enddo
     endif
   enddo
-  call write_int(output_determinants,N_det_cas, 'Number of determinants in the CAS')
+  call write_int(6,N_det_cas, 'Number of determinants in the CAS')
 
 END_PROVIDER
 
diff --git a/src/Determinants/zmq.irp.f b/src/Determinants/zmq.irp.f
index b0ef55d6..5e42cff5 100644
--- a/src/Determinants/zmq.irp.f
+++ b/src/Determinants/zmq.irp.f
@@ -87,6 +87,7 @@ integer function zmq_get_$X(zmq_to_qp_run_socket, worker_id)
   integer                        :: rc
   character*(256)                :: msg
 
+  PROVIDE zmq_state
   zmq_get_$X = 0
   if (mpi_master) then
     write(msg,'(A,1X,I8,1X,A200)') 'get_data '//trim(zmq_state), worker_id, '$X'
@@ -272,6 +273,7 @@ integer function zmq_get_psi_det(zmq_to_qp_run_socket, worker_id)
   integer*8                      :: rc8
   character*(256)                :: msg
 
+  PROVIDE zmq_state
   zmq_get_psi_det = 0
   if (mpi_master) then
     write(msg,'(A,1X,I8,1X,A200)') 'get_data '//trim(zmq_state), worker_id, 'psi_det'
@@ -319,6 +321,7 @@ integer function zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id)
   integer*8                      :: rc8
   character*(256)                :: msg
 
+  PROVIDE zmq_state psi_det_size
   zmq_get_psi_coef = 0
   if (mpi_master) then
     write(msg,'(A,1X,I8,1X,A200)') 'get_data '//trim(zmq_state), worker_id, 'psi_coef'
diff --git a/src/Ezfio_files/output.irp.f b/src/Ezfio_files/output.irp.f
index 25f862bd..c8e50d48 100644
--- a/src/Ezfio_files/output.irp.f
+++ b/src/Ezfio_files/output.irp.f
@@ -8,36 +8,6 @@
   call wall_time(output_wall_time_0)
 END_PROVIDER
  
-BEGIN_SHELL [ /bin/bash ]
- 
- for NAME in $(\ls -d ${QP_ROOT}/src/*/)
- do
- NAME=$(basename ${NAME})
- cat << EOF
- BEGIN_PROVIDER [ integer, output_$NAME ]
-  implicit none
-  BEGIN_DOC
-  ! Output file for $NAME
-  END_DOC
-  PROVIDE output_wall_time_0 output_cpu_time_0 ezfio_filename
-!  integer                        :: getUnitAndOpen
-!  call ezfio_set_output_empty(.False.)
-IRP_IF COARRAY
-  if (this_image() == 1) then
-    output_$NAME = 6 !getUnitAndOpen(trim(ezfio_filename)//'/output/'//'$NAME.rst','a')
-  else
-    output_$NAME = getUnitAndOpen('/dev/null','w')
-  endif
-IRP_ELSE
-  output_$NAME = 6 !getUnitAndOpen(trim(ezfio_filename)//'/output/'//'$NAME.rst','a')
-IRP_ENDIF
-  write(output_$NAME,'(A)')                                          &
-      '--------------------------------------------------------------------------------'
- END_PROVIDER
-EOF
- done
- 
-END_SHELL
  
 subroutine write_time(iunit)
   implicit none
diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f
index c39dba36..9d884456 100644
--- a/src/Integrals_Bielec/mo_bi_integrals.irp.f
+++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f
@@ -122,17 +122,17 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
     endif
     
   else
-!    call add_integrals_to_map(full_ijkl_bitmask_4)
+    call add_integrals_to_map(full_ijkl_bitmask_4)
 
 !     call four_index_transform_zmq(ao_integrals_map,mo_integrals_map, &
 !         mo_coef, size(mo_coef,1),                                      &
 !         1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num,                    &
 !         1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
 !
-     call four_index_transform_block(ao_integrals_map,mo_integrals_map, &
-         mo_coef, size(mo_coef,1),                                      &
-         1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num,                    &
-         1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
+!     call four_index_transform_block(ao_integrals_map,mo_integrals_map, &
+!         mo_coef, size(mo_coef,1),                                      &
+!         1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num,                    &
+!         1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
 !
 !     call four_index_transform(ao_integrals_map,mo_integrals_map, &
 !         mo_coef, size(mo_coef,1),                                      &
diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f
index b0400f67..5d01efb9 100644
--- a/src/MO_Basis/ao_ortho_canonical.irp.f
+++ b/src/MO_Basis/ao_ortho_canonical.irp.f
@@ -129,7 +129,7 @@ END_PROVIDER
     enddo
 
     ao_ortho_canonical_num = ao_cart_to_sphe_num
-    call ortho_canonical (ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
+    call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
       ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num)
 
     call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, &
diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f
index 4806582b..212d2b6f 100644
--- a/src/MO_Basis/utils.irp.f
+++ b/src/MO_Basis/utils.irp.f
@@ -55,7 +55,7 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output)
   double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
   
-  call write_time(output_mo_basis)
+  call write_time(6)
   if (m /= mo_tot_num) then
     print *, irp_here, ': Error : m/= mo_tot_num'
     stop 1
@@ -78,12 +78,12 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output)
   
   call lapack_diag(eigvalues,R,A,n,m)
   if (output) then
-    write (output_mo_basis,'(A)')  'MOs are now **'//trim(label)//'**'
-    write (output_mo_basis,'(A)') ''
-    write (output_mo_basis,'(A)')  'Eigenvalues'
-    write (output_mo_basis,'(A)') '-----------'
-    write (output_mo_basis,'(A)')  ''
-    write (output_mo_basis,'(A)') '======== ================'
+    write (6,'(A)')  'MOs are now **'//trim(label)//'**'
+    write (6,'(A)') ''
+    write (6,'(A)')  'Eigenvalues'
+    write (6,'(A)') '-----------'
+    write (6,'(A)')  ''
+    write (6,'(A)') '======== ================'
   endif
   if (sign == -1) then
     do i=1,m
@@ -92,15 +92,15 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output)
   endif
   if (output) then
     do i=1,m
-      write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
+      write (6,'(I8,1X,F16.10)')  i,eigvalues(i)
     enddo
-    write (output_mo_basis,'(A)') '======== ================'
-    write (output_mo_basis,'(A)')  ''
+    write (6,'(A)') '======== ================'
+    write (6,'(A)')  ''
   endif
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
   deallocate(A,mo_coef_new,R,eigvalues)
-  call write_time(output_mo_basis)
+  call write_time(6)
   
   mo_label = label
 end
@@ -115,7 +115,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   double precision, allocatable  :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
   
-  call write_time(output_mo_basis)
+  call write_time(6)
   if (m /= mo_tot_num) then
     print *, irp_here, ': Error : m/= mo_tot_num'
     stop 1
@@ -132,22 +132,22 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   
   call svd(A,lda,U,lda,D,Vt,lda,m,n)
 
-  write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**'
-  write (output_mo_basis,'(A)')  ''
-  write (output_mo_basis,'(A)') 'Eigenvalues'
-  write (output_mo_basis,'(A)')  '-----------'
-  write (output_mo_basis,'(A)') ''
-  write (output_mo_basis,'(A)')  '======== ================'
+  write (6,'(A)') 'MOs are now **'//trim(label)//'**'
+  write (6,'(A)')  ''
+  write (6,'(A)') 'Eigenvalues'
+  write (6,'(A)')  '-----------'
+  write (6,'(A)') ''
+  write (6,'(A)')  '======== ================'
 
   do i=1,m
-    write (output_mo_basis,'(I8,1X,F16.10)')  i,D(i)
+    write (6,'(I8,1X,F16.10)')  i,D(i)
   enddo
-  write (output_mo_basis,'(A)')  '======== ================'
-  write (output_mo_basis,'(A)')  ''
+  write (6,'(A)')  '======== ================'
+  write (6,'(A)')  ''
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
   deallocate(A,mo_coef_new,U,Vt,D)
-  call write_time(output_mo_basis)
+  call write_time(6)
   
   mo_label = label
 end
@@ -162,7 +162,7 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,
   integer,allocatable :: iorder(:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
   
-  call write_time(output_mo_basis)
+  call write_time(6)
   if (m /= mo_tot_num) then
     print *, irp_here, ': Error : m/= mo_tot_num'
     stop 1
@@ -213,21 +213,21 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,
    print*,''
   enddo
 
-  write (output_mo_basis,'(A)')  'MOs are now **'//trim(label)//'**'
-  write (output_mo_basis,'(A)')  ''
-  write (output_mo_basis,'(A)')  'Eigenvalues'
-  write (output_mo_basis,'(A)')  '-----------'
-  write (output_mo_basis,'(A)')  ''
-  write (output_mo_basis,'(A)')  '======== ================'
+  write (6,'(A)')  'MOs are now **'//trim(label)//'**'
+  write (6,'(A)')  ''
+  write (6,'(A)')  'Eigenvalues'
+  write (6,'(A)')  '-----------'
+  write (6,'(A)')  ''
+  write (6,'(A)')  '======== ================'
   do i = 1, m
-   write (output_mo_basis,'(I8,1X,F16.10)')  i,eigvalues(i)
+   write (6,'(I8,1X,F16.10)')  i,eigvalues(i)
   enddo
-  write (output_mo_basis,'(A)')  '======== ================'
-  write (output_mo_basis,'(A)')  ''
+  write (6,'(A)')  '======== ================'
+  write (6,'(A)')  ''
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
   deallocate(mo_coef_new,R,eigvalues)
-  call write_time(output_mo_basis)
+  call write_time(6)
   
   mo_label = label
   SOFT_TOUCH mo_coef mo_label
@@ -261,12 +261,12 @@ subroutine mo_sort_by_observable(observable,label)
    enddo
   enddo
 
-  write (output_mo_basis,'(A)')  'MOs are now **'//trim(label)//'**'
-  write (output_mo_basis,'(A)')  ''
+  write (6,'(A)')  'MOs are now **'//trim(label)//'**'
+  write (6,'(A)')  ''
 
   
   deallocate(mo_coef_new,value)
-! call write_time(output_mo_basis)
+! call write_time(6)
   
   mo_label = label
   SOFT_TOUCH mo_coef mo_label
diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f
index 910e9167..3528bf50 100644
--- a/src/Nuclei/nuclei.irp.f
+++ b/src/Nuclei/nuclei.irp.f
@@ -31,26 +31,26 @@ BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
      character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
      double precision, parameter    :: a0= 0.529177249d0
      
-     call write_time(output_Nuclei)
-     write(output_Nuclei,'(A)') ''
-     write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)'
-     write(output_Nuclei,'(A)') '==============================='
-     write(output_Nuclei,'(A)') ''
-     write(output_Nuclei,ft)                                         &
+     call write_time(6)
+     write(6,'(A)') ''
+     write(6,'(A)') 'Nuclear Coordinates (Angstroms)'
+     write(6,'(A)') '==============================='
+     write(6,'(A)') ''
+     write(6,ft)                                         &
          '================','============','============','============','============'
-     write(output_Nuclei,*)                                          &
+     write(6,*)                                          &
          '     Atom          Charge          X            Y            Z '
-     write(output_Nuclei,ft)                                         &
+     write(6,ft)                                         &
          '================','============','============','============','============'
      do i=1,nucl_num
-       write(output_Nuclei,f) nucl_label(i), nucl_charge(i),         &
+       write(6,f) nucl_label(i), nucl_charge(i),         &
            nucl_coord(i,1)*a0,                                       &
            nucl_coord(i,2)*a0,                                       &
            nucl_coord(i,3)*a0
      enddo
-     write(output_Nuclei,ft)                                         &
+     write(6,ft)                                         &
          '================','============','============','============','============'
-     write(output_Nuclei,'(A)') ''
+     write(6,'(A)') ''
      
    endif
    
@@ -189,10 +189,9 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
      enddo
      nuclear_repulsion *= 0.5d0
    end if
-   
-   call write_time(output_Nuclei)
-   call write_double(output_Nuclei,nuclear_repulsion,                &
-       'Nuclear repulsion energy')
+
+   call write_time(6)
+   call write_double(6,nuclear_repulsion,'Nuclear repulsion energy')
    
    if (disk_access_nuclear_repulsion.EQ.'Write') then
      if (mpi_master) then
@@ -203,8 +202,8 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
    
 END_PROVIDER
 
- BEGIN_PROVIDER [ character*(4), element_name, (0:128)]
-&BEGIN_PROVIDER [ double precision, element_mass, (0:128) ]
+ BEGIN_PROVIDER [ character*(4), element_name, (0:127)]
+&BEGIN_PROVIDER [ double precision, element_mass, (0:127) ]
    BEGIN_DOC
    ! Array of the name of element, sorted by nuclear charge (integer)
    END_DOC
@@ -216,7 +215,7 @@ END_PROVIDER
      filename = trim(filename)//'/data/list_element.txt'
      iunit = getUnitAndOpen(filename,'r')
      element_mass(:) = 0.d0
-     do i=0,128
+     do i=0,127
        write(element_name(i),'(I4)') i
      enddo
      character*(80)                 :: buffer, dummy
@@ -232,11 +231,11 @@ END_PROVIDER
  IRP_IF MPI
   include 'mpif.h'
   integer                        :: ierr
-  call MPI_BCAST( element_name, size(element_name)*4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
+  call MPI_BCAST( element_name, 128*4, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
   if (ierr /= MPI_SUCCESS) then
     stop 'Unable to read element_name with MPI'
   endif
-  call MPI_BCAST( element_mass, size(element_mass), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+  call MPI_BCAST( element_mass, 128, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
   if (ierr /= MPI_SUCCESS) then
     stop 'Unable to read element_name with MPI'
   endif
diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f
index 6e1b9565..29d8784f 100644
--- a/src/Utils/LinearAlgebra.irp.f
+++ b/src/Utils/LinearAlgebra.irp.f
@@ -81,10 +81,11 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
 
   call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
 
+  D(:) = dsqrt(D(:))
   m=n
   do i=1,n
     if ( D(i) >= 1.d-6 ) then
-      D(i) = 1.d0/dsqrt(D(i))
+      D(i) = 1.d0/D(i)
     else
       m = i-1
       print *,  'Removed Linear dependencies below:', 1.d0/D(m)
@@ -101,25 +102,19 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
     endif
   enddo
 
-
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(S,U,D,Vt,n,C,m) &
-  !$OMP PRIVATE(i,j)
-
-  !$OMP DO
   do j=1,n
     do i=1,n
       S(i,j) = U(i,j)*D(j)
     enddo
+  enddo
+
+  do j=1,n
     do i=1,n
       U(i,j) = C(i,j)
     enddo
   enddo
-  !$OMP END DO
   
-  !$OMP END PARALLEL
-
-  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
+  call dgemm('N','N',n,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
   deallocate (U, Vt, D, S)
   
 end
diff --git a/src/ZMQ/put_get.irp.f b/src/ZMQ/put_get.irp.f
index 4086b8ed..40d2e881 100644
--- a/src/ZMQ/put_get.irp.f
+++ b/src/ZMQ/put_get.irp.f
@@ -51,6 +51,7 @@ integer function zmq_get_dvector(zmq_to_qp_run_socket, worker_id, name, x, size_
   integer*8                      :: rc8
   character*(256)                :: msg
 
+  PROVIDE zmq_state
   ! Success
   zmq_get_dvector = 0