subroutine build_pfon_density(pdmat, mo_a, mo_b, occ_a, occ_b, scf_type, nbf, nelec_a, nelec_b)
use precision, only: dp
use mathlib, only: pack_matrix
implicit none
real(kind=dp), intent(inout) :: pdmat(:,:)
real(kind=dp), intent(in) :: mo_a(:,:), mo_b(:,:)
real(kind=dp), intent(in) :: occ_a(:), occ_b(:)
integer, intent(in) :: nbf, scf_type, nelec_a, nelec_b
real(kind=dp), allocatable :: dtmp(:,:)
integer :: i, mu, nu
integer :: n_double, n_single
real(kind=dp) :: occ_factor
allocate(dtmp(nbf, nbf), source=0.0_dp)
pdmat(:,:) = 0.0_dp
select case(scf_type)
case(1) ! RHF
do i = 1, nbf
if (occ_a(i) > 1.0e-14_dp) then
do mu = 1, nbf
do nu = 1, nbf
dtmp(mu,nu) = dtmp(mu,nu) + occ_a(i) * mo_a(mu,i)*mo_a(nu,i)
end do
end do
end if
end do
call pack_matrix(dtmp, pdmat(:,1))
case(2) ! UHF
do i = 1, nbf
if (occ_a(i) > 1.0e-14_dp) then
do mu = 1, nbf
do nu = 1, nbf
dtmp(mu,nu) = dtmp(mu,nu) + occ_a(i) * mo_a(mu,i)*mo_a(nu,i)
end do
end do
end if
end do
call pack_matrix(dtmp, pdmat(:,1))
dtmp(:,:) = 0.0_dp
do i = 1, nbf
if (occ_b(i) > 1.0e-14_dp) then
do mu = 1, nbf
do nu = 1, nbf
dtmp(mu,nu) = dtmp(mu,nu) + occ_b(i) * mo_b(mu,i)*mo_b(nu,i)
end do
end do
end if
end do
call pack_matrix(dtmp, pdmat(:,2))
case(3) ! ROHF
n_double = nelec_b
n_single = nelec_a - nelec_b
dtmp(:,:) = 0.0_dp
do i = 1, nbf
if (occ_a(i) > 1.0e-14_dp) then
if (i <= n_double) then
occ_factor = occ_a(i)
else if (i <= n_double + n_single) then
occ_factor = 1.0_dp
else
occ_factor = occ_a(i) ! Virtual orbitals
endif
do mu = 1, nbf
do nu = 1, nbf
dtmp(mu,nu) = dtmp(mu,nu) + occ_factor * mo_a(mu,i)*mo_a(nu,i)
end do
end do
end if
end do
call pack_matrix(dtmp, pdmat(:,1))
dtmp(:,:) = 0.0_dp
do i = 1, nbf
if (occ_b(i) > 1.0e-14_dp) then
if (i <= n_double) then
occ_factor = occ_b(i)
else
occ_factor = 0.0_dp
endif
do mu = 1, nbf
do nu = 1, nbf
dtmp(mu,nu) = dtmp(mu,nu) + occ_factor * mo_a(mu,i)*mo_a(nu,i)
end do
end do
end if
end do
call pack_matrix(dtmp, pdmat(:,2))
end select
deallocate(dtmp)
end subroutine build_pfon_density