anchor.diff:    integer                        :: iter2
     double precision, allocatable  :: W(:,:,:),  U(:,:,:), R(:,:)
     double precision, allocatable  :: y(:,:,:,:), h(:,:,:,:), lambda(:)
  +  double precision, allocatable  :: c(:), H_small(:,:)
     double precision               :: diag_h_mat_elem
     double precision, allocatable  :: residual_norm(:)
     character*(16384)              :: write_buffer
  
positive.code:       ! Express eigenvectors of h in the determinant basis
        ! --------------------------------------------------
        
       !$OMP PARALLEL DEFAULT(NONE)                                     &
          !$OMP PRIVATE(k,i,l,iter2) &
          !$OMP SHARED(U,W,R,y,iter,lambda,N_st_diag,sze,to_print, &
          !$OMP residual_norm,nuclear_repulsion,N_st)
        do k=1,N_st_diag
          !$OMP DO
          do i=1,sze
            U(i,k,iter+1) = 0.d0
            W(i,k,iter+1) = 0.d0
          enddo
          !$OMP END DO
           do iter2=1,iter
            do l=1,N_st_diag
              !$OMP DO
              do i=1,sze
                U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
                W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
              enddo
              !$OMP END DO NOWAIT
            enddo
          enddo
        
          ! Compute residual vector
          ! -----------------------
        
          !$OMP DO
          do i=1,sze
            R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1)
          enddo
          !$OMP END DO
          !$OMP SINGLE
          if (k <= N_st) then
            residual_norm(k) = u_dot_u(R(1,k),sze)
            to_print(1,k) = lambda(k) + nuclear_repulsion
            to_print(2,k) = residual_norm(k)
          endif
          !$OMP END SINGLE
        enddo
        !$OMP END PARALLEL
        
        write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))')  iter, to_print(:,1:N_st)
        call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
  
negative.code:       y(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max),      &
        residual_norm(N_st_diag),                                      &
        overlap(N_st_diag,N_st_diag),                                  &
        c(N_st_diag*davidson_sze_max),                                 &
        H_small(N_st_diag,N_st_diag),                                  &
        lambda(N_st_diag*davidson_sze_max))
    
    ASSERT (N_st > 0)
  