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