subroutine CalcUEGMP2()
use SymExcitDataMod, only: kPointToBasisFn
use SystemData, only: ElecPairs, NMAXX, NMAXY, NMAXZ, OrbECutOff, &
tMP2UEGRestrict, kiRestrict, kiMsRestrict, kjRestrict, kjMsRestrict, &
Madelung, tMadelung, tUEGFreeze, FreezeCutoff, kvec, tUEG2
use Determinants, only: GetH0Element4, get_helement_excit
integer :: Ki(3), Kj(3), Ka(3), LowLoop, HighLoop, X, i, Elec1Ind, Elec2Ind, K, Orbi, Orbj
integer :: iSpn, FirstA, nJ(NEl), a_loc, Ex(2, maxExcit), kx, ky, kz, OrbB
integer :: ki2, kj2
logical :: tParity
real(dp) :: Ranger, length
HElement_t(dp) :: hel, H0tmp, mp2, mp2all
!Divvy up the ij pairs
Ranger = real(ElecPairs, dp) / real(nProcessors, dp)
LowLoop = int(iProcIndex * Ranger) + 1
Highloop = int((iProcIndex + 1) * Ranger)
if ((iProcIndex + 1) == nProcessors) Highloop = ElecPairs
if (iProcIndex == 0) then
if (lowLoop /= 1) write(stderr, *) "Error here!"
end if
write(stdout, *) "Total ij pairs: ", ElecPairs
write(stdout, *) "Considering ij pairs from: ", LowLoop, " to ", HighLoop
! write(stdout,*) "HFDet: ",HFDet(:)
do i = LowLoop, HighLoop !Looping over electron pairs on this processor
X = ElecPairs - i
K = INT((SQRT(8.0_dp * REAL(X, dp) + 1.0_dp) - 1.0_dp) / 2.0_dp)
Elec1Ind = NEl - 1 - K
Elec2Ind = NEl - X + ((K * (K + 1)) / 2)
Orbi = HFDet(Elec1Ind)
Orbj = HFDet(Elec2Ind)
Ki = G1(Orbi)%k
Kj = G1(Orbj)%k
!=======================================
if (tUEG2) then
Ki = kvec(Orbi, 1:3)
Kj = kvec(Orbj, 1:3)
end if
!=======================================
if (tUEGFreeze) then
ki2 = ki(1)**2 + ki(2)**2 + ki(3)**2
kj2 = kj(1)**2 + kj(2)**2 + kj(3)**2
if (.not. (ki2 > FreezeCutoff .and. kj2 > FreezeCutoff)) cycle
end if
if (tMP2UEGRestrict) then
if (.not. ( &
(kiRestrict(1) == ki(1) .and. kiRestrict(2) == ki(2) .and. kiRestrict(3) == ki(3) .and. &
kjRestrict(1) == kj(1) .and. kjRestrict(2) == kj(2) .and. kjRestrict(3) == kj(3) .and. &
kjMsRestrict == G1(Orbi)%Ms .and. kiMsRestrict == G1(Orbj)%Ms) .or. &
! the other way round
(kiRestrict(1) == kj(1) .and. kiRestrict(2) == kj(2) .and. kiRestrict(3) == kj(3) .and. &
kjRestrict(1) == ki(1) .and. kjRestrict(2) == ki(2) .and. kjRestrict(3) == ki(3) .and. &
kiMsRestrict == G1(Orbi)%Ms .and. kjMsRestrict == G1(Orbj)%Ms)) &
) cycle
write(stdout, *) "Restricting calculation to i,j pair: ", Orbi, Orbj
end if
IF ((G1(Orbi)%Ms) * (G1(Orbj)%Ms) == -1) THEN
!We have an alpha beta pair of electrons.
iSpn = 2
ELSE
IF (G1(Orbi)%Ms == 1) THEN
!We have an alpha alpha pair of electrons.
iSpn = 3
ELSE
!We have a beta beta pair of electrons.
iSpn = 1
end if
end if
! write(stdout,*) "ijpair: ",Orbi,Orbj
if ((iSpn == 3) .or. (iSpn == 1)) then
if (iSpn == 3) then
FirstA = 2 !Loop over alpha
else
FirstA = 1 !Loop over beta
end if
do a_loc = FirstA, nBasis, 2
!Loop over all a
!Reject if a is occupied
if (IsOcc(iLutHF, a_loc)) cycle
Ka = G1(a_loc)%k
!=======================================
if (tUEG2) then
Ka = kvec(a_loc, 1:3)
end if
!=======================================
!Find k labels of b
kx = Ki(1) + Kj(1) - Ka(1)
if (abs(kx) > NMAXX) cycle
ky = Ki(2) + Kj(2) - Ka(2)
if (abs(ky) > NMAXY) cycle
kz = Ki(3) + Kj(3) - Ka(3)
if (abs(kz) > NMAXZ) cycle
length = real((kx**2) + (ky**2) + (kz**2), dp)
if (length > OrbECutoff) cycle
!Find the actual k orbital
if (iSpn == 3) then
!want alpha
OrbB = kPointToBasisFn(kx, ky, kz, 2)
else
!want beta
OrbB = kPointToBasisFn(kx, ky, kz, 1)
end if
!Reject k orbital if it is occupied or gt a
if (IsOcc(iLutHF, OrbB)) cycle
if (OrbB >= a_loc) cycle
!Find det
call make_double(HFDet, nJ, elec1ind, elec2ind, a_loc, &
orbB, ex, tParity)
!Sum in mp2 contrib
hel = get_helement_excit(HFDet, nJ, 2, Ex, tParity)
H0tmp = getH0Element4(nJ, HFDet)
H0tmp = Fii - H0tmp
if (tMadelung) then
H0tmp = H0tmp + 2.0_dp * Madelung
end if
mp2 = mp2 + (hel**2) / H0tmp
end do
else if (iSpn == 2) then
do a_loc = 1, nBasis
!Loop over all a_loc
!Reject if a is occupied
if (IsOcc(iLutHF, a_loc)) cycle
Ka = G1(a_loc)%k
!=======================================
if (tUEG2) then
Ka = kvec(a_loc, 1:3)
end if
!=======================================
!Find k labels of b
kx = Ki(1) + Kj(1) - Ka(1)
if (abs(kx) > NMAXX) cycle
ky = Ki(2) + Kj(2) - Ka(2)
if (abs(ky) > NMAXY) cycle
kz = Ki(3) + Kj(3) - Ka(3)
if (abs(kz) > NMAXZ) cycle
length = real((kx**2) + (ky**2) + (kz**2), dp)
if (length > OrbECutoff) cycle
!Find the actual k orbital
if (is_beta(a_loc)) then
!want alpha b orbital
OrbB = kPointToBasisFn(kx, ky, kz, 2)
else
!want beta
OrbB = kPointToBasisFn(kx, ky, kz, 1)
end if
!Reject k orbital if it is occupied or gt a
if (IsOcc(iLutHF, OrbB)) cycle
if (OrbB >= a_loc) cycle
!Find det
call make_double(HFDet, nJ, elec1ind, elec2ind, a_loc, &
orbB, ex, tParity)
!Sum in mp2 contrib
hel = get_helement_excit(HFDet, nJ, 2, Ex, tParity)
H0tmp = getH0Element4(nJ, HFDet)
H0tmp = Fii - H0tmp
if (tMadelung) then
H0tmp = H0tmp + 2.0_dp * Madelung
end if
mp2 = mp2 + (hel**2) / H0tmp
end do
end if
end do
mp2all = 0.0_dp
!Sum contributions across nodes.
call MPISumAll(mp2, mp2all)
write(stdout, "(A,2G25.15)") "MP2 energy calculated: ", MP2All, MP2All + Hii
call neci_flush(stdout)
end subroutine CalcUEGMP2