@brief Compute "energy weighted density matrix"
Note
This quantity is actually the Lagrangian matrix,
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=dp) | :: | eps(:) | ||||
integer | :: | nbf | ||||
type(information), | intent(inout) | :: | infos |
subroutine eijden(eps, nbf, infos) use oqp_tagarray_driver use mathlib, only: orthogonal_transform_sym use messages, only: show_message, with_abort implicit none character(len=*), parameter :: subroutine_name = "eijden" type(information), intent(inout) :: infos integer :: nbf integer :: i, j, ij, ne, ok real(kind=dp) :: eps(:) real(kind=dp), allocatable :: c(:,:), scr(:),tempd(:) ! tagarray real(kind=dp), contiguous, pointer :: & mo_energy_a(:), mo_a(:,:), & fock_a(:), fock_b(:), dmat_a(:), dmat_b(:) character(len=*), parameter :: tags_alpha(2) = (/ character(len=80) :: & OQP_E_MO_A, OQP_VEC_MO_A /) character(len=*), parameter :: tags_beta(4) = (/ character(len=80) :: & OQP_FOCK_A, OQP_DM_A, OQP_FOCK_B, OQP_DM_B /) if (infos%control%scftype>1) then allocate(c(nbf,nbf), scr(8*nbf), tempd(nbf*(nbf+1)/2), stat=ok) end if ne = infos%mol_prop%nelec/2 select case (infos%control%scftype) ! RHF case case (1) call data_has_tags(infos%dat, tags_alpha, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_energy_a) call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a) ij = 0 do i = 1, nbf do j = 1, i ij = ij+1 eps(ij) = -2*sum(mo_energy_a(1:ne)& *mo_a(i,1:ne)& *mo_a(j,1:ne)) end do end do ! U/ROHF case case (2:) call data_has_tags(infos%dat, tags_beta, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_FOCK_A, fock_a) call tagarray_get_data(infos%dat, OQP_FOCK_B, fock_b) call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a) call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b) ! Alpha part call unpack_matrix(dmat_a,c,nbf,'U') call orthogonal_transform_sym(nbf, nbf, fock_a, c, nbf, tempd) ! Beta part call unpack_matrix(dmat_b,c,nbf,'U') call orthogonal_transform_sym(nbf, nbf, fock_b, c, nbf, eps) eps = -eps - tempd ! Half the diagonal ij = 0 do i = 1, nbf ij = ij+i eps(ij) = 0.5d0*eps(ij) end do end select end subroutine eijden