#include "macros.h" ! Routines to perform final operations on RDMs, which includes calculating ! 1-RDMs from 2-RDMs, apply symmetries, and writing RDMs out. module rdm_finalising use bit_rep_data, only: NIfTot use constants use Parallel_neci, only: iProcIndex, nProcessors use rdm_data, only: rdm_list_t, rdm_spawn_t, one_rdm_t use rdm_data_utils, only: calc_separate_rdm_labels, extract_sign_rdm, add_to_rdm_spawn_t use rdm_data_utils, only: communicate_rdm_spawn_t_wrapper use util_mod use CalcData, only: tAdaptiveShift use RotateOrbsMod, only: FourIndInts use SystemData, only: tGUGA, nSpatorbs use LoggingData, only: tWriteSpinFreeRDM, t_print_molcas_rdms, t_print_hdf5_rdms use matrix_util, only: print_matrix use guga_bitRepOps, only: extract_2_rdm_ind use guga_rdm, only: output_molcas_rdms implicit none real(dp), public, protected :: RDM_energy contains 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 subroutine output_2rdm_wrapper(rdm_defs, est, rdm, rdm_recv, rdm_recv_2, spawn) ! Call routines to output RDMs in all requested forms. ! We also calculate the Hermitian errors here, since this is something ! we typically want to do at the same point (the very end of a ! simulation usually), and this requires large parallel ! communications, as does the printing. use LoggingData, only: tWrite_normalised_RDMs, tWrite_RDMs_to_read use rdm_data, only: rdm_definitions_t, rdm_estimates_t, rdm_list_t, rdm_spawn_t, tOpenShell use rdm_estimators, only: calc_hermitian_errors use hash, only: clear_hash_table type(rdm_definitions_t), intent(in) :: rdm_defs type(rdm_estimates_t), intent(inout) :: est type(rdm_list_t), intent(in) :: rdm type(rdm_list_t), intent(inout) :: rdm_recv, rdm_recv_2 type(rdm_spawn_t), intent(inout) :: spawn ! Print the RDM popsfile, which requires no communication, first. This ! is in case we run out of memory during one of the other routines ! called here, which could crash if there isn't enough memory for ! communication! Then at least we can read the RDM back in... if (tWrite_RDMs_to_read) call print_rdm_popsfile(rdm) call calc_hermitian_errors(rdm, rdm_recv, spawn, est%norm, & est%max_error_herm, est%sum_error_herm) if (tGUGA) then spawn%free_slots = spawn%init_free_slots(0:nProcessors - 1) call clear_hash_table(spawn%rdm_send%hash_table) call make_hermitian_rdm(rdm, rdm_defs%nrdms_standard, spawn, rdm_recv) call print_spinfree_2rdm(rdm_defs, rdm_recv, est%norm) if (t_print_molcas_rdms) then call output_molcas_rdms(rdm_defs, rdm_recv, est%norm) end if else if (tWriteSpinFreeRDM) & call print_spinfree_2rdm_wrapper(rdm_defs, rdm, rdm_recv, spawn, est%norm) if (tWrite_Normalised_RDMs) & call print_rdms_spin_sym_wrapper(rdm_defs, rdm, rdm_recv, rdm_recv_2, & spawn, est%norm, tOpenShell) end if end subroutine output_2rdm_wrapper subroutine calc_1rdms_from_spinfree_2rdms(one_rdms, two_rdms, rdm_trace) ! For each 2-RDM in two_rdms, calculate the corresponding spin-free ! 1-RDM: ! ! \gamma^{spinfree}_{p,q} = \frac{1}{N-1} \sum_a \Gamma^{spinfree}_{pa,qa} ! ! Here, p, q and a are spatial labels. N is the number of electrons. ! The spinfree 1-RDM is defined in terms of the spinned 1-RDM by: ! ! \gamma^{spinfree}_{p,q} = \gamma_{p\alpha,q\alpha) + \gamma_{p\beta,q\beta) ! IMPORTANT: This routine should *only* be used by taking *spin-free* ! 2-RDMs as the input. Specifically, it takes spin-free RDMs as returned ! by the create_spinfree_2rdm routine, which does *not* restrict the ! labels allowed. Inputting 2-RDMs in other will give incorrect results! ! The output 1-RDM elements are sorted in the standard form: elements ! are indexed using the SymLabelListInv_rot array, so that the 1-RDMs ! will be in block-diagonal form, with elements within each symmetry ! block stored together. use Parallel_neci, only: MPISumAll use RotateOrbsData, only: SymLabelListInv_rot use SystemData, only: nel use UMatCache, only: spatial type(one_rdm_t), intent(inout) :: one_rdms(:) type(rdm_list_t), intent(in) :: two_rdms real(dp), intent(in) :: rdm_trace(:) integer(int_rdm) :: pqrs integer :: pq, rs, p, q, r, s ! spatial orbitals integer :: ielem, irdm, ierr real(dp) :: rdm_sign(two_rdms%sign_length) real(dp), allocatable :: temp_rdm(:, :) do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix = 0.0_dp end do ! Loop over all elements of the 2-RDM, \Gamma^{spinfree}_{pq,rs}, where ! p, q, r and s are spatial labels. If at least two spatial indices are ! the same then we have a contribution to the 1-RDM. do ielem = 1, two_rdms%nelements pqrs = two_rdms%elements(0, ielem) ! Obtain spin orbital labels and the RDM element. if (tGUGA) then call extract_2_rdm_ind(pqrs, p, q, r, s) else call calc_separate_rdm_labels(pqrs, pq, rs, r, s, q, p) end if call extract_sign_rdm(two_rdms%elements(:, ielem), rdm_sign) associate(ind => SymLabelListInv_rot) if (tGUGA) then if (r == s) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(p), ind(q)) = & one_rdms(irdm)%matrix(ind(p), ind(q)) + rdm_sign(irdm) / 2.0_dp end do end if ! if I count both, I do not need the factor 2.. if (p == q) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(r), ind(s)) = & one_rdms(irdm)%matrix(ind(r), ind(s)) + rdm_sign(irdm) / 2.0_dp end do end if else ! An element of the form \Gamma_{pa,ra}. if (q == s) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(p), ind(r)) = & one_rdms(irdm)%matrix(ind(p), ind(r)) + rdm_sign(irdm) end do end if end if end associate end do ! Allocate a temporary array in which to receive the MPI communication. allocate(temp_rdm(size(one_rdms(1)%matrix, 1), size(one_rdms(1)%matrix, 2)), stat=ierr) ! Perform a sum over all processes, for each 1-RDM being sampled. do irdm = 1, size(one_rdms) call MPISumAll(one_rdms(irdm)%matrix, temp_rdm) ! Copy summed RDM back to the main array, and normalise. one_rdms(irdm)%matrix = temp_rdm / (rdm_trace(irdm) * real(nel - 1, dp)) end do deallocate(temp_rdm, stat=ierr) end subroutine calc_1rdms_from_spinfree_2rdms subroutine calc_1rdms_from_2rdms(rdm_defs, one_rdms, two_rdms, rdm_trace, open_shell) ! For each 2-RDM in two_rdms, if open_shell is true then calculate the ! full spinned 1-RDM, otherwise calculate the spinfree 1-RDM. The ! former case is defined by: ! ! \gamma_{i,j} = \frac{1}{N-1} \sum_k \Gamma_{ik,jk} ! ! Here, i, j and k are spatial labels. N is the number of electrons. ! ! The spinfree case is then a contraction over the spin labels of the ! spinned 1-RDM: ! ! \gamma^{spinfree}_{p,q} = \gamma_{p\alpha,q\alpha) + \gamma_{p\beta,q\beta) ! ! where p and q are spatial labels. ! The output 1-RDM elements are sorted in the standard form: elements ! are indexed using the SymLabelListInv_rot array, so that the 1-RDMs ! will be in block-diagonal form, with elements within each symmetry ! block stored together. use Parallel_neci, only: MPISumAll use rdm_data, only: rdm_definitions_t use RotateOrbsData, only: SymLabelListInv_rot use SystemData, only: nel use UMatCache, only: spatial type(rdm_definitions_t), intent(in) :: rdm_defs type(one_rdm_t), intent(inout) :: one_rdms(:) type(rdm_list_t), intent(in) :: two_rdms real(dp), intent(in) :: rdm_trace(:) logical, intent(in) :: open_shell integer(int_rdm) :: ijkl integer :: ielem, irdm, ierr integer :: ij, kl, i, j, k, l ! spin orbitals integer :: p, q, r, s ! spatial orbitals real(dp) :: rdm_sign(two_rdms%sign_length) real(dp), allocatable :: temp_rdm(:, :) logical :: is_transition_rdm #ifdef DEBUG_ character(*), parameter :: this_routine = "calc_1rdms_from_2rdms" #endif ASSERT(.not. tGUGA) do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix = 0.0_dp end do ! Loop over all elements of the 2-RDM, \Gamma_{pq,rs}, where p, q, r ! and s are spatial labels. If at least two spatial indices are the ! same then we have a contribution to the 1-RDM. do ielem = 1, two_rdms%nelements ijkl = two_rdms%elements(0, ielem) ! Obtain spin orbital labels and the RDM element. call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l) ! For closed shell systems we work with spatial orbitals, to ! calculate spin-free 1RDMs. if (open_shell) then p = i; q = j; r = k; s = l; else p = spatial(i); q = spatial(j); r = spatial(k); s = spatial(l); end if call extract_sign_rdm(two_rdms%elements(:, ielem), rdm_sign) ! If abba or baab term - swap last two indices and sign. if (.not. same_spin(i, k)) then if (open_shell) then r = l; s = k; else r = spatial(l); s = spatial(k); end if rdm_sign = -rdm_sign end if associate(ind => SymLabelListInv_rot) ! An element of the form \Gamma_{aq,as}. if (p == r) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(q), ind(s)) = one_rdms(irdm)%matrix(ind(q), ind(s)) + rdm_sign(irdm) end do end if ! An element of the form \Gamma_{pa,ra}. if (q == s) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(p), ind(r)) = one_rdms(irdm)%matrix(ind(p), ind(r)) + rdm_sign(irdm) end do end if ! The below cases give contributions by swapping one pair of ! indices. Only include these contributions if we have aaaa or ! bbbb terms. This because if we had a term with spin signature ! abab (for example), then swapping as below would give abba ! or baab terms, which don't contribute to the 1-RDM. if (same_spin(k, l)) then ! An element of the form \Gamma_{pa,as}. if (p == s) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(q), ind(r)) = one_rdms(irdm)%matrix(ind(q), ind(r)) - rdm_sign(irdm) end do end if ! An element of the form \Gamma_{aq,ra}. if (q == r) then do irdm = 1, size(one_rdms) one_rdms(irdm)%matrix(ind(p), ind(s)) = one_rdms(irdm)%matrix(ind(p), ind(s)) - rdm_sign(irdm) end do end if end if end associate end do ! Allocate a temporary RDM array. allocate(temp_rdm(size(one_rdms(1)%matrix, 1), size(one_rdms(1)%matrix, 2)), stat=ierr) ! Make non-transition RDMs symmetric. do irdm = 1, size(one_rdms) is_transition_rdm = rdm_defs%state_labels(1, irdm) /= rdm_defs%state_labels(2, irdm) if (.not. is_transition_rdm) then ! Use temp_rdm as temporary space for the transpose, to (hopefully) ! prevent a temporary array being created in the sum below. temp_rdm = transpose(one_rdms(irdm)%matrix) one_rdms(irdm)%matrix = (one_rdms(irdm)%matrix + temp_rdm) / 2.0_dp end if end do ! Perform a sum over all processes, for each 1-RDM being sampled. do irdm = 1, size(one_rdms) call MPISumAll(one_rdms(irdm)%matrix, temp_rdm) ! Copy summed RDM back to the main array, and normalise. one_rdms(irdm)%matrix = temp_rdm / (rdm_trace(irdm) * real(nel - 1, dp)) end do deallocate(temp_rdm, stat=ierr) end subroutine calc_1rdms_from_2rdms subroutine make_hermitian_rdm(rdm, nrdms_standard, spawn, rdm_recv) ! Take the RDM in the rdm object, and output a new RDM object which is ! the same but with hermiticy applied to the non-transition RDMs, i.e., ! the elements above and below the diagonal are averaged appropriately. ! Transition RDM signs are set to zero, so transition RDMs are not ! included in the output RDM. use rdm_data_utils, only: annihilate_rdm_list type(rdm_list_t), intent(in) :: rdm integer, intent(in) :: nrdms_standard type(rdm_spawn_t), intent(inout) :: spawn type(rdm_list_t), intent(inout) :: rdm_recv integer(int_rdm) :: ijkl integer :: ielem, ij, kl, i, j, k, l real(dp) :: rdm_sign(rdm%sign_length) logical :: nearly_full, finished, all_finished #ifdef DEBUG_ character(*), parameter :: this_routine = "make_hermitian_rdm" #endif integer(int_rdm) :: ij_, kl_ ! If we're about to fill up the spawn list, perform a communication. nearly_full = .false. ! Have we finished adding RDM elements to the spawned list? finished = .false. rdm_recv%nelements = 0 do ielem = 1, rdm%nelements ! If the spawned list is nearly full, perform a communication. if (nearly_full) then call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) nearly_full = .false. end if ijkl = rdm%elements(0, ielem) call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign) ! Set sign for transition RDMs to 0. rdm_sign(nrdms_standard + 1:) = 0.0_dp if (tGUGA) then call extract_2_rdm_ind(ijkl, i, j, k, l, ij_, kl_) ij = int(ij_) kl = int(kl_) rdm_sign = rdm_sign / 2.0_dp call add_to_rdm_spawn_t(spawn, k, l, i, j, rdm_sign, .true., nearly_full) call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .true., nearly_full) else ! Obtain spin orbital labels and the RDM element. call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l) ! Factor of a half to account for prevent double-counting, and ! instead average elements from above and below the diagonal. if (ij /= kl) rdm_sign = 0.5_dp * rdm_sign ! If in the lower half of the RDM, reflect to the upper half. if (ij > kl) then call add_to_rdm_spawn_t(spawn, k, l, i, j, rdm_sign, .false., nearly_full) else call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .false., nearly_full) end if end if end do finished = .true. ! Keep performing communications until all RDM spawnings on every ! processor have been communicated. do call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) if (all_finished) exit end do call annihilate_rdm_list(rdm_recv) end subroutine make_hermitian_rdm subroutine apply_symmetries_for_output(rdm, rdm_recv, spawn, open_shell) ! This routine will take in rdm, and output a new rdm which will have ! all appropriate symmetries applied so that the latter RDM can be ! passed to the routine to write RDMs. ! The input RDM must already have hermiticy symmetry applied to it. ! WARNING: To clarify potential confusion, we point out that this ! routine also applies hermiticy again, but *only* for the spatial ! labels, not the full spin labels. It does so specifically using a ! particular legacy ordering. This is not a mistake - for open shell ! systems we do need to apply full hermiticy, but also need to apply spatial ! hermiticy because spatial labels are written within each output file. use rdm_data_utils, only: annihilate_rdm_list type(rdm_list_t), intent(in) :: rdm type(rdm_list_t), intent(inout) :: rdm_recv type(rdm_spawn_t), intent(inout) :: spawn logical, intent(in) :: open_shell integer(int_rdm) :: ijkl integer :: ielem integer :: ij, kl, i, j, k, l ! spin orbitals integer :: p, q, r, s ! spatial orbitals integer :: pq_legacy, rs_legacy ! spatial orbitals real(dp) :: rdm_sign(rdm%sign_length) logical :: nearly_full, finished, all_finished #ifdef DEBUG_ character(*), parameter :: this_routine = "apply_symmetries_for_output" #endif ASSERT(.not. tGUGA) ! If we're about to fill up the spawn list, perform a communication. nearly_full = .false. ! Have we finished adding RDM elements to the spawned list? finished = .false. rdm_recv%nelements = 0 do ielem = 1, rdm%nelements ! If the spawned list is nearly full, perform a communication. if (nearly_full) then call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) nearly_full = .false. end if ijkl = rdm%elements(0, ielem) ! Obtain spin orbital labels and the RDM element. call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l) call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign) ! When there are two elements which are guaranteed to be exactly the ! same, we usually only want to print one of them (and to average ! over the equal terms). This function returns the labels we want. call apply_legacy_output_ordering(i, j, k, l, rdm_sign, pq_legacy, rs_legacy) ! For closed shell systems, want bbbb -> aaaa, baba -> abab, ! baab -> abba, by flipping all spins, which is a symmetry for such ! systems. If the first label has beta spin then we definitely want ! to move this RDM element to the equivalent flipped term so go ahead ! and flip all the spins. if (.not. open_shell) then if (is_beta(i)) then ! The ab_pair macro swaps alpha and beta spins of a label ! while keeping the spatial orbital unchanged. i = ab_pair(i) j = ab_pair(j) k = ab_pair(k) l = ab_pair(l) end if ! If the spatial parts of i and j are the same, and the spatial ! parts of k and l are *also* the same, then the RDM element won't ! have been added into both equivalent spin-flipped arrays ! because i<j and k<l is enforced), so we don't count twice. if (.not. (is_in_pair(i, j) .and. is_in_pair(k, l))) then ! Also, if (i,j) and (k,l) have the same spatial parts, but ! different spin parts ((alpha,beta) and (beta,alpha), or ! vice versa) then they only occur once, again because we ! enforce i<j and k<l for all stored RDM elements. if (.not. (pq_legacy == rs_legacy .and. (.not. ij == kl))) then rdm_sign = rdm_sign * 0.5_dp end if end if end if call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .false., nearly_full) if (open_shell) then ! For open shell systems, if i and j have the same spatial parts, ! and k and l do too, then we only have baba spin signature, ! (because we enforce i<j, k<l) but we'd like to print out abab too. if (is_in_pair(i, j) .and. is_in_pair(k, l)) then call add_to_rdm_spawn_t(spawn, j, i, l, k, rdm_sign, .false., nearly_full) end if ! Because we enforce hermiticy symmetry in the output, we would ! only print the following terms with baab. We want to print it ! with abba too here, so do that. if (pq_legacy == rs_legacy .and. (.not. ij == kl)) then call add_to_rdm_spawn_t(spawn, k, l, i, j, rdm_sign, .false., nearly_full) end if end if end do finished = .true. ! Keep performing communications until all RDM spawnings on every ! processor have been communicated. do call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) if (all_finished) exit end do call annihilate_rdm_list(rdm_recv) end subroutine apply_symmetries_for_output pure subroutine apply_legacy_output_ordering(i, j, k, l, rdm_sign, pq_legacy, rs_legacy) ! Enforce the symmetries of RDMs to only keep certain combinations of ! i, j, k, l spin labels, where a redundancy exists. Whenever we have ! an unused combination, flip/swap labels (and the sign if necessary). ! For example the 2-RDM is hermitian, so if ij /= kl, then we only need ! to print either \Gamma_{ij,kl} or \Gamma{kl,ij}, but not both. Which ! combinations we decide to print is decided below, which is purely a ! legacy decision (as far as I know!). See comments below for definitions ! of what we keep. use SystemData, only: nbasis use UMatCache, only: spatial integer, intent(inout) :: i, j, k, l ! spin orbitals real(dp), intent(inout) :: rdm_sign(:) integer, intent(out) :: pq_legacy, rs_legacy ! spatial orbitals integer :: i_temp, j_temp ! spin orbitals integer :: p, q, r, s ! spatial_orbitals ! RDMs are output in files labelled by their spin signatures: ! aaaa, abab, abba, bbbb, baba or baab. ! Within each file, therefore, only spatial orbital labels are printed. ! Thus, we need to use spatial orbitals to determine which RDM elements ! are to kept, and which transformed. p = spatial(i); q = spatial(j); r = spatial(k); s = spatial(l); ! When we calculate the combined labels, pq and rs, we would ! usually have p and q swapped below, and similarly with r and s. ! However, the old RDM files prints only RDM elements with pq < rs, ! where pq and rs are defined as follows. pq_legacy = (q - 1) * nbasis + p rs_legacy = (s - 1) * nbasis + r ! Apply symmetry (for *real* RDMs), to only print elements from one ! half of the RDM, using the legacy ordering. if (pq_legacy > rs_legacy) then i_temp = i; j_temp = j; i = k; j = l; k = i_temp; l = j_temp; end if ! If either i and j have the same spatial part, of k and l have the ! same spatial part, and we have a spin signature with 2 alphas and ! 2 betas, then the convention is to output it as either abab or ! baba, but *not* as abba or baab. If we have abba or baab in this ! case then we have to swap two indices and introduce a minus sign. ! Because we enforce i<j and k<l in all RDM elements, there are only ! two possibilities to consider: if (is_in_pair(i, j) .and. is_beta(i) .and. is_alpha(j) .and. & is_alpha(k) .and. is_beta(l)) then i = ab_pair(i) j = ab_pair(j) rdm_sign = -rdm_sign else if (is_in_pair(k, l) .and. is_alpha(i) .and. is_beta(j) .and. & is_beta(k) .and. is_alpha(l)) then k = ab_pair(k) l = ab_pair(l) rdm_sign = -rdm_sign end if end subroutine apply_legacy_output_ordering subroutine print_rdms_spin_sym_wrapper(rdm_defs, rdm, rdm_recv, rdm_recv_2, spawn, rdm_trace, open_shell) ! Compress the full spinned-RDMs by summing over spin-equivalent terms ! (i.e. aaaa and bbbb rdms), and also applying symmetry of (*real*) ! RDMs. The result will be stored in rdm_recv_2. Then, print it out to a ! file. ! IMPORTANT: Although the rdm object has inout status, it will *not* ! be modified. The inout status is to allow for the optional possibility ! of updating the first argument of make_hermitian_rdm, which is not ! used here. use hash, only: clear_hash_table use rdm_Data, only: rdm_definitions_t use LoggingData, only: t_calc_double_occ use double_occ_mod, only: calc_double_occ_from_rdm type(rdm_definitions_t), intent(in) :: rdm_defs type(rdm_list_t), intent(in) :: rdm type(rdm_list_t), intent(inout) :: rdm_recv, rdm_recv_2 type(rdm_spawn_t), intent(inout) :: spawn real(dp), intent(in) :: rdm_trace(rdm%sign_length) logical, intent(in) :: open_shell integer :: nrdms_to_print spawn%free_slots = spawn%init_free_slots(0:nProcessors - 1) call clear_hash_table(spawn%rdm_send%hash_table) call make_hermitian_rdm(rdm, rdm_defs%nrdms_standard, spawn, rdm_recv) call apply_symmetries_for_output(rdm_recv, rdm_recv_2, spawn, open_shell) ! Only print non-transition RDMs, for now. nrdms_to_print = rdm_defs%nrdms_standard call print_rdms_with_spin(rdm_defs, nrdms_to_print, rdm_recv_2, rdm_trace, open_shell) ! intermediate hack: if (t_calc_double_occ) then call calc_double_occ_from_rdm(rdm_recv_2, rdm_trace) end if end subroutine print_rdms_spin_sym_wrapper subroutine create_spinfree_2rdm(rdm, nrdms_standard, spawn, rdm_recv) ! Take an standard (spinned) 2-RDM, stored in rdm, and output the ! spinfree version of it to the spawn%rdm_recv object. ! The input RDM has elements equal to: ! ! \Gamma_{ij,kl} = < a^+_i a^+_j a_l a_k > ! ! where i, j, k and l are spin orbital labels, and the output spinfree ! RDM has elements equal to: ! ! \Gamma^{spinfree}_{pq,rs} = \sum_{x,y} < a^+_{p,x} a^+_{q,y} a_{s,y} a_{r,x} > ! ! where p, q, r, s are spatial orbital labels, and x and y are spin ! labels (alpha or beta) which are summed over. ! Thus, all terms with spin signature aaaa, abab, bbbb or baba are ! summed together. Terms with spin signature abba or baab have their ! final two spin orbital labels swapped (introducing a minus sign), so ! that they give a contribution to the resulting spinfree RDM element. use rdm_data_utils, only: annihilate_rdm_list use SystemData, only: nbasis use UMatCache, only: spatial type(rdm_list_t), intent(in) :: rdm integer, intent(in) :: nrdms_standard type(rdm_spawn_t), intent(inout) :: spawn type(rdm_list_t), intent(inout) :: rdm_recv integer(int_rdm) :: ijkl integer :: ielem integer :: ij, kl, i, j, k, l, k_orig, l_orig ! spin orbitals integer :: pq, rs, p, q, r, s ! spatial orbitals real(dp) :: rdm_sign(rdm%sign_length) logical :: nearly_full, finished, all_finished ! If we're about to fill up the spawn list, perform a communication. nearly_full = .false. ! Have we finished adding RDM elements to the spawned list? finished = .false. rdm_recv%nelements = 0 do ielem = 1, rdm%nelements ! If the spawned list is nearly full, perform a communication. if (nearly_full) then call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) nearly_full = .false. end if ijkl = rdm%elements(0, ielem) ! Obtain spin orbital labels and the RDM element. call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l) call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign) ! Store the original labels, before we possibly swap them. k_orig = k; l_orig = l; ! If this term is abba or baab then we can make it abab or baba by ! swapping the last two indices, which introduces a minus sign. ! It will then contribute to a spinfree 2-RDM element. if (.not. same_spin(i, k)) then l = k_orig k = l_orig rdm_sign = -rdm_sign end if ! Get the spatial orbital labels from the spin orbital ones. p = spatial(i); q = spatial(j); r = spatial(k); s = spatial(l); ! The 'combined' labels. pq = (p - 1) * nbasis + q rs = (r - 1) * nbasis + s ! If the RDM is not symmetrised then the same term will be added ! from both above below the diagonal, so in this case we want a ! factor of a half to average and not double count. However, ! we do not apply hermiticty to transition RDMs, as they are not ! hermitian. So only apply this to non-transition RDMs. if (pq /= rs) rdm_sign(1:nrdms_standard) = rdm_sign(1:nrdms_standard) * 0.5_dp ! Due to the fact that RDM elements are only stored with p < q and ! r < s, the following terms are only stored with baba spin, never ! with abab. Double this term to make up for it. if (p == q .and. r == s) rdm_sign = 2.0_dp * rdm_sign ! Add all spinfree 2-RDM elements corresponding to these labels. call add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full) ! If this is an aaaa or bbbb term then *minus* this RDM element will ! be equal to the equivalent RDM element with the last two labels ! swapped. So, add this contribution into that RDM element. We ! don't have to do this, but doing so applies some extra averaging. ! Want to apply all the averaging possible over equivalent elements. if (same_spin(i, j)) then ! Re-extract sign in case it has been modified. call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign) ! Swap the spatial labels. r = spatial(l_orig); s = spatial(k_orig); rs = (r - 1) * nbasis + s ! Half the sign of non-transition RDMs, where we apply ! hermiticity to average above and below the diagonal, to ! prevent double counting. if (pq /= rs) rdm_sign(1:nrdms_standard) = rdm_sign(1:nrdms_standard) * 0.5_dp rdm_sign = -rdm_sign call add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full) end if end do finished = .true. ! Keep performing communications until all RDM spawnings on every ! processor have been communicated. do call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished) if (all_finished) exit end do call annihilate_rdm_list(rdm_recv) contains subroutine add_rdm_elements(p, q, r, s, nrdms_standard, rdm_sign, spawn, nearly_full) ! Add in the single contribution rdm_sign to the following elements ! of the spinfree 2-RDM: ! ! \Gamma^{spinfree}_{pq,rs} = \sum_{x,y} < a^+_{p,x} a^+_{q,y} a_{s,y} a_{r,x} > ! \Gamma^{spinfree}_{qp,sr} = \sum_{x,y} < a^+_{q,x} a^+_{p,y} a_{r,y} a_{s,x} > ! ! And for non-transition RDMs, which are hermitian, also add in ! the following elements: ! ! \Gamma^{spinfree}_{rs,pq} = \sum_{x,y} < a^+_{r,x} a^+_{s,y} a_{q,y} a_{p,x} > ! \Gamma^{spinfree}_{sr,qp} = \sum_{x,y} < a^+_{s,x} a^+_{r,y} a_{p,y} a_{q,x} > ! ! where x and y are spin labels which are summed over in the final ! result. ! ! For a *REAL* spinfree 2-RDM, all of these elements are rigorously ! equal, so it is appropriate that we add all contributions in ! together like this (the last two aren't equal to the first two for ! transition RDMs, so don't add these for transition RDMs). ! ! The if-statements in here prevent adding to the same RDM element ! twice. integer, intent(in) :: p, q, r, s ! spatial orbitals integer, intent(in) :: nrdms_standard real(dp), intent(in) :: rdm_sign(:) type(rdm_spawn_t), intent(inout) :: spawn logical, intent(inout) :: nearly_full real(dp) :: rdm_sign_temp(size(rdm_sign)) rdm_sign_temp = 0.0_dp ! RDM element \Gamma_{pq,rs}. call add_to_rdm_spawn_t(spawn, p, q, r, s, rdm_sign, .true., nearly_full) ! RDM element \Gamma_{qp,sr}. if (.not. (p == q .and. r == s)) then call add_to_rdm_spawn_t(spawn, q, p, s, r, rdm_sign, .true., nearly_full) end if if (pq /= rs) then ! We only want to apply hermiticity symmetry to non-transition ! RDMs, since transition RDMs are not hermitian. rdm_sign_temp(1:nrdms_standard) = rdm_sign(1:nrdms_standard) ! RDM element \Gamma_{rs,pq}. call add_to_rdm_spawn_t(spawn, r, s, p, q, rdm_sign_temp, .true., nearly_full) ! RDM element \Gamma_{sr,qp}. if (.not. (p == q .and. r == s)) then call add_to_rdm_spawn_t(spawn, s, r, q, p, rdm_sign_temp, .true., nearly_full) end if end if end subroutine add_rdm_elements end subroutine create_spinfree_2rdm subroutine print_spinfree_2rdm_wrapper(rdm_defs, rdm, rdm_recv, spawn, rdm_trace) ! A wrapper function to take in an RDM object, in a form as accumulated ! during an FCIQMC simulation, and print it in a form consistent for ! MPQC to read in. To do this, it must call the additional routine, ! create_spinfree_2rdm, to make the spinfree 2-RDM itself (see that ! routine for a definition). The spinfree 2-RDM is stored in rdm_recv, ! which is then printed to a file. use hash, only: clear_hash_table use rdm_data, only: rdm_definitions_t type(rdm_definitions_t), intent(in) :: rdm_defs type(rdm_list_t), intent(in) :: rdm type(rdm_list_t), intent(inout) :: rdm_recv type(rdm_spawn_t), intent(inout) :: spawn real(dp), intent(in) :: rdm_trace(rdm%sign_length) spawn%free_slots = spawn%init_free_slots(0:nProcessors - 1) call clear_hash_table(spawn%rdm_send%hash_table) call create_spinfree_2rdm(rdm, rdm_defs%nrdms_standard, spawn, rdm_recv) call print_spinfree_2rdm(rdm_defs, rdm_recv, rdm_trace) end subroutine print_spinfree_2rdm_wrapper subroutine print_rdm_popsfile(rdm) ! Print the RDM object stored in rdm to a file called RDM_POPSFILE. ! This is done in a binary form, which can be quickly read in a ! subsequent NECI simulation (see read_2rdm_popsfile). The RDM is ! printed directly in its int_rdm form, rather than decoding signs ! to their real(dp) form. use Parallel_neci, only: MPIBarrier use util_mod, only: get_free_unit type(rdm_list_t), intent(in) :: rdm integer :: ielem, iproc, ierr, pops_unit do iproc = 0, nProcessors - 1 if (iproc == iProcIndex) then ! Let the first processor clear the file to start with. if (iproc == 0) then pops_unit = get_free_unit() open(pops_unit, file='RDM_POPSFILE', status='replace', form='unformatted') ! Let the first processor start by printing the number of ! RDMs being sampled. write (pops_unit) rdm%sign_length else pops_unit = get_free_unit() open(pops_unit, file='RDM_POPSFILE', status='old', position='append', form='unformatted') end if do ielem = 1, rdm%nelements write(pops_unit) rdm%elements(:, ielem) end do close(pops_unit) end if ! Wait for the current processor to finish printing its RDM elements. call MPIBarrier(ierr) end do end subroutine print_rdm_popsfile 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 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 ! ------- Routines for finalising 1-RDMs --------------------------------- subroutine finalise_1e_rdm(rdm_defs, one_rdms, norm_1rdm) ! This routine takes the 1-RDM (matrix), sums it across processors, ! normalises it, makes it hermitian, and prints out the 'popsfile' ! 1-RDM, i.e. the 1-RDM before it is normalised or made Hermitian. use LoggingData, only: twrite_RDMs_to_read, tForceCauchySchwarz use LoggingData, only: RDMExcitLevel use Parallel_neci, only: iProcIndex, MPISumAll use rdm_data, only: rdm_definitions_t, one_rdm_t use RotateOrbsData, only: NoOrbs type(rdm_definitions_t), intent(in) :: rdm_defs type(one_rdm_t), intent(inout) :: one_rdms(:) real(dp), intent(out) :: norm_1rdm(size(one_rdms)) integer :: irdm, ierr real(dp) :: SumN_Rho_ii(size(one_rdms)) real(dp), allocatable :: AllNode_one_rdm(:, :) if (RDMExcitLevel == 1) then allocate(AllNode_one_rdm(NoOrbs, NoOrbs), stat=ierr) do irdm = 1, size(one_rdms) call MPISumAll(one_rdms(irdm)%matrix, AllNode_one_rdm) one_rdms(irdm)%matrix = AllNode_one_rdm end do deallocate(AllNode_one_rdm) end if ! Find the normalisation. call calc_1e_norms(rdm_defs, one_rdms, norm_1rdm) if (iProcIndex == 0) then do irdm = 1, size(one_rdms) ! Write out the unnormalised, non-hermitian OneRDM_POPS. if (twrite_RDMs_to_read) call write_1rdm(rdm_defs, one_rdms(irdm)%matrix, irdm, norm_1rdm(irdm), .false.) end do if (RDMExcitLevel == 1) then ! Only non-transition RDMs should be hermitian and obey the ! Cauchy-Schwarz inequality. do irdm = 1, rdm_defs%nrdms_standard call make_1e_rdm_hermitian(one_rdms(irdm)%matrix, norm_1rdm(irdm)) if (tForceCauchySchwarz) call Force_Cauchy_Schwarz(one_rdms(irdm)%matrix) end do end if end if end subroutine finalise_1e_rdm subroutine Force_Cauchy_Schwarz(matrix) use RotateOrbsData, only: SymLabelListInv_rot use SystemData, only: nbasis real(dp), intent(inout) :: matrix(:, :) integer :: i, j real(dp) :: UpperBound write(stdout, '("Ensuring that Cauchy--Schwarz inequality holds.")') associate(ind => SymLabelListInv_rot) do i = 1, nbasis do j = 1, nbasis UpperBound = sqrt(matrix(ind(i), ind(i)) * matrix(ind(j), ind(j))) if (abs(matrix(ind(i), ind(j))) > UpperBound) then if (matrix(ind(i), ind(j)) < 0.0_dp) then matrix(ind(i), ind(j)) = -UpperBound else if (matrix(ind(i), ind(j)) > 0.0_dp) then matrix(ind(i), ind(j)) = UpperBound end if write(stdout, '("Changing element:")') i, j else cycle end if end do end do end associate end subroutine Force_Cauchy_Schwarz subroutine calc_1e_norms(rdm_defs, one_rdms, norm_1rdm) ! Calculate and return the factor by which the input 1-RDMs need to be ! divided by, in order to correctly normalise them. The transition ! RDMs need to be normalised differently - see comments below. use rdm_data, only: rdm_definitions_t use RotateOrbsData, only: NoOrbs use SystemData, only: nel type(rdm_definitions_t), intent(in) :: rdm_defs type(one_rdm_t), intent(in) :: one_rdms(:) real(dp), intent(out) :: norm_1rdm(:) integer :: irdm, i real(dp) :: trace_1rdm(size(one_rdms)) trace_1rdm = 0.0_dp norm_1rdm = 0.0_dp do irdm = 1, size(one_rdms) do i = 1, NoOrbs trace_1rdm(irdm) = trace_1rdm(irdm) + one_rdms(irdm)%matrix(i, i) end do end do ! The non-transition RDMs must be normalised to have a trace of nel. do irdm = 1, rdm_defs%nrdms_standard norm_1rdm(irdm) = real(nel, dp) / trace_1rdm(irdm) end do ! The transition RDMs should be normalised using the normalisation of ! the contributing wave functions, which can be calculated from the ! traces of the corresponding RDMs. do irdm = rdm_defs%nrdms_standard + 1, size(one_rdms) norm_1rdm(irdm) = sqrt(norm_1rdm(rdm_defs%state_labels(1, irdm)) * norm_1rdm(rdm_defs%state_labels(2, irdm))) end do end subroutine calc_1e_norms subroutine calc_rho_ii_and_sum_n(one_rdms, norm_1rdm, SumN_Rho_ii) ! Calculate the rho_ii arrays, which hold the diagonals of the 1-RDMs, ! appropriately normalised and sorted in terms of the energy ordering ! of the orbitals. Then, calculate SumN_Rho_ii, which holds the sum of ! this diagonal for orbitals occupied in the HF determinant. use FciMCData, only: HFDet_True use LoggingData, only: tDiagRDM use rdm_data, only: tOpenShell use RotateOrbsData, only: SymLabelListInv_rot, NoOrbs use SystemData, only: BRR, nel use UMatCache, only: gtID type(one_rdm_t), intent(inout) :: one_rdms(:) real(dp), intent(in) :: norm_1rdm(:) real(dp), intent(out) :: SumN_Rho_ii(:) integer :: irdm, i, HFDet_ID, BRR_ID ! Want to sum the diagonal elements of the 1-RDM for the HF orbitals. ! Given the HF orbitals, SymLabelListInv_rot tells us their position ! in the 1-RDM. SumN_Rho_ii = 0.0_dp do irdm = 1, size(one_rdms) do i = 1, NoOrbs ! rho_ii is the diagonal elements of the 1-RDM. We want this ! ordered according to the energy of the orbitals. Brr has the ! orbital numbers in order of energy... i.e Brr(2) = the orbital ! index with the second lowest energy. BRR is always in spin ! orbitals. i gives the energy level, BRR gives the orbital, ! SymLabelListInv_rot gives the position of this orbital in ! one_rdm. associate(ind => SymLabelListInv_rot) if (tDiagRDM) then if (tOpenShell) then one_rdms(irdm)%rho_ii(i) = one_rdms(irdm)%matrix(ind(BRR(i)), ind(BRR(i))) * norm_1rdm(irdm) else BRR_ID = gtID(BRR(2 * i)) one_rdms(irdm)%rho_ii(i) = one_rdms(irdm)%matrix(ind(BRR_ID), ind(BRR_ID)) * norm_1rdm(irdm) end if end if if (i <= nel) then if (tOpenShell) then SumN_Rho_ii(irdm) = SumN_Rho_ii(irdm) + & (one_rdms(irdm)%matrix(ind(HFDet_True(i)), ind(HFDet_True(i))) * norm_1rdm(irdm)) else HFDet_ID = gtID(HFDet_True(i)) SumN_Rho_ii(irdm) = SumN_Rho_ii(irdm) + & (one_rdms(irdm)%matrix(ind(HFDet_ID), ind(HFDet_ID)) * norm_1rdm(irdm)) / 2.0_dp end if end if end associate end do end do end subroutine calc_rho_ii_and_sum_n subroutine make_1e_rdm_hermitian(matrix, norm_1rdm) ! Simply average the 1-RDM(i,j) and 1-RDM(j,i) elements which should ! be equal in a perfect world. use RotateOrbsData, only: SymLabelListInv_rot, NoOrbs real(dp), intent(inout) :: matrix(:, :) real(dp), intent(in) :: norm_1rdm real(dp) :: max_error_herm, sum_error_herm integer :: i, j real(dp) :: temp max_error_herm = 0.0_dp sum_error_herm = 0.0_dp associate(ind => SymLabelListInv_rot) do i = 1, NoOrbs do j = i, NoOrbs if ((abs((matrix(ind(i), ind(j)) * norm_1rdm) - (matrix(ind(j), ind(i)) * norm_1rdm))) > max_error_herm) then max_error_herm = abs(matrix(ind(i), ind(j)) * norm_1rdm - matrix(ind(j), ind(i)) * norm_1rdm) end if sum_error_herm = sum_error_herm + abs(matrix(ind(i), ind(j)) * norm_1rdm - matrix(ind(j), ind(i)) * norm_1rdm) temp = (matrix(ind(i), ind(j)) + matrix(ind(j), ind(i))) / 2.0_dp matrix(ind(i), ind(j)) = temp matrix(ind(j), ind(i)) = temp end do end do end associate ! Output the hermiticity errors. write(stdout, '(1X,"MAX ABS ERROR IN 1-RDM HERMITICITY",F20.13)') max_error_herm write(stdout, '(1X,"MAX SUM ERROR IN 1-RDM HERMITICITY",F20.13)') sum_error_herm end subroutine make_1e_rdm_hermitian 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 end module rdm_finalising