diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c
index 76b64dd0..3510db37 100644
--- a/src/csf/cfgCI_utils.c
+++ b/src/csf/cfgCI_utils.c
@@ -255,23 +255,24 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
 
 void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
 
-void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
-  int i,j;
-  //for(j=0;j<cols;++j){
-  //  for(i=0;i<rows;++i){
-  //    printf(" %3.2f ",overlapMatrix[j*rows + i]);
-  //  }
-  //  printf("\n");
-  //}
-  // Call the function ortho_qr from qp
-  ortho_qr_csf(overlapMatrix, rows, orthoMatrix, rows, cols);
-  //for(j=0;j<cols;++j){
-  //  for(i=0;i<rows;++i){
-  //    printf(" %3.2f ",orthoMatrix[j*rows + i]);
-  //  }
-  //  printf("\n");
-  //}
-}
+// QR to orthogonalize CSFs does not work
+//void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
+//  int i,j;
+//  //for(j=0;j<cols;++j){
+//  //  for(i=0;i<rows;++i){
+//  //    printf(" %3.2f ",overlapMatrix[j*rows + i]);
+//  //  }
+//  //  printf("\n");
+//  //}
+//  // Call the function ortho_qr from qp
+//  ortho_qr_csf(overlapMatrix, rows, orthoMatrix, rows, cols);
+//  //for(j=0;j<cols;++j){
+//  //  for(i=0;i<rows;++i){
+//  //    printf(" %3.2f ",orthoMatrix[j*rows + i]);
+//  //  }
+//  //  printf("\n");
+//  //}
+//}
 
 void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
 
diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f
index 21c19aaa..3fa437a5 100644
--- a/src/csf/sigma_vector.irp.f
+++ b/src/csf/sigma_vector.irp.f
@@ -112,7 +112,7 @@
   !  endif
   !enddo
   n_CSF = 0
-  print *," -9(((((((((((((( NSOMOMin=",NSOMOMin
+  !print *," -9(((((((((((((( NSOMOMin=",NSOMOMin
   ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1 
   if(ncfgprev.eq.-1)then
     ncfgprev=1
@@ -145,10 +145,10 @@
       endif
     endif
     n_CSF += ncfg*dimcsfpercfg
-    print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " ncfgprev=",ncfgprev, " senor=",cfg_seniority_index(i)
+    !print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " ncfgprev=",ncfgprev, " senor=",cfg_seniority_index(i)
     ncfgprev = cfg_seniority_index(i+2)
   end do
-  print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
+  !print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
 
 END_PROVIDER
 
@@ -240,7 +240,7 @@ end subroutine get_phase_qp_to_cfg
   ! initialization
   psi_coef_config = 0.d0
   DetToCSFTransformationMatrix(0,:,:) = 1.d0
-  do i = 2-iand(MS,1), NSOMOMax,2
+  do i = 2-iand(MS,1_8), NSOMOMax,2
     Isomo = IBSET(0_8, i) - 1_8
     ! rows = Ncsfs
     ! cols = Ndets
@@ -342,8 +342,8 @@ end subroutine get_phase_qp_to_cfg
   nsomomin = elec_alpha_num-elec_beta_num
   rowsmax = 0
   colsmax = 0
-  print *,"NSOMOMax = ",NSOMOMax
-  print *,"NSOMOMin = ",NSOMOMin
+  !print *,"NSOMOMax = ",NSOMOMax
+  !print *,"NSOMOMin = ",NSOMOMin
   !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
   ! Type
   ! 1. SOMO -> SOMO
@@ -525,7 +525,7 @@ end subroutine get_phase_qp_to_cfg
         end do
      end do
   end do
-  print *,"Rowsmax=",rowsmax," Colsmax=",colsmax
+  !print *,"Rowsmax=",rowsmax," Colsmax=",colsmax
   END_PROVIDER
 
   BEGIN_PROVIDER [ real*8, AIJpqContainer, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)]
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index 3463370f..16f67182 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1133,102 +1133,103 @@ subroutine ortho_svd(A,LDA,m,n)
 
 end
 
-subroutine ortho_qr_withB(A,LDA,B,m,n)
-  implicit none
-  BEGIN_DOC
-  ! Orthogonalization using Q.R factorization
-  !
-  ! A : Overlap Matrix
-  !
-  ! LDA : leftmost dimension of A
-  !
-  ! m : Number of rows of A
-  !
-  ! n : Number of columns of A
-  !
-  ! B : Output orthogonal basis
-  !
-  END_DOC
-  integer, intent(in)            :: m,n, LDA
-  double precision, intent(inout) :: A(LDA,n)
-  double precision, intent(inout) :: B(LDA,n)
+! QR to orthonormalize CSFs does not work :-(
+!subroutine ortho_qr_withB(A,LDA,B,m,n)
+!  implicit none
+!  BEGIN_DOC
+!  ! Orthogonalization using Q.R factorization
+!  !
+!  ! A : Overlap Matrix
+!  !
+!  ! LDA : leftmost dimension of A
+!  !
+!  ! m : Number of rows of A
+!  !
+!  ! n : Number of columns of A
+!  !
+!  ! B : Output orthogonal basis
+!  !
+!  END_DOC
+!  integer, intent(in)            :: m,n, LDA
+!  double precision, intent(inout) :: A(LDA,n)
+!  double precision, intent(inout) :: B(LDA,n)
+!
+!  integer                        :: LWORK, INFO
+!  integer, allocatable           :: jpvt(:)
+!  double precision, allocatable  :: TAU(:), WORK(:)
+!  double precision, allocatable  :: C(:,:)
+!  double precision :: norm
+!  integer :: i,j
+!
+!  allocate (TAU(min(m,n)), WORK(1))
+!  allocate (jpvt(n))
+!  !print *," In function ortho"
+!  B = A
+!
+!  jpvt(1:n)=1
+!
+!  LWORK=-1
+!  call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
+!
+!  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+!  LWORK=max(n,int(WORK(1)))
+!
+!  deallocate(WORK)
+!  allocate(WORK(LWORK))
+!  call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
+!  print *,A
+!  print *,jpvt
+!  deallocate(WORK,TAU)
+!  !stop
+!
+!  !LWORK=-1
+!  !call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
+!  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+!  !LWORK=max(n,int(WORK(1)))
+!
+!  !deallocate(WORK)
+!  !allocate(WORK(LWORK))
+!  !call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
+!
+!  !LWORK=-1
+!  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
+!  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+!  !LWORK=max(n,int(WORK(1)))
+!
+!  !deallocate(WORK)
+!  !allocate(WORK(LWORK))
+!  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
+!  !
+!  !allocate(C(LDA,n))
+!  !call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
+!  !norm = 0.0d0
+!  !B = 0.0d0
+!  !!print *,C
+!  !do i=1,m
+!  !  norm = 0.0d0
+!  !  do j=1,n
+!  !    norm = norm + C(j,i)*C(j,i)
+!  !  end do
+!  !  norm = 1.0d0/dsqrt(norm)
+!  !  do j=1,n
+!  !    B(j,i) = C(j,i)
+!  !  end do
+!  !end do
+!  !print *,B
+!
+!
+!  !deallocate(WORK,TAU)
+!end
 
-  integer                        :: LWORK, INFO
-  integer, allocatable           :: jpvt(:)
-  double precision, allocatable  :: TAU(:), WORK(:)
-  double precision, allocatable  :: C(:,:)
-  double precision :: norm
-  integer :: i,j
-
-  allocate (TAU(min(m,n)), WORK(1))
-  allocate (jpvt(n))
-  !print *," In function ortho"
-  B = A
-
-  jpvt(1:n)=1
-
-  LWORK=-1
-  call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
-
-  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
-  LWORK=max(n,int(WORK(1)))
-
-  deallocate(WORK)
-  allocate(WORK(LWORK))
-  call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
-  print *,A
-  print *,jpvt
-  deallocate(WORK,TAU)
-  !stop
-
-  !LWORK=-1
-  !call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
-  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
-  !LWORK=max(n,int(WORK(1)))
-
-  !deallocate(WORK)
-  !allocate(WORK(LWORK))
-  !call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
-
-  !LWORK=-1
-  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
-  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
-  !LWORK=max(n,int(WORK(1)))
-
-  !deallocate(WORK)
-  !allocate(WORK(LWORK))
-  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
-  !
-  !allocate(C(LDA,n))
-  !call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
-  !norm = 0.0d0
-  !B = 0.0d0
-  !!print *,C
-  !do i=1,m
-  !  norm = 0.0d0
-  !  do j=1,n
-  !    norm = norm + C(j,i)*C(j,i)
-  !  end do
-  !  norm = 1.0d0/dsqrt(norm)
-  !  do j=1,n
-  !    B(j,i) = C(j,i)
-  !  end do
-  !end do
-  !print *,B
-
-
-  !deallocate(WORK,TAU)
-end
-
-subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
-  use iso_c_binding
-  integer(c_int32_t), value      :: LDA
-  integer(c_int32_t), value      :: m
-  integer(c_int32_t), value      :: n
-  integer(c_int16_t)             :: A(LDA,n)
-  integer(c_int16_t)             :: B(LDA,n)
-  call ortho_qr_withB(A,LDA,B,m,n)
-end subroutine ortho_qr_csf
+!subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
+!  use iso_c_binding
+!  integer(c_int32_t), value      :: LDA
+!  integer(c_int32_t), value      :: m
+!  integer(c_int32_t), value      :: n
+!  integer(c_int16_t)             :: A(LDA,n)
+!  integer(c_int16_t)             :: B(LDA,n)
+!  call ortho_qr_withB(A,LDA,B,m,n)
+!end subroutine ortho_qr_csf
 
 subroutine ortho_qr(A,LDA,m,n)
   implicit none