subroutine finalise_rdms(rdm_defs, one_rdms, two_rdms, rdm_recv, rdm_recv_2, en_pert, spawn, rdm_estimates, tInitsRDMs)
! Wrapper routine, called at the end of a simulation, which in turn
! calls all required finalisation routines.
use SystemData, only: called_as_lib
use LoggingData, only: tBrokenSymNOs, occ_numb_diff, RDMExcitLevel, tExplicitAllRDM
use LoggingData, only: tPrint1RDM, tDiagRDM, tDumpForcesInfo
use LoggingData, only: tDipoles, tWrite_normalised_RDMs
use Parallel_neci, only: iProcIndex, MPIBarrier, MPIBCast
use rdm_data, only: tRotatedNos, FinaliseRDMs_Time, tOpenShell, print_2rdm_est
use rdm_data, only: rdm_list_t, rdm_spawn_t, one_rdm_t, rdm_estimates_t
use rdm_data, only: rdm_definitions_t, en_pert_t
use rdm_estimators, only: calc_2rdm_estimates_wrapper, write_rdm_estimates
use rdm_nat_orbs, only: find_nat_orb_occ_numbers, BrokenSymNo
use rdm_hdf5, only: write_rdms_hdf5
use timing_neci, only: set_timer, halt_timer
type(rdm_definitions_t), intent(in) :: rdm_defs
type(one_rdm_t), intent(inout) :: one_rdms(:)
type(rdm_list_t), intent(in) :: two_rdms
type(rdm_list_t), intent(inout) :: rdm_recv, rdm_recv_2
type(en_pert_t), intent(in) :: en_pert
type(rdm_spawn_t), intent(inout) :: spawn
type(rdm_estimates_t), intent(inout) :: rdm_estimates
logical, intent(in) :: tInitsRDMs
integer :: irdm, ierr
real(dp) :: norm_1rdm(size(one_rdms)), SumN_Rho_ii(size(one_rdms))
character(255) :: RDMName
call set_timer(FinaliseRDMs_Time)
if (tInitsRDMS) then
RDMName = 'RDMs (Initiators)'
else if (tAdaptiveShift) then
RDMName = 'RDMs (Lagrangian)'
else
RDMName = 'RDMs'
end if
if (tExplicitAllRDM) then
write(stdout, '(/,"****","'//trim(RDMName)//'"," CALCULATED EXPLICITLY ****",1X,/)')
else
write(stdout, '(/,"****","'//trim(RDMName)//'"," CALCULATED STOCHASTICALLY ****",1X,/)')
end if
! Combine the 1- or 2-RDM from all processors, etc.
! Stuff using the 2-RDMs:
if (RDMExcitLevel /= 1) then
! We always want to calculate one final RDM energy, whether or not we're
! calculating the energy throughout the calculation.
! Unless of course, only the 1-RDM is being calculated.
! Calculate the RDM estmimates from the final few iterations,
! since it was last calculated. But only do this if they're
! actually to be printed.
call calc_2rdm_estimates_wrapper(rdm_defs, rdm_estimates, two_rdms, en_pert)
! Output the final 2-RDMs themselves, in all forms desired.
call output_2rdm_wrapper(rdm_defs, rdm_estimates, two_rdms, rdm_recv, rdm_recv_2, spawn)
! Calculate the 1-RDMs from the 2-RDMS, if required.
if (RDMExcitLevel == 3 .or. tDiagRDM .or. tPrint1RDM .or. tDumpForcesInfo .or. tDipoles) then
if (tGUGA) then
call calc_1rdms_from_spinfree_2rdms(one_rdms, two_rdms, rdm_estimates%norm)
else
call calc_1rdms_from_2rdms(rdm_defs, one_rdms, two_rdms, rdm_estimates%norm, tOpenShell)
end if
! The 1-RDM will have been constructed to be normalised already.
norm_1rdm = 1.0_dp
end if
end if
! Stuff using the 1-RDMs:
if (RDMExcitLevel == 1 .or. RDMExcitLevel == 3) then
! Output banner for start of 1-RDM section in the output.
write(stdout, '(1x,2("="),1x,"INFORMATION FOR FINAL 1-","'//trim(RDMName)//'",1x,57("="))')
if (RDMExcitLevel == 1) call finalise_1e_rdm(rdm_defs, one_rdms, norm_1rdm)
if (iProcIndex == 0) then
call calc_rho_ii_and_sum_n(one_rdms, norm_1rdm, SumN_Rho_ii)
do irdm = 1, size(one_rdms)
write(stdout, '(/,1x,"INFORMATION FOR 1-RDM",1x,'//int_fmt(irdm)//',":")') irdm
write(stdout, '(/,1X,"SUM OF 1-RDM(i,i) FOR THE N LOWEST ENERGY HF ORBITALS:",1X,F20.13)') SumN_Rho_ii(irdm)
if (RDMExcitLevel == 1 .or. tPrint1RDM) then
! Write out the final, normalised, hermitian OneRDM.
call write_1rdm(rdm_defs, one_rdms(irdm)%matrix, irdm, norm_1rdm(irdm), .true., tInitsRDMs)
end if
end do
end if
call MPIBarrier(ierr)
do irdm = 1, rdm_defs%nrdms_standard
! Call the routines from NatOrbs that diagonalise the one electron
! reduced density matrix.
tRotatedNOs = .false. ! Needed for BrokenSymNo routine
if (tDiagRDM) call find_nat_orb_occ_numbers(one_rdms(irdm), irdm)
! After all the NO calculations are finished we'd like to do another
! rotation to obtain symmetry-broken natural orbitals.
if (tBrokenSymNOs) then
call BrokenSymNO(one_rdms(irdm)%evalues, occ_numb_diff)
end if
end do
! Output banner for the end of the 1-RDM section.
write(stdout, '(/,1x,89("="),/)')
end if
if (t_print_hdf5_rdms .and. tWriteSpinFreeRDM) then
call create_spinfree_2rdm(two_rdms, rdm_defs%nrdms_standard, spawn, rdm_recv)
call calc_1rdms_from_2rdms(rdm_defs, one_rdms, two_rdms, rdm_estimates%norm, tOpenShell)
call write_rdms_hdf5(rdm_defs, rdm_recv, rdm_estimates%norm, one_rdms)
end if
! Write the final instantaneous 2-RDM estimates, and also the final
! report of the total 2-RDM estimates.
if (RDMExcitLevel /= 1 .and. iProcIndex == 0) then
call write_rdm_estimates(rdm_defs, rdm_estimates, .true., print_2rdm_est, &
tInitsRDMs)
end if
if (print_2rdm_est .and. called_as_lib) then
RDM_energy = rdm_estimates%energy_num(1) / rdm_estimates%norm(1)
call MPIBCast(RDM_energy)
write(stdout, *) 'RDM_energy at rdm_finalising.F90 ', RDM_energy
end if
! this is allocated in find_nat_orb_occ_numbers and used later in
! brokensymno, ugh. Have to deallocate it somewhere though
if (allocated(FourIndInts)) deallocate(FourIndInts)
call halt_timer(FinaliseRDMs_Time)
end subroutine finalise_rdms