#include "macros.h"
module LMat_freeze
use IntegralsData, only: nFrozen, UMat, umat_win
use UMatCache, only: numBasisIndices, UMatInd
use constants
use util_mod, only: operator(.div.), custom_findloc
use OneEInts, only: TMat2D
use SystemData, only: ECore
use Parallel_neci
use MPI_wrapper
use mpi, only: MPI_DOUBLE_PRECISION
implicit none
private
public :: freeze_lmat, t_freeze, map_indices, init_freeze_buffers, &
flush_freeze_buffers, add_core_en
!> Parameter for number of double excitations a single LMat entry can contribute to
!! (counting permutations)
integer, parameter :: num_ex = 4
integer, parameter :: num_inds = 6
! Step corresponding to conjugation of an index
integer, parameter :: step = (num_inds / 2)
! Local storage for ECore and TMat
real(dp) :: ECore_local, ECore_tot
HElement_t(dp), allocatable :: TMat_local(:, :), TMat_total(:, :)
! Direct orbital excitation = pair of indices in an index array which is apart by step
contains
!> Initialize the local storage for the diagonal and one-electron terms (two-electron
!! terms are shared memory). This is called before reading in the 6-index integrals
subroutine init_freeze_buffers()
! Setup the buffers for accumulating the correction per processor
if(nFrozen > 0) then
ECore_local = 0.0_dp
if(allocated(TMat_local)) deallocate(TMat_local)
allocate(TMat_local(size(TMat2D, dim=1), size(TMat2D, dim=2)))
if(allocated(TMat_total)) deallocate(TMat_total)
allocate(TMat_total(size(TMat2D, dim=1), size(TMat2D, dim=2)))
TMat_local = 0.0_dp
TMat_total = 0.0_dp
end if
end subroutine init_freeze_buffers
!------------------------------------------------------------------------------------------!
!> Sum the locally accumulated corrections to the diagonal and one-body terms up and
!! add them to the global terms, then deallocate temporaries. This is called
!! after reading in the 6-index integrals
subroutine flush_freeze_buffers()
integer(MPIArg) :: tmat_size
integer(MPIArg) :: ierr
if(nFrozen > 0) then
if(.not. allocated(TMat_local) .or. .not. allocated(TMat_total)) &
call stop_all("flush_freeze_buffers", &
"Buffers for freezing three-body interaction not allocated")
! Sum the core energy from each proc on this node
call MPI_Allreduce(ECore_local, ECore_tot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
mpi_comm_intra, ierr)
ECore = ECore + ECore_tot
tmat_size = size(TMat_local)
! Sum the 1-body terms from each proc on this node
call MPI_Allreduce(TMat_local, TMat_total, tmat_size, MPI_DOUBLE_PRECISION, MPI_SUM, &
mpi_comm_intra, ierr)
TMat2D = TMat2D+TMat_total
! Deallocate the temporary for the 1-body integrals
deallocate(TMat_local)
deallocate(TMat_total)
! The 2-body terms are shared memory, no operation is required at this point
end if
end subroutine flush_freeze_buffers
!------------------------------------------------------------------------------------------!
!> Checks if an entry is zeroed due to frozen orbitals being included
!> @param[in] indices array of size 6 containing the indices of the entry in question in the frozen orbital numbering
!> @return t_freeze true if the entry is zeroed
pure function t_freeze(indices)
integer(int64), intent(in) :: indices(num_inds)
logical :: t_freeze
t_freeze = any(indices < 1)
end function t_freeze
!------------------------------------------------------------------------------------------!
!> Maps a set of six indices from pre-freeze to post-freeze orbital indexing
!> @param indices on entry: array of indices in pre-freeze indexing, on return: same array in post-freeze indexing
subroutine map_indices(indices)
integer(int64), intent(inout) :: indices(:)
indices = indices - numBasisIndices(nFrozen)
end subroutine map_indices
!------------------------------------------------------------------------------------------!
!> Checks if the entry is neglected due to frozen orbitals being included and absorbs entries
!! into the lower order matrix elements if required
subroutine freeze_lmat(matel, indices)
HElement_t(dp), intent(inout) :: matel
integer(int64), intent(inout) :: indices(num_inds)
! Offset the orbital indexing
call map_indices(indices)
call add_core_en(matel, indices)
end subroutine freeze_lmat
!------------------------------------------------------------------------------------------!
!> Absorb entries with repeated frozen orbitals into the corresponding lower-order
!! terms.
subroutine add_core_en(matel, indices)
HElement_t(dp), intent(inout) :: matel
integer(int64), intent(in) :: indices(num_inds)
integer :: counts
integer(int64) :: idx(num_ex), ct
real(dp) :: prefactor, delta
integer(MPIArg) :: ierr
if(t_freeze(indices)) then
! Count the number of different indices appearing
counts = count_frozen_inds(indices)
! How many pairs of frozen orbitals do we have
select case(counts)
! 1 => double excitation
case(1)
idx = frozen_double_entry(indices, prefactor)
! There are up to eight potential double excitations with four given indices
! to which the LMat entry contributes - this is because LMat is hermitian
! but UMat is not
do ct = 1, num_ex
! If the index is assigned, there is a duplicate frozen orb
if(idx(ct) > 0) then
! Absorb the matrix element into UMat (with the corresponding prefactor
! according to permutation and number of possible spin terms
! UMat(idx(ct)) = UMat(idx(ct)) + prefactor*matel
delta = prefactor * matel
! Use a remote update of UMat
call MPI_Accumulate(delta, 1, MPI_DOUBLE_PRECISION, 0, &
idx(ct) - 1, 1, MPI_DOUBLE_PRECISION, MPI_SUM, umat_win, ierr)
endif
end do
! 2 => single excitation
case(2)
idx(1:2) = frozen_single_entry(indices, prefactor)
if(idx(1) > 0) then
! Absorb the matrix element into TMat
! TMat2D is indexed with spin orbs
call add_to_tmat(spatToSpinAlpha(idx(1)), spatToSpinAlpha(idx(2)))
call add_to_tmat(spatToSpinBeta(idx(1)), spatToSpinBeta(idx(2)))
endif
! 3 => diagonal element
case(3)
call frozen_diagonal_entry(indices, prefactor)
ECore_local = ECore_local + prefactor * matel
end select
! Zero the matrix element for further usage (i.e. will not turn up anymore)
matel = 0.0_dp
endif
contains
subroutine add_to_tmat(ind1, ind2)
integer(int64), intent(in) :: ind1, ind2
TMat_local(ind1, ind2) = TMat_local(ind1, ind2) + prefactor * matel
! TMat is hermitian, as is the 3-body term
if(ind1 /= ind2) &
TMat_local(ind2, ind1) = TMat_local(ind2, ind1) + prefactor * matel
end subroutine add_to_tmat
end subroutine add_core_en
!------------------------------------------------------------------------------------------!
!> Count how many frozen indices there are, determine if the entry has to be added
!! to a lower-order term and if yes, which one.
!> @param[in] indices array of size 6 with the orbital indices of the entry
!> @return level indicates to which lower order term this entry has to be added
!! and therefore, which subroutine has to be called in the following
!! 0 - none
!! 1 - one frozen index pair -> two-body term
!! 2 - two frozen index pairs -> one-body term
!! 3 - three frozen index pairs -> diagonal term
pure function count_frozen_inds(indices) result(level)
integer(int64), intent(in) :: indices(num_inds)
integer :: level
integer :: counts
integer :: ct
level = 0
! If there is an unpaired frozen orbital, the entry is discarded, it does not
! contribute to any contraction
do ct = 1, num_inds
if(count(indices(ct) == indices) == 1 .and. indices(ct) < 1) return
end do
! Else, it contributes to a contraction according to the number of frozen indices
counts = count(indices < 1)
! An even number of frozen orbital indices is required
if(modulo(counts, 2) /= 0) return
level = counts.div.2
end function count_frozen_inds
!------------------------------------------------------------------------------------------!
!> Get the index of the UMat entry to which the LMat entry with given indices shall be
!! added if it is frozen.
!> @param[in] indices array of lenght 6 with the orbital indices of the LMat entry
!> @param[ou] t_par flag indicating if the matrix element enters UMat with a -1
!> @return index index of the UMat entry to add the LMat entry to, 0 if entry is not frozen
function frozen_double_entry(indices, prefactor) result(idx)
integer(int64), intent(in) :: indices(num_inds)
real(dp), intent(out) :: prefactor
integer(int64) :: idx(num_ex)
integer :: f_one, f_two, ct, counts, i, rs, marks(3)
logical :: unfrozen(num_inds)
integer(int64) :: uf_idx(4)
idx = 0
unfrozen = indices > 0
! Mark where the unfrozen indices are
marks = 0
! Identify the direct unfrozen pair
! Look for the direct pair of unfrozen indices (these shall stay direct)
do ct = 1, step
if(unfrozen(ct) .and. unfrozen(ct + step)) then
uf_idx(1) = indices(ct)
uf_idx(3) = indices(ct + step)
! Mark the extracted indices so they are not taken again
marks(1) = ct
marks(2) = ct + step
exit
end if
end do
! Now, find the other two (direct or exchange) unfrozen indices
rs = 2
do ct = 1, num_inds
if(unfrozen(ct) .and. .not. any(ct == marks)) then
uf_idx(rs) = indices(ct)
marks(3) = ct
rs = rs + 2
endif
end do
! Now, get the required UMat indices, including all transposed index pairs
idx = permute_umat_inds(int(uf_idx(1)), int(uf_idx(2)), int(uf_idx(3)), int(uf_idx(4)))
! Check the permutation and possible factor of two due to spin
! First, get the two frozen orbitals
f_one = custom_findloc(unfrozen, .false., back=.false.)
f_two = custom_findloc(unfrozen, .false., back=.true.)
! There are two spin configurations in a close-shell frozen scenario for terms
! with a direct frozen orbital (i.e. no permutation required)
! All others enter with a prefactor of -1
if(is_direct(f_one, f_two)) then
prefactor = 2.0_dp
else
prefactor = -1.0_dp
! If this is a diagonal term (only three different indices), permutational symmetry
! gives an extra factor of 2 for the exchange term. If only two different
! indices appear, also the direct term gets this prefactor
counts = 1
! Count the number of different indices
do ct = 2, num_inds
if(.not. any(indices(ct) == indices(1:ct - 1))) counts = counts + 1
end do
if(counts == 2) then
prefactor = 2.0_dp * prefactor
elseif(counts == 3) then
! If there are three different indices, only add the extra factor if there is no direct unfrozen pair
do ct = 1, step
if(is_repeated_pair(indices, ct) .and. unfrozen(ct)) return
end do
prefactor = 2.0_dp * prefactor
end if
end if
end function frozen_double_entry
!------------------------------------------------------------------------------------------!
!> Returns the UMatInd values of all possible permutations of the input indices
!> @param[in] a,b,c,d orbital indices of a two-body element
!> @return inds array of size 4 containing the indices of UMat corresponding to these
!! four orbitals (in this order) and their hermitian conjugates
!! entries of 0 indicate that no position in UMat has to be addressed
function permute_umat_inds(a, b, c, d) result(inds)
integer, intent(in) :: a, b, c, d
integer(int64) :: inds(num_ex)
integer :: ct
inds = 0
! These are adjoint to each other, they are the same in LMat, but not UMat
inds(1) = UMatInd(a, b, c, d)
inds(2) = UMatInd(c, b, a, d)
inds(3) = UMatInd(a, d, c, b)
inds(4) = UMatInd(c, d, a, b)
! All terms are only counted if they are different from a previously occuring index
! This approach might be less performant than competing ways of doing this check,
! but it is very clear and less prone to bugs. If this proves to be a performance
! bottleneck, it can be changed to a faster implementation later on
do ct = 2, num_ex
! If the index already appeared, delete it
if(any(inds(ct) == inds(1:ct - 1))) inds(ct) = 0
end do
end function permute_umat_inds
!------------------------------------------------------------------------------------------!
!> For a given set of 6 orbital indices with two pairs of frozen orbitals,
!! returns the orbital indices of the corresponding single excitation and the prefactor
!! for the matrix element due to spin
!> @param[in] indices array of size 6 with orbital indices, four of which have to be
!! repeated frozen indices
!> @param[out] prefactor on return, the prefactor of the matrix element when added to
!! the one-body terms
!> @return orbs the two non-frozen orbitals
function frozen_single_entry(indices, prefactor) result(orbs)
integer(int64), intent(in) :: indices(num_inds)
real(dp), intent(out) :: prefactor
integer :: orbs(2)
integer :: f_orbs(num_inds - 2), ct
integer :: uf_one, uf_two, directs
logical :: unfrozen(num_inds)
orbs = 0
unfrozen = indices > 0
f_orbs = int(pack(indices,.not. unfrozen))
! For the single entry, there are only two unfrozen orbs, find these and return them
orbs = int(pack(indices, unfrozen))
! Again, two things to check: The number of direct orbitals and the permutation
! At this point, there is no unpaired frozen orbital, so there are two options:
! a) there are two different orbs
prefactor = 1.0_dp
if(count(f_orbs(1) == f_orbs) < size(f_orbs)) then
! We have to count how many potential spin configurations contribute
! The parity is counting the number of permutations required to end up
! with a completely direct expression
! -> The parity is odd if and only if the number of direct excitations is exactly one
! (then, one permutation is required to fix the two non-direct excits, else,
! (either zero or two are needed)
! -> count every direct excitation with a factor -1
directs = 0
! Count the number of direct excitations
do ct = 1, step
! Each direct repeated orbital doubles the number of spin configs
! (the case of four identical orbs is excluded)
if(is_repeated_pair(indices, ct) .or. (unfrozen(ct) .and. unfrozen(ct + step))) then
directs = directs + 1
endif
end do
select case(directs)
! If there are none, there is no freedom for spin choice, but the LMat entry
! appears twice in the matrix element evaluation due to symmetry if the
! unfrozen entry is repeated
case(0)
if(orbs(1) == orbs(2)) prefactor = 2.0_dp
! If there is one, two possible spin configs contribute, and we have odd parity
case(1)
prefactor = -2.0_dp
case(3)
! If there are three (just two is not possible), there are two spin degrees of
! freedom -> factor of 4
prefactor = 4.0_dp
end select
else
! b) All frozen orbs are the same -> only one spin config is allowed, only parity of
! the permutation counts
! Here, only a direct exctiation has even parity (same rule as above)
! Get the posisitons of the unfrozen orbs
uf_one = custom_findloc(unfrozen, .true., back=.false.)
uf_two = custom_findloc(unfrozen, .true., back=.true.)
if(.not. is_direct(uf_one, uf_two)) prefactor = -1.0_dp * prefactor
end if
end function frozen_single_entry
!------------------------------------------------------------------------------------------!
!> Determine the prefactor for a diagonal contribution (i.e. three pairs of frozen orbitals)
!> @param[in] indices orbital indices of the frozen entry - need to belong to a contribution
! to a diagonal matrix element of only frozen orbitals
!> @param[out] prefactor prefactor with which this LMat entry enters the matrix element
pure subroutine frozen_diagonal_entry(indices, prefactor)
integer(int64), intent(in) :: indices(num_inds)
real(dp), intent(out) :: prefactor
integer :: ct, directs
logical :: t_quad
! Only one info is needed here, that is the parity and number of contributing spin
! configs.
! If there are five or more repeated indices, the LMat entry does not contribute
! (it would require three same-spin electrons
t_quad = .false.
do ct = 1, 3
if(count(indices(ct) == indices) > 4) then
prefactor = 0.0_dp
return
end if
if(count(indices(ct) == indices) > 3) t_quad = .true.
end do
! Both relate to the number of different direct excits, so count these
directs = 0
do ct = 1, step
if(is_repeated_pair(indices, ct)) directs = directs + 1
end do
! If there is a quadruple index, the spin of it is fixed (both have to be occupied), then,
! the prefactor is always +/- 2, depending on if there is more than one direct pair
if(t_quad) then
! The only options with a quadruple index are
! a) one direct excitation => -2
! b) three direct exctiations => +2
prefactor = merge(-2.0_dp, 2.0_dp, directs == 1)
else
! In the other case, there are three relevant cases:
! a) All excitations are direct => prefactor 8 = 2**3 (all spins are free)
! b) One excitation is direct => prefactor -4 (from the spin freedom of the direct orb + overall spin)
! c) No excitation is direct => prefactor 4
! the factor of 4 comes from two sources: a factor of 2 from the overall spin
! and another factor of two from permutational symmetry: both indirect terms
! are identical in this case, so this entry has to be counted as both
if(directs == 0) then
prefactor = 4.0_dp
elseif(directs == 1) then
prefactor = -4.0_dp
else ! directs == 3
prefactor = 8.0_dp
endif
endif
end subroutine frozen_diagonal_entry
!------------------------------------------------------------------------------------------!
!> For two positions of indices in a 6-index set, return if these are a direct pair
!> @param[in] one, two integers between 1 and 6
!> @return t_dir true if the two integers correspond to positions that are direct, i.e.
!! for which the LMat entry is symmetric under exchange of the values
pure function is_direct(one, two) result(t_dir)
integer, intent(in) :: one, two
logical :: t_dir
t_dir = (one + step == two)
end function is_direct
!------------------------------------------------------------------------------------------!
!> For a set of 6 orbital indices, return whether the pair at a given position is repeated
!> @param[in] indices array of size 6 containing a set of indices indexing an LMat entry
!> @param[in] ct integer between 1 and 3 labeling a position in the index set
!> @return t_dir true if the two indices at the position ct (i.e. ct and ct+step in indices)
!! are the same
pure function is_repeated_pair(indices, ct) result(t_dir)
integer(int64), intent(in) :: indices(num_inds)
integer, intent(in) :: ct
logical :: t_dir
t_dir = (indices(ct) == indices(ct + step))
end function is_repeated_pair
!------------------------------------------------------------------------------------------!
end module LMat_freeze