LMat_freeze.F90 Source File


Source Code

Source Code

#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

    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


    !> 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

            ! 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
                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)
                end do
                ! 2 => single excitation
                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)))
                ! 3 => diagonal element
                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


        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
            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
        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
            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
            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
                if(orbs(1) == orbs(2)) prefactor = 2.0_dp
                ! If there is one, two possible spin configs contribute, and we have odd parity
                prefactor = -2.0_dp
                ! 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
            ! 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
            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)
            ! 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

    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