LOGICAL FUNCTION GETCACHEDUMATEL(IDI, IDJ, IDK, IDL, UMATEL, ICACHE, ICACHEI, A, B, ITYPE)
! In:
! IDI,IDJ,IDK,IDL: orbitals indices of the integral(s). These are
! rearranged, as disucssed below.
! Out:
! UMatEl: <ij|u|kl> integral (if found in the cache).
! A,B: Cache indices of (I,K) and (J,L) pairs (see GetCacheIndex).
! ICache,ICacheI: location of integral within the cache (or where it
! should be stored if not found).
! iType: gives information (bit-wise) on where the requested
! integral lies within the cache slot. Useful if the
! integral has to be computed (using the reordered indices).
! Lookup in the cache to see if there's a stored element. If not, return TRUE,
! together with the information on how to compute the integral (re-ordered
! indices, cache indices and a "type" index).
! If the integral exists in the cache, return false and return the stored element in UMatEl.
! NOTE: This will rearrange IDI,IDJ,IDK,IDL into the correct order
! (i,k)<=(j,l) and i<=k, j<=l. ICACHE corresponds to the pair (i,j), and
! ICACHEI is the index in that cache where the cache should be located.
! Note: (i,k)>(j,l) := (k>l) || ((k==l)&&(i>j))
! use SystemData, only : nBasis,G1
INTEGER IDI, IDJ, IDK, IDL, ICACHE, ICACHEI
INTEGER ICACHEI1, ICACHEI2
HElement_t(dp) UMATEL
INTEGER A, B, ITYPE, ISTAR, ISWAP
! LOGICAL tDebug
! IF(IDI.eq.14.and.IDJ.eq.17.and.IDK.eq.23.and.IDL.eq.6) THEN
! write(stdout,*) "Setting tDebug!"
! tDebug=.true.
! ELSE
! tDebug=.false.
! end if
! write(stdout,"(A,4I5)") "GCUI",IDI,IDJ,IDK,IDL
IF (NSLOTS == 0) THEN
!We don't have a cache so signal failure.
GETCACHEDUMATEL = .TRUE.
ICACHE = 0
ITYPE = 0
RETURN
end if
! First ensure the indices are in the correct order
IF (tSmallUMat) THEN
!tSmallUMat is set if we have nStates slots per pair for storing the <ik|jk> integrals.
ITYPE = 0
IF (IDI == IDK) THEN
B = IDI
IF (IDL < IDJ) THEN
CALL SWAP(IDJ, IDL)
ITYPE = 2
end if
CALL GETCACHEINDEX(IDJ, IDL, A)
else if (IDJ == IDL) THEN
B = IDJ
IF (IDK < IDI) THEN
CALL SWAP(IDI, IDK)
ITYPE = 2
end if
CALL GETCACHEINDEX(IDI, IDK, A)
ELSE !Can't find a pair the same
!We don't have a cache for this type of element so signal failure.
GETCACHEDUMATEL = .TRUE.
ICACHE = 0
ITYPE = 0
RETURN
end if
ELSE
!Otherwise Normal caching
ITYPE = 0
ISTAR = 0
ISWAP = 0
IF (IDK < IDI) THEN
CALL SWAP(IDI, IDK)
ISTAR = IOR(ISTAR, 1)
end if
IF (IDL < IDJ) THEN
CALL SWAP(IDJ, IDL)
ISTAR = IOR(ISTAR, 2)
end if
CALL GETCACHEINDEX(IDI, IDK, A)
CALL GETCACHEINDEX(IDJ, IDL, B)
IF (A > B) THEN
CALL SWAP(A, B)
CALL SWAP(IDI, IDJ)
CALL SWAP(IDK, IDL)
ISWAP = 1
end if
IF (HElement_t_size == 1) THEN
! Eight integrals from ijkl are the same.
ITYPE = 0
ELSE
! Complex orbitals and integrals, so we need to consider different types
! Using notation abcd rather than ijkl. 6/2/07 and 19/2/06
!
! <> mean swap pairs <ab|cd> -> <ba|dc>
! *. means complex conjugate of first codensity i.e. <ab|cd> -> <cb|ad>
! .* for second and ** for both.
!
! abcd -> badc <>
! | |
! | \|/ |-> cdab ** -> dcba **<>
! | cbad *. -|
! | |-> bcda *.<>
! \|/
! adcb .* -> dabc .*<>
!Now consider what must occur to the other integral in the slot to recover pair
!(abcd,cbad). 0 indicates 1st in slot, 1 means 2nd. * indicated conjg.
!
! .. abcd cbad 0 1
! *. cbad abcd 1 0
! .* adcb cbad 1* 0*
! ** cdab adcb 0* 1*
! <> badc dabc 0 1*
! *.<> bcda dcba 1* 0
! .*<> dabc badc 1 0*
! **<> dcba bcda 0* 1
! Of the type, bit zero indicates which of the two integrals in a slot to use.
!Bit 1 is set if the integral should be complex conjugated.
! Bit 2 is set if the other integral in the slot should be complex conjugated
!if we are to have the structure (<ij|kl>,<kj|il>) in the slot.
! This is only used by CacheUMatEl when adding a slot to the cache.
!
! <ij|u|kl> and <kj|u|il> (which are distinct when complex orbitals are used).
! TYPE 0 TYPE 1
!
!
IF (ISTAR == 0) THEN
IF (ISWAP == 0) THEN
ITYPE = 0 !0 1
ELSE
ITYPE = 4 !0 1*
end if
else if (ISTAR == 1) THEN
! If we star the first pair, that corresponds to the plain TYPE 1. If we swap too, then we complex conj.
IF (ISWAP == 0) THEN
ITYPE = 1 !1 0
ELSE
ITYPE = 3 !1* 0
end if
else if (ISTAR == 2) THEN
! If we star the second pair, that corresponds to TYPE 1.
! If there's no swap, it's complex conjugated, otherwise it's not.
IF (ISWAP == 0) THEN
ITYPE = 7 !1* 0*
ELSE
ITYPE = 1 !1 0*
end if
else if (ISTAR == 3) THEN
! We've starred both pairs
! We complex conjg setting bit 1 but using type 0
IF (ISWAP == 0) THEN
ITYPE = 6 !0* 1*
ELSE
ITYPE = 2 !0* 1
end if
end if
end if
end if
! IF(tDebug) write(stdout,"(A,8I5)") "GCUE",IDI,IDJ,IDK,IDL,A,B,iType,UMatCacheFlag
ICACHE = A
IF (NSLOTS == NPAIRS .OR. tSmallUMat) THEN
! we've a small enough system to store everything.
ICACHEI = B
ELSE
IF (UMATCACHEFLAG == 1) THEN
! UMATCACHEFLAG=1 means we are storing a sequence of cache elements,
! in a blank cache
! We store them linearly in the cache, and distribute them around later
!Find the last value in the cache
CALL BINARYSEARCH(NPAIRS + 1, UMATLABELS(1:NSLOTS, A), 1, NSLOTS, ICACHEI, ICACHEI1, ICACHEI2)
ICACHEI = ICACHEI1
ICACHEI2 = ICACHEI1
IF (UMatLabels(iCacheI, A) /= 0) iCacheOvCount = iCacheOvCount + 1
!write(stdout,*) "Cache Overwrite", A,B
! write(stdout,*) IDI,IDJ,IDK,IDL
! write(stdout,*) A,B,NSLOTS,NPAIRS
! write(stdout,*) ICACHEI1,ICACHEI2
! write(stdout,*) ICACHEI
ELSE
! IF(tDebug) THEN
! CALL DumpUMatCache(nBasis)
! write(8,*) B,NSLOTS
! write(8,*) UMATLABELS(1:NSLOTS,A)
! end if
CALL BINARYSEARCH(B, UMATLABELS(1:NSLOTS, A), 1, NSLOTS, ICACHEI, ICACHEI1, ICACHEI2)
! IF(tDebug) write(8,*) "***",UMATLABELS(ICACHEI,A),ICACHEI,ICACHEI1,ICACHEI2
end if
end if
IF (UMATLABELS(ICACHEI, ICACHE) == B) THEN
!write(stdout,*) "C",IDI,IDJ,IDK,IDL,ITYPE,UMatCacheData(0:nTypes-1,ICACHEI,ICACHE)
UMATEL = UMatCacheData(IAND(ITYPE, 1), ICACHEI, ICACHE)
#ifdef CMPLX_
IF (BTEST(ITYPE, 1)) UMATEL = CONJG(UMATEL) ! Bit 1 tells us whether we need to complex conjg the integral
#endif
! signal success
GETCACHEDUMATEL = .FALSE.
ELSE
! signal failure
GETCACHEDUMATEL = .TRUE.
! write(stdout8,*) A,B,ICACHEI1,ICACHEI2,ICACHEI
end if
RETURN
END FUNCTION GETCACHEDUMATEL