diff --git a/src/csf/cfgCI_interface.f90 b/src/csf/cfgCI_interface.f90
index b701f0ec..73bd600d 100644
--- a/src/csf/cfgCI_interface.f90
+++ b/src/csf/cfgCI_interface.f90
@@ -46,6 +46,24 @@ module cfunctions
            real   (kind=C_DOUBLE ),intent(out) :: csftodetmatrix(rowsmax,colsmax)
          end subroutine getCSFtoDETTransformationMatrix
       end interface
+       interface
+         subroutine gramSchmidt(A, m, n, B) bind(C, name='gramSchmidt')
+           import C_INT32_T, C_INT64_T, C_DOUBLE
+           integer(kind=C_INT32_T),value,intent(in) :: m
+           integer(kind=C_INT32_T),value,intent(in) :: n
+           real (kind=C_DOUBLE ),intent(in) :: A(m,n)
+           real (kind=C_DOUBLE ),intent(out) :: B(m,n)
+         end subroutine gramSchmidt
+      end interface
+       interface
+         subroutine gramSchmidt_qp(A, m, n, B) bind(C, name='gramSchmidt_qp')
+           import C_INT32_T, C_INT64_T, C_DOUBLE
+           integer(kind=C_INT32_T),value,intent(in) :: m
+           integer(kind=C_INT32_T),value,intent(in) :: n
+           real (kind=C_DOUBLE ),intent(in) :: A(m,n)
+           real (kind=C_DOUBLE ),intent(out) :: B(m,n)
+         end subroutine gramSchmidt_qp
+      end interface
 end module cfunctions
 
 subroutine f_dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) &
diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c
index 6094f5bb..64593159 100644
--- a/src/csf/cfgCI_utils.c
+++ b/src/csf/cfgCI_utils.c
@@ -252,6 +252,26 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
     buildTreeDriver(bftree, *NSOMO, MS, NBF);
 }
 
+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");
+  //}
+}
+
 void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
 
     // vector
diff --git a/src/csf/tree_utils.h b/src/csf/tree_utils.h
index 5a788149..42fce59c 100644
--- a/src/csf/tree_utils.h
+++ b/src/csf/tree_utils.h
@@ -47,6 +47,7 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
 void getSetBits(int64_t n, int *nsetbits);
 void getOverlapMatrix(int64_t Isomo, int64_t MS, double **overlapMatrixptr, int *rows, int *cols, int *NSOMOout);
 void getOverlapMatrix_withDet(double *bftodetmatrixI, int rowsbftodetI, int colsbftodetI, int64_t Isomo, int64_t MS, double **overlapMatrixI, int *rowsI, int *colsI, int *NSOMO);
+void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
 void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
 
 
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index 405d2d20..3463370f 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1133,6 +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)
+
+  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(A,LDA,m,n)
   implicit none
   BEGIN_DOC