@brief Solve \f$ F C = \eps S C \f$
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp) | :: | fock_m(*) | ||||
real(kind=dp) | :: | orbitals(:,:) | ||||
real(kind=dp) | :: | orbital_e(*) | ||||
real(kind=dp) | :: | q_m(:,:) |
subroutine get_ab_initio_orbital(fock_m, orbitals, orbital_e, q_m) use messages, only: show_message, WITH_ABORT use mathlib, only: orthogonal_transform use eigen, only: diag_symm_full use mathlib, only: unpack_matrix, pack_matrix implicit none ! real(kind=dp) :: orbital_e(*) real(kind=dp) :: fock_m(*) real(kind=dp) :: orbitals(:,:), q_m(:,:) real(kind=dp), allocatable :: wrk(:,:), f2(:,:), tfock(:,:) integer :: info, ok, nbf nbf = ubound(orbitals, 1) allocate(wrk(nbf,nbf), f2(nbf,nbf), tfock(nbf, nbf), stat=ok) if (ok/=0) call show_message('Cannot allocate memory', WITH_ABORT) call unpack_matrix(fock_m, f2) ! F = S^{-1/2} * F * S^{-1/2} ! can use orthogonal_transform subroutine, because S^{-1/2} is symmetric call orthogonal_transform('n', nbf, q_m, f2, tfock, wrk) deallocate(f2, wrk) ! Get orbitals C' = S^{1/2} C call diag_symm_full(1, nbf, tfock, nbf, orbital_e, info) ! Recover orbitals in AO basis ! C = S^{-1/2} * C' call dgemm('n', 'n', nbf, nbf, nbf, & 1.0_dp, q_m, nbf, & tfock, nbf, & 0.0_dp, orbitals, nbf) end subroutine get_ab_initio_orbital