subroutine print_spinfree_2rdm(rdm_defs, rdm, rdm_trace)
! Print all the RDM elements stored in rdm to a single file (for each
! state being sampled).
! The stem of the filenames is "spinfree_TwoRDM", and a final line
! required by MPQC to read spinfree 2-RDMs is also printed. This
! routine also assumes that the RDM element labels are already in
! spatial form, performing no transformation from spin to spatial
! form. This routine is therefore appropriate for printing spinfree
! 2-RDMs.
use Parallel_neci, only: MPIBarrier
use rdm_data, only: rdm_definitions_t
use sort_mod, only: sort
use util_mod, only: get_free_unit
implicit none
type(rdm_definitions_t), intent(in) :: rdm_defs
type(rdm_list_t), intent(in) :: rdm
real(dp), intent(in) :: rdm_trace(rdm%sign_length)
character(*), parameter :: this_routine = "print_spinfree_2rdm"
integer(int_rdm) :: pqrs
integer :: ielem, irdm, iunit, iproc, ierr
integer :: pq, rs, p, q, r, s ! spatial orbitals
integer(n_int) :: pq_, rs_
real(dp) :: rdm_sign(rdm%sign_length)
character(40) :: rdm_filename
associate(state_labels => rdm_defs%state_labels, repeat_label => rdm_defs%repeat_label)
do iproc = 0, nProcessors - 1
if (iproc == iProcIndex) then
! Loop over all RDMs beings sampled.
do irdm = 1, rdm_defs%nrdms
if (state_labels(1, irdm) == state_labels(2, irdm)) then
write(rdm_filename, '("spinfree_","'//trim(rdm_defs%output_file_prefix)//'",".",' &
//int_fmt(state_labels(1, irdm), 0)//')') irdm
else
write(rdm_filename, '("spinfree_","'//trim(rdm_defs%output_file_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)
end if
! Open the file to be written to.
iunit = get_free_unit()
! Let the first process clear the file, if it already exist.
if (iproc == 0) then
open(iunit, file=rdm_filename, status='replace')
else
open(iunit, file=rdm_filename, status='old', position='append')
end if
do ielem = 1, rdm%nelements
pqrs = rdm%elements(0, ielem)
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
! Normalise.
rdm_sign = rdm_sign / rdm_trace
if (tGUGA) then
! Obtain spatial orbital labels.
call extract_2_rdm_ind(pqrs, p, q, r, s, pq_, rs_)
pq = int(pq_)
rs = int(rs_)
else
! Obtain spin orbital labels.
call calc_separate_rdm_labels(pqrs, pq, rs, p, s, q, r)
end if
if (abs(rdm_sign(irdm)) > 1.e-12_dp) then
if (p >= q .and. pq >= rs .and. p >= r .and. p >= s) then
write(iunit, "(4I6, G25.17)") p, q, r, s, rdm_sign(irdm)
end if
end if
end do
close(iunit)
end do
end if
! Wait for the current processor to finish printing its RDM elements.
call MPIBarrier(ierr)
end do
end associate
end subroutine print_spinfree_2rdm