anchor.diff:        ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
         ! -------------------------------------------
   
  -      !$OMP PARALLEL PRIVATE(k,l,iter2) SHARED(h,U,W,sze,iter)
  -      do l=1,N_st_diag
  -        !$OMP DO 
  -        do k=1,N_st_diag
  -          do iter2=1,iter-1
  -            h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze)
  -            h(k,iter,l,iter2) = h(k,iter2,l,iter)
  -          enddo
  -        enddo
  -        !$OMP END DO
  -      enddo
  -      do l=1,N_st_diag
  -        !$OMP DO 
  -        do k=1,l
  -          h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze)
  -          h(l,iter,k,iter) = h(k,iter,l,iter)
  -        enddo
  -        !$OMP END DO 
  -      enddo
  -      !$OMP END PARALLEL 
  -
  -      !DEBUG H MATRIX
  -      !do i=1,iter
  -      !  print '(10(x,F16.10))',  h(1,i,1,1:i)
  -      !enddo
  -      !print *,  ''
  -      !END
  -      
  +      call dgemm('T','N', N_st_diag*iter, N_st_diag, sze,               &
  +          1.d0, U, size(U,1), W(1,1,iter), size(W,1),   &
  +          0.d0, h(1,1,1,iter), size(h,1)*size(h,2))
  +
         ! Diagonalize h
         ! -------------
         call lapack_diag(lambda,y,h,N_st_diag*davidson_sze_max,N_st_diag*iter)
  
positive.code:       ! Gram-Schmidt
        ! ------------
        
        double precision               :: c(N_st_diag)
        do k=1,N_st_diag
          do iter2=1,iter
            call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter2),size(U,1),  &
                U(1,k,iter+1),1,0.d0,c,1)
            call dgemv('N',sze,N_st_diag,-1.d0,U(1,1,iter2),size(U,1), &
                c,1,1.d0,U(1,k,iter+1),1)
          enddo
          call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter+1),size(U,1),   &
              U(1,k,iter+1),1,0.d0,c,1)
          call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1),        &
              c,1,1.d0,U(1,k,iter+1),1)
          call normalize( U(1,k,iter+1), sze )
        enddo
  
      
        !DEBUG : CHECK OVERLAP
        !print *,  '==='
        !do k=1,iter+1
        !  do l=1,k
        !  c = u_dot_v(U(1,1,k),U(1,1,l),sze)
        !  print *,  k,l, c
        !  enddo
        !enddo
        !print *,  '==='
        !pause
        !END DEBUG
  
  
      enddo
  
  
negative.code:       enddo
      enddo
      
      ! Compute overlap matrix Ut.U to orthonormalize U
      ! -----------------------------------------------
  
      call DGEMM('T','N', N_st_diag, N_st_diag, sze,                   &
          1.d0, U(1,1,1), size(U,1), U(1,1,1), size(U,1),              &
          0.d0, overlap, size(overlap,1))
      
      call ortho_lowdin(overlap,size(overlap,1),N_st_diag,U,size(U,1),sze)
  
      do iter=1,davidson_sze_max-1
        
        ! Compute |W_k> = \sum_i |i><i|H|u_k>
        ! -----------------------------------------
        
        call H_u_0_nstates(W(1,1,iter),U(1,1,iter),H_jj,sze,dets_in,Nint,N_st_diag,sze_8)
        
  