@brief Compute matrix inverse square root using SVD and removing
linear dependency
@detail This subroutine is used to obtain set of canonical orbitals
by diagonalization of the basis set overlap matrix
Q = S^{-1/2}, Q^T * S * Q = I
@param[in] s Overlap matrix, symmetric packed format
@param[out] q Matrix inverse square root, square matrix
@param[in] nbf Dimeension of matrices S and Q, basis set size
@param[out] qrnk Rank of matrix Q
@param[in] tol optional, tolerance to remove linear dependency,
default = 1.0e-8
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp), | intent(in) | :: | s(*) | |||
real(kind=dp), | intent(out) | :: | q(nbf,*) | |||
integer, | intent(in) | :: | nbf | |||
integer, | intent(out), | optional | :: | qrnk | ||
real(kind=dp), | intent(in), | optional | :: | tol |
subroutine matrix_invsqrt(s, q, nbf, qrnk, tol) use messages, only: show_message, with_abort use eigen, only: diag_symm_packed implicit none real(kind=dp), intent(in) :: s(*) real(kind=dp), intent(out) :: q(nbf,*) integer, intent(in) :: nbf real(kind=dp), optional, intent(in) :: tol integer, optional, intent(out) :: qrnk real(kind=dp), parameter :: deftol = 1.0d-08 real(kind=dp), allocatable :: tmp(:), eig(:) real(kind=dp) :: rtol integer :: nbf2, ok, i, j rtol = deftol if (present(tol)) rtol = tol nbf2 = nbf*(nbf+1)/2 allocate(tmp(nbf2), & eig(nbf), & stat=ok) if (ok/=0) call show_message('Cannot allocate memory', WITH_ABORT) tmp(:) = s(1:nbf2) ! Compute SVD call diag_symm_packed(1, nbf, nbf, nbf, tmp, eig, q, ok) ! Compute Q = S^{-1/2}, eliminating eigenvectors corresponding ! to small eigenvalues j = 0 do i = 1, nbf if (eig(i) >= rtol) then j = j+1 q(:,j) = q(:,i) / sqrt(eig(i)) end if end do q(:,j+1:nbf) = 0 if (present(qrnk)) qrnk = j end subroutine matrix_invsqrt