subroutine write_evales_and_transform_mat(rdm, irdm, SumDiag)
use LoggingData, only: tNoNOTransform
use rdm_data, only: tOpenShell, one_rdm_t
use SystemData, only: nbasis, nel, BRR
use UMatCache, only: gtID
use util_mod, only: get_free_unit, int_fmt
type(one_rdm_t), intent(in) :: rdm
integer, intent(in) :: irdm
real(dp), intent(in) :: SumDiag
integer :: i, j, Evalues_unit, NatOrbs_unit, jSpat, jInd, NO_Number
integer :: i_no, i_normal
real(dp) :: Corr_Entropy, Norm_Evalues, SumN_NO_Occ
logical :: tNegEvalue, tWrittenEvalue
character(20) :: evals_filename, transform_filename
if (tOpenShell) then
Norm_Evalues = SumDiag / real(NEl, dp)
else
Norm_Evalues = 2.0_dp * (SumDiag / real(NEl, dp))
end if
! Write out normalised Evalues to file and calculate the correlation
! entropy.
Corr_Entropy = 0.0_dp
Evalues_unit = get_free_unit()
write(evals_filename, '("NO_OCC_NUMEBRS.",'//int_fmt(irdm, 0)//')') irdm
open(Evalues_unit, file=trim(evals_filename), status='unknown')
write(Evalues_unit, '(A)') '# NOs (natural orbitals) ordered by occupation number'
write(Evalues_unit, '(A)') '# MOs (HF orbitals) ordered by energy'
write(Evalues_unit, '(A1,A5,A30,A20,A30)') '#', 'NO', 'NO OCCUPATION NUMBER', 'MO', 'MO OCCUPATION NUMBER'
tNegEvalue = .false.
SumN_NO_Occ = 0.0_dp
NO_Number = 1
do i = 1, NoOrbs
if (tOpenShell) then
write(Evalues_unit, '(I6,G35.17,I15,G35.17)') i, rdm%evalues(i) / Norm_Evalues, &
BRR(i), rdm%Rho_ii(i)
if (rdm%evalues(i) > 0.0_dp) then
Corr_Entropy = Corr_Entropy - (abs(rdm%evalues(i) / Norm_Evalues) &
* LOG(abs(rdm%evalues(i) / Norm_Evalues)))
else
tNegEvalue = .true.
end if
if (i <= NEl) SumN_NO_Occ = SumN_NO_Occ + (rdm%evalues(i) / Norm_Evalues)
else
write(Evalues_unit, '(I6,G35.17,I15,G35.17)') (2 * i) - 1, rdm%evalues(i) / Norm_Evalues, &
BRR((2 * i) - 1), rdm%Rho_ii(i) / 2.0_dp
if (rdm%evalues(i) > 0.0_dp) then
Corr_Entropy = Corr_Entropy - (2.0_dp * (abs(rdm%evalues(i) / Norm_Evalues) &
* LOG(abs(rdm%evalues(i) / Norm_Evalues))))
else
tNegEvalue = .true.
end if
write(Evalues_unit, '(I6,G35.17,I15,G35.17)') 2 * i, rdm%evalues(i) / Norm_Evalues, &
BRR(2 * i), rdm%Rho_ii(i) / 2.0_dp
if (i <= (NEl / 2)) SumN_NO_Occ = SumN_NO_Occ + (2.0_dp * (rdm%evalues(i) / Norm_Evalues))
end if
end do
close(Evalues_unit)
write(stdout, '(1X,A45,F30.20)') 'SUM OF THE N LARGEST NO OCCUPATION NUMBERS: ', SumN_NO_Occ
write(stdout, '(1X,A20,F30.20)') 'CORRELATION ENTROPY', Corr_Entropy
write(stdout, '(1X,A33,F30.20)') 'CORRELATION ENTROPY PER ELECTRON', Corr_Entropy / real(NEl, dp)
if (tNegEvalue) write(stdout, '(1X,"WARNING: Negative NO occupation numbers found.")')
! Write out the evectors to file.
! This is the matrix that transforms the molecular orbitals into the
! natural orbitals. rdm%evalues(i) corresponds to Evector NatOrbsMat(1:nBasis,i)
! We just want the rdm%evalues in the same order as above, but the
! 1:nBasis part (corresponding to the molecular orbitals), needs to
! refer to the actual orbital labels. Want these orbitals to preferably
! be in order, run through the orbital, need the position to find the
! corresponding NatOrbs element, use rdm%sym_list_inv_no.
associate(arr_ind => rdm%sym_list_inv_no)
if (.not. tNoNOTransform) then
NatOrbs_unit = get_free_unit()
write(transform_filename, '("NO_TRANSFORM.",'//int_fmt(irdm, 0)//')') irdm
open(NatOrbs_unit, file=trim(transform_filename), status='unknown')
write(NatOrbs_unit, '(2A6,2A30)') '# MO', 'NO', 'Transform Coeff', 'NO OCC NUMBER'
! write out in terms of spin orbitals, all alpha then all beta.
NO_Number = 1
do i_normal = 1, NoOrbs
i_no = arr_ind(i_normal)
tWrittenEvalue = .false.
do j = 1, nBasis
! Here i corresponds to the natural orbital, and j to the
! molecular orbital. i is actually the spin orbital in
! this case.
if (tOpenShell) then
jInd = j
else
if (mod(j, 2) /= 0) then
jInd = gtID(j)
else
cycle
end if
end if
if (tWrittenEvalue) then
if (.not. near_zero(rdm%matrix(arr_ind(jInd), i_no))) &
write(NatOrbs_unit, '(2I6,G35.17)') j, NO_Number, rdm%matrix(arr_ind(jInd), i_no)
else
if (.not. near_zero(rdm%matrix(arr_ind(jInd), i_no))) then
write(NatOrbs_unit, '(2I6,2G35.17)') j, NO_Number, rdm%matrix(arr_ind(jInd), i_no), &
rdm%evalues(i_no) / Norm_Evalues
tWrittenEvalue = .true.
end if
end if
end do
NO_Number = NO_Number + 1
if (.not. tOpenShell) then
tWrittenEvalue = .false.
do j = 2, nBasis, 2
! Here i corresponds to the natural orbital, and j to
! the molecular orbital. i is actually the spin orbital
! in this case.
jSpat = gtID(j)
if (tWrittenEvalue) then
if (.not. near_zero(rdm%matrix(arr_ind(jSpat), i_no))) &
write(NatOrbs_unit, '(2I6,G35.17)') j, NO_Number, &
rdm%matrix(arr_ind(jSpat), i_no)
else
if (.not. near_zero(rdm%matrix(arr_ind(jSpat), i_no))) then
write(NatOrbs_unit, '(2I6,2G35.17)') j, NO_Number, rdm%matrix(arr_ind(jSpat), i_no), &
rdm%evalues(i_no) / Norm_Evalues
tWrittenEvalue = .true.
end if
end if
end do
NO_Number = NO_Number + 1
end if
end do
close(NatOrbs_unit)
end if
end associate
end subroutine write_evales_and_transform_mat