calc_rdm_spin Subroutine

public subroutine calc_rdm_spin(rdm, rdm_norm, rdm_spin)


Type IntentOptional Attributes Name
type(rdm_list_t), intent(in) :: rdm
real(kind=dp), intent(in) :: rdm_norm(rdm%sign_length)
real(kind=dp), intent(out) :: rdm_spin(rdm%sign_length)


Source Code

Source Code

    subroutine calc_rdm_spin(rdm, rdm_norm, rdm_spin)

        ! Return the (unnormalised) estimate of <S^2> from the instantaneous
        ! 2-RDM estimates.

        use rdm_data, only: rdm_spawn_t
        use SystemData, only: nel
        use UMatCache, only: spatial

        type(rdm_list_t), intent(in) :: rdm
        real(dp), intent(in) :: rdm_norm(rdm%sign_length)
        real(dp), intent(out) :: rdm_spin(rdm%sign_length)

        integer(int_rdm) :: ijkl
        integer :: ielem
        integer :: ij, kl, i, j, k, l ! spin orbitals
        integer :: p, q, r, s ! spatial orbitals
        real(dp) :: rdm_sign(rdm%sign_length)

        rdm_spin = 0.0_dp

        ! Loop over all RDM elements.
        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)
            ! Obtain spatial orbital labels.
            p = spatial(i); q = spatial(j);
            r = spatial(k); s = spatial(l);
            ! Note to the reader for the following code: if mod(i,2) == 1 then
            ! i is a beta (b) orbital, if mod(i,2) == 0 then it is an alpha (a)
            ! obrital.

            ! The following if-statement allows IJIJ spatial combinations.
            if (p == r .and. q == s) then
                ! If we get to this point then we definitely have a contribution
                ! to add in, so extract the sign.
                call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

                ! If all labels have the same spatial part (IIII):
                if (p == q) then
                    if (is_beta(i) .and. is_alpha(j) .and. is_beta(k) .and. is_alpha(l)) then
                        rdm_spin = rdm_spin - 6.0_dp * rdm_sign
                    end if

                    ! We only get here if the spatial parts obey IJIJ, for I /= J:

                    ! The following if-statement allows the following spin combinations:
                    ! aaaa, bbbb, abab and baba.
                    if (mod(i, 2) == mod(k, 2) .and. mod(j, 2) == mod(l, 2)) then

                        if (mod(i, 2) == mod(j, 2)) then
                            ! aaaa and bbbb.
                            rdm_spin = rdm_spin + 2.0_dp * rdm_sign
                            ! abab and baba.
                            rdm_spin = rdm_spin - 2.0_dp * rdm_sign
                        end if
                        ! We only get here if the spin parts are abba or baab.
                        rdm_spin = rdm_spin + 4.0_dp * rdm_sign
                    end if

                end if
            end if
        end do

        rdm_spin = rdm_spin + 3.0_dp * real(nel, dp) * rdm_norm
        rdm_spin = rdm_spin / 4.0_dp

    end subroutine calc_rdm_spin