subroutine write_1rdm(rdm_defs, one_rdm, irdm, norm_1rdm, tNormalise, tInitsRDM)
! This routine writes out the OneRDM. If tNormalise is true, we are
! printing the normalised, hermitian matrix. Otherwise, norm_1rdm is
! ignored and we print both 1-RDM(i,j) and 1-RDM(j,i) (in binary)
! for the OneRDM_POPS file to be read in in a restart calculation.
use rdm_data, only: tOpenShell, rdm_definitions_t
use RotateOrbsData, only: SymLabelListInv_rot
use SystemData, only: nbasis
use UMatCache, only: gtID
use util_mod, only: get_free_unit, int_fmt
type(rdm_definitions_t), intent(in) :: rdm_defs
real(dp), intent(in) :: one_rdm(:, :)
integer, intent(in) :: irdm
real(dp), intent(in) :: norm_1rdm
logical, intent(in) :: tNormalise
logical, intent(in), optional :: tInitsRDM
integer :: i, j, iSpat, jSpat, one_rdm_unit, one_rdm_unit_spinfree
logical :: is_transition_rdm
character(20) :: filename
character(20) :: filename_prefix
character(*), parameter :: default_prefix = "OneRDM."
if (present(tInitsRDM)) then
if (tInitsRDM) then
filename_prefix = "InitsOneRDM."
else
filename_prefix = default_prefix
end if
else
filename_prefix = default_prefix
end if
associate (state_labels => rdm_defs%state_labels, &
repeat_label => rdm_defs%repeat_label, &
ind => SymLabelListInv_rot)
is_transition_rdm = state_labels(1, irdm) /= state_labels(2, irdm)
if (tWriteSpinFreeRDM .or. tGUGA) then
one_rdm_unit_spinfree = get_free_unit()
open(one_rdm_unit_spinfree, file="spinfree-1-RDM")
end if
if (.not. tGUGA) then
if (tNormalise) then
! Haven't got the capabilities to produce multiple 1-RDMs yet.
write(stdout, '(1X,"Writing out the *normalised* 1 electron density matrix to file")')
call neci_flush(stdout)
one_rdm_unit = get_free_unit()
if (is_transition_rdm) then
write(filename, '("'//trim(filename_prefix)//'",' &
//int_fmt(state_labels(1, irdm), 0)//',"_",' &
//int_fmt(state_labels(2, irdm), 0)//',".",i1)') &
state_labels(1, irdm), state_labels(2, irdm), repeat_label(irdm)
else
write(filename, '("'//trim(filename_prefix)//'",' &
//int_fmt(state_labels(1, irdm), 0)//')') irdm
end if
open(one_rdm_unit, file=trim(filename), status='unknown')
else
! Only every write out 1 of these at the moment.
write(stdout, '(1X,"Writing out the *unnormalised* 1 electron density matrix to file for reading in")')
call neci_flush(stdout)
one_rdm_unit = get_free_unit()
if (is_transition_rdm) then
write(filename, '("OneRDM_POPS.",'//int_fmt(state_labels(1, irdm), 0)//',"_",' &
//int_fmt(state_labels(2, irdm), 0)//',".",i1)') &
state_labels(1, irdm), state_labels(2, irdm), repeat_label(irdm)
else
write(filename, '("OneRDM_POPS.",'//int_fmt(state_labels(1, irdm), 0)//')') irdm
end if
open(one_rdm_unit, file=trim(filename), status='unknown', form='unformatted')
end if
end if
if (tGUGA) then
do i = 1, nSpatorbs
do j = 1, nSpatorbs
if (abs(one_rdm(ind(i), ind(j))) > EPS) then
if (tNormalise) then
if (i <= j) then
write(one_rdm_unit_spinfree, "(2i6, g25.17)") i, j, &
(one_rdm(ind(i), ind(j)) * norm_1rdm)
end if
else
write(one_rdm_unit_spinfree) i, j, one_rdm(ind(i), ind(j))
end if
end if
end do
end do
else
! Currently always printing 1-RDM in spin orbitals.
do i = 1, nbasis
do j = 1, nbasis
if (tOpenShell) then
if (abs(one_rdm(ind(i), ind(j))) > 1.0e-12_dp) then
if (tNormalise .and. (i <= j .or. is_transition_rdm)) then
write(one_rdm_unit, "(2I6,G25.17)") i, j, one_rdm(ind(i), ind(j)) * norm_1rdm
else if (.not. tNormalise) then
! For the pops, we haven't made the 1-RDM hermitian yet,
! so print both the 1-RDM(i,j) and 1-RDM(j,i) elements.
! This is written in binary.
write(one_rdm_unit) i, j, one_rdm(ind(i), ind(j))
end if
end if
else
iSpat = gtID(i)
jSpat = gtID(j)
if (abs(one_rdm(ind(iSpat), ind(jSpat))) > 1.0e-12_dp) then
if (tNormalise .and. (i <= j .or. is_transition_rdm)) then
if ((mod(i, 2) == 0 .and. mod(j, 2) == 0) .or. &
(mod(i, 2) /= 0 .and. mod(j, 2) /= 0)) then
write(one_rdm_unit, "(2I6,G25.17)") i, j, &
(one_rdm(ind(iSpat), ind(jSpat)) * norm_1rdm) / 2.0_dp
end if
else if (.not. tNormalise) then
! The popsfile can be printed in spatial orbitals.
if (mod(i, 2) == 0 .and. mod(j, 2) == 0) then
write(one_rdm_unit) iSpat, jSpat, one_rdm(ind(iSpat), ind(jSpat))
end if
end if
end if
end if
end do
end do
! if we want spin-free RDMs also print out the spin-free 1-RDM
if (tWriteSpinFreeRDM .and. tNormalise) then
do i = 1, nbasis / 2
do j = 1, nbasis / 2
if (abs(one_rdm(ind(i), ind(j))) > EPS) then
if (tNormalise) then
if (i <= j .or. is_transition_rdm) then
write(one_rdm_unit_spinfree, "(2I6,G25.17)") i, j, &
one_rdm(ind(i), ind(j)) * norm_1rdm
end if
else
write(one_rdm_unit_spinfree) i, j, one_rdm(ind(i), ind(j))
end if
end if
end do
end do
end if
end if ! tGUGA
if (.not. tGUGA) close(one_rdm_unit)
if (tWriteSpinFreeRDM .or. tGUGA) close(one_rdm_unit_spinfree)
end associate
end subroutine write_1rdm