subroutine print_rdms_with_spin(rdm_defs, nrdms_to_print, rdm, rdm_trace, open_shell)
! Print the RDM stored in rdm to files, normalised by rdm_trace.
! This routine will print out *all* the spin cobminations separately,
! including both aaaa and bbbb arrays, and all other combinations.
! The files are called 'TwoRDM_aaaa', 'TwoRDM_abab', 'TwoRDM_abba', etc...
use Parallel_neci, only: MPIBarrier
use rdm_data, only: rdm_definitions_t
use sort_mod, only: sort
use UMatCache, only: spatial
use util_mod, only: get_free_unit
type(rdm_definitions_t), intent(in) :: rdm_defs
integer, intent(in) :: nrdms_to_print
type(rdm_list_t), intent(inout) :: rdm
real(dp), intent(in) :: rdm_trace(rdm%sign_length)
logical, intent(in) :: open_shell
integer(int_rdm) :: ijkl
integer :: ij, kl, i, j, k, l ! spin orbitals
integer :: p, q, r, s ! spatial orbitals
integer :: ielem, irdm, ierr, iproc, write_unit
integer :: iunit_aaaa, iunit_abab, iunit_abba
integer :: iunit_bbbb, iunit_baba, iunit_baab
real(dp) :: rdm_sign(rdm%sign_length)
character(12) :: sgn_len, suffix
character(len=*), parameter :: t_r = 'print_rdms_with_spin'
associate(state_labels => rdm_defs%state_labels, repeat_label => rdm_defs%repeat_label)
! Store rdm%sign_length as a string, for the formatting string.
write(sgn_len, '(i3)') rdm%sign_length
call sort(rdm%elements(:, 1:rdm%nelements))
do iproc = 0, nProcessors - 1
do irdm = 1, nrdms_to_print
if (state_labels(1, irdm) == state_labels(2, irdm)) then
write(suffix, '('//int_fmt(state_labels(1, irdm), 0)//')') irdm
else
write(suffix, '('//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)
end if
if (iproc == iProcIndex) then
! Open all the files to be written to:
! Let the first processor clear all the files to start with.
if (iproc == 0) then
iunit_aaaa = get_free_unit()
open(iunit_aaaa, file=trim(rdm_defs%output_file_prefix)// &
'_aaaa.'//trim(suffix), status='replace')
iunit_abab = get_free_unit()
open(iunit_abab, file=trim(rdm_defs%output_file_prefix)// &
'_abab.'//trim(suffix), status='replace')
iunit_abba = get_free_unit()
open(iunit_abba, file=trim(rdm_defs%output_file_prefix)// &
'_abba.'//trim(suffix), status='replace')
if (open_shell) then
iunit_bbbb = get_free_unit()
open(iunit_bbbb, file=trim(rdm_defs%output_file_prefix)// &
'_bbbb.'//trim(suffix), status='replace')
iunit_baba = get_free_unit()
open(iunit_baba, file=trim(rdm_defs%output_file_prefix)// &
'_baba.'//trim(suffix), status='replace')
iunit_baab = get_free_unit()
open(iunit_baab, file=trim(rdm_defs%output_file_prefix)// &
'_baab.'//trim(suffix), status='replace')
end if
else
iunit_aaaa = get_free_unit()
open(iunit_aaaa, file=trim(rdm_defs%output_file_prefix)// &
'_aaaa.'//trim(suffix), status='old', position='append')
iunit_abab = get_free_unit()
open(iunit_abab, file=trim(rdm_defs%output_file_prefix)// &
'_abab.'//trim(suffix), status='old', position='append')
iunit_abba = get_free_unit()
open(iunit_abba, file=trim(rdm_defs%output_file_prefix)// &
'_abba.'//trim(suffix), status='old', position='append')
if (open_shell) then
iunit_bbbb = get_free_unit()
open(iunit_bbbb, file=trim(rdm_defs%output_file_prefix)// &
'_bbbb.'//trim(suffix), status='old', position='append')
iunit_baba = get_free_unit()
open(iunit_baba, file=trim(rdm_defs%output_file_prefix)// &
'_baba.'//trim(suffix), status='old', position='append')
iunit_baab = get_free_unit()
open(iunit_baab, file=trim(rdm_defs%output_file_prefix)// &
'_baab.'//trim(suffix), status='old', position='append')
end if
end if
do ielem = 1, rdm%nelements
ijkl = rdm%elements(0, ielem)
! Obtain spin orbital labels.
call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
! Normalise.
rdm_sign = rdm_sign / rdm_trace
p = spatial(i); q = spatial(j);
r = spatial(k); s = spatial(l);
if ((.not. open_shell) .and. is_beta(i)) then
call stop_all(t_r, "This is a closed shell system but we have an open shell type RDM element.&
& An error must have occured.")
end if
! Find out what the spin labels are, and print the RDM
! element to the appropriate file.
if (is_alpha(i) .and. is_alpha(j) .and. is_alpha(k) .and. is_alpha(l)) then
write_unit = iunit_aaaa
else if (is_alpha(i) .and. is_beta(j) .and. is_alpha(k) .and. is_beta(l)) then
write_unit = iunit_abab
else if (is_alpha(i) .and. is_beta(j) .and. is_beta(k) .and. is_alpha(l)) then
write_unit = iunit_abba
else if (is_beta(i) .and. is_beta(j) .and. is_beta(k) .and. is_beta(l)) then
write_unit = iunit_bbbb
else if (is_beta(i) .and. is_alpha(j) .and. is_beta(k) .and. is_alpha(l)) then
write_unit = iunit_baba
else if (is_beta(i) .and. is_alpha(j) .and. is_alpha(k) .and. is_beta(l)) then
write_unit = iunit_baab
end if
if (abs(rdm_sign(irdm)) > 1.e-12_dp) then
write(write_unit, '(4i6,'//trim(sgn_len)//'g25.17)') p, q, r, s, rdm_sign(irdm)
end if
end do
close(iunit_aaaa); close (iunit_abab); close (iunit_abba);
if (open_shell) then
close(iunit_bbbb); close (iunit_baba); close (iunit_baab);
end if
end if
end do
! Wait for the current processor to finish printing its RDM elements.
call MPIBarrier(ierr)
end do
end associate
end subroutine print_rdms_with_spin