rpanewb Subroutine

  subroutine rpanewb(ndsr,bvec,q,novec,nvec,ick,tamm_dancoff)
    use precision, only: dp

    implicit none

    integer, intent(in) :: ndsr
    real(kind=dp), intent(out), dimension(:,:) :: bvec
    real(kind=dp), intent(inout), dimension(:,:) :: q
    integer, intent(inout) :: novec, nvec
    integer, intent(out) :: ick
    logical, intent(in) :: tamm_dancoff

    real(kind=dp) :: bq, fnorm
    integer :: istat, k, ms, ndsrt, mxvec
    real(kind=dp), parameter :: norm_threshold = 1.0D-09

    mxvec = ubound(bvec, 2)

!   Save nvec as novec
    novec = nvec
    ick = 0

!   Modified Gram-Schmidt orthonormalization
    ndsrt = ndsr*2
    if (tamm_dancoff) ndsrt = ndsr

    do k = 1, ndsrt
!     MGS: orthonormalize next vector w.r.t. all
!     previous vectors
      do istat = 1, nvec
        bq = dot_product(bvec(:,istat),q(:,k))
        q(:,k) = q(:,k) - bq*bvec(:,istat)
      end do

      fnorm = norm2(q(:,k))

!     Possible linear dependency, skip this vector
      if (fnorm<norm_threshold) cycle

      if (nvec==mxvec) then
!       Error termination, no space left for new vectors
        ick = 2
        return
      end if

!     Append new b vector
      nvec = nvec+1
      bvec(:,nvec) = q(:,k)/fnorm
    end do

!   Error termination, no vectors added
    ms = nvec-novec
    if (ms==0) ick = 3
  end subroutine rpanewb