IntInit Subroutine

public subroutine IntInit(iCacheFlag)

Uses

C.. Real space hubbard C.. we pre-compute the 2-e integrals C.. Generate the 2e integrals (UMAT) C.. we pre-compute the 2-e integrals C.. Generate the 2e integrals (UMAT) C.. Non-periodic hubbard (mom space) C.. Most normal Hubbards C.. The UEG doesn’t store coul integrals C.. We need to init the arrays regardless of whether we’re storing H C..Need to initialise the Fourier arrays C..

C.. We’re doing a 2D simulation C.. we pre-compute the 2-e integrals C.. Generate the 2e integrals (UMAT) C.. we need to generate TMAT - Now setup in individual routines

C..Cube multiplier C.. the UEG has k=2pi n/L rather than pi n/L, so we need to multiply the C.. KE up by 4

Arguments

Type IntentOptional Attributes Name
integer :: iCacheFlag

Source Code

    Subroutine IntInit(iCacheFlag)
!who knows what for
        Use global_utilities
        Use OneEInts, only: SetupTMat, SetupPropInts, OneEPropInts, PropCore, SetupFieldInts, OneEFieldInts, FieldCore, UpdateOneEInts
        USE UMatCache, only: FreezeTransfer, CreateInvBRR, GetUMatSize, SetupUMat2D_df
        Use UMatCache, only: SetupUMatCache
        use SystemData, only: nBasisMax, Alpha, BHub, BRR, nmsh, nEl
        use SystemData, only: G1, iSpinSkip, nBasis, nMax, nMaxZ
        use SystemData, only: Omega, tAlpha, TBIN, tCPMD, tDFread, tHFOrder
        use SystemData, only: thub, tpbc, treadint, ttilt, TUEG, tPickVirtUniform
        use SystemData, only: uhub, arr, alat, treal, tReltvy
        use SymExcitDataMod, only: tBuildOccVirtList, tBuildSpinSepLists
        use LoggingData, only: tCalcPropEst, iNumPropToEst, EstPropFile
        use CalcData, only: nFields_it, tCalcWithField, FieldFiles_it
        use Parallel_neci, only: iProcIndex, MPIBcast
        use MemoryManager, only: TagIntType
        use sym_mod, only: GenSymStatePairs
        use read_fci
        use LoggingData, only: t_umat_output
        use real_space_hubbard, only: init_tmat
        use k_space_hubbard, only: init_tmat_kspace
        use lattice_mod, only: lat

        INTEGER iCacheFlag
        complex(dp), ALLOCATABLE :: ZIA(:)
        INTEGER(TagIntType), SAVE :: tagZIA = 0
        INTEGER i!,j,k,l,idi,idj,idk,idl,Index1
        INTEGER TmatInt
        integer(int64) :: UMatInt, ii
        real(dp) :: UMatMem
        integer iErr, IntSize
        character(25), parameter :: this_routine = 'IntInit'

        FREEZETRANSFER = .false.

        IF(THFBASIS) THEN
            write(stdout, *) "Using Hartree-Fock Basis"
            IF(.NOT. THFCALC) write(stdout, *) "Reading Hartree-Fock Basis"
        ENDIF

!I_VMAX.EQ.1.AND..NOT.TENERGY.AND.
        IF(.not. tNeedsVirts .and. NTFROZEN == 0) THEN
            write(stdout, *) "MaxVertices=1 and NTFROZEN=0."
            IF(ABS(ARR(NEL, 1) - ARR(NEL + 1, 1)) > 1.0e-3_dp) THEN
                write(stdout, *) "NEL spinorbitals completely fills up degenerate set."
                write(stdout, *) "Only calculating vertex sum, so freezing all virtuals."
                NTFROZEN = NBASIS - NEL
            ELSE
                write(stdout, *) "NEL spinorbitals does not completely fill up degenerate set."
                write(stdout, *) "NOT freezing all virtuals."
            end if
        end if

!     IF(.NOT.(TREAD.OR.TREADRHO)) THEN
        IF(NTFROZEN < 0) THEN
            I = NEL + (-NTFROZEN)
        ELSE
            I = nBasis - NTFROZEN
        ENDIF
        ! Initialize the umat-index function
        call setup_UMatInd()
        IF(TCPMD) THEN
!.. We don't need to do init any 4-index integrals, but we do need to init the 2-index
            write(stdout, *) " *** INITIALIZING CPMD 2-index integrals ***"
            call shared_allocate_mpi(umat_win, umat,(/1_int64/))
            !allocate(UMat(1), stat=ierr)
            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
            CALL GENSymStatePairs(nBasis / 2, .false.)
            CALL SetupTMAT(nBasis, 2, TMATINT)
            CALL CPMDINIT2INDINT(nBasis, I, G1, NEL, ECORE, tHFOrder, ARR, BRR, iCacheFlag)
        else if(tVASP) THEN
            call shared_allocate_mpi(umat_win, umat,(/1_int64/))
            !allocate(UMat(1), stat=ierr)
            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
            CALL GENSymStatePairs(nBasis / 2, .false.)
            CALL SetupTMAT(nBasis, 2, TMATINT)
            CALL VASPInitIntegrals(I, ECore, tHFOrder)
!      else if(TREADINT.AND.TReadCacheInts)...
!set up dummy umat, read in TMAT, <ij|ij> and <ii|jj>
!Work out size needed for cache
!Initialise cache
!read in integral and put in cache
!change flag to read integrals from cache
        else if(TREADINT .AND. TDFREAD) THEN
            call shared_allocate_mpi(umat_win, umat,(/1_int64/))
            !allocate(UMat(1), stat=ierr)
            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
            CALL SetupTMAT(nBasis, 2, TMATINT)
            Call ReadDalton1EIntegrals(G1, nBasis, ECore)
            Call ReadDF2EIntegrals(nBasis, I)
        else if(TREADINT .AND. tRIIntegrals) THEN
            call shared_allocate_mpi(umat_win, umat,(/1_int64/))
            !allocate(UMat(1), stat=ierr)
            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
!Why is this called twice here?!
            CALL SetupTMAT(nBasis, 2, TMATINT)
            CALL SetupTMAT(nBasis, iSpinSkip, TMATINT)
            Call ReadRIIntegrals(nBasis, I)
            CALL READFCIINT(UMAT, umat_win, NBASIS, ECORE)
            NBASISMAX(2, 3) = 0
            write(stdout, *) ' ECORE=', ECORE
        else if(TREADINT) THEN
            write(stdout, '(A)') '*** READING PRIMITIVE INTEGRALS FROM FCIDUMP ***'
!.. Generate the 2e integrals (UMAT)
            ISPINSKIP = NBasisMax(2, 3)
            IF(ISPINSKIP <= 0) call stop_all(this_routine, 'NBASISMAX(2,3) ISpinSkip unset')
!nBasisMax(2,3) is iSpinSkip = 1 if UHF and 2 if RHF/ROHF
            CALL GetUMatSize(nBasis, UMATINT)
            write(stdout, *) "UMatSize: ", UMATINT
            UMatMem = REAL(UMatInt, dp) * REAL(HElement_t_sizeB, dp) * (9.536743164e-7_dp)
            write(stdout, "(A,G20.10,A)") "Memory required for integral storage: ", UMatMem, " Mb/Shared Memory"
            call neci_flush(stdout)
            call shared_allocate_mpi(umat_win, umat,(/UMatInt/))
            !allocate(UMat(UMatInt), stat=ierr)
            LogAlloc(ierr, 'UMat', int(UMatInt), HElement_t_SizeB, tagUMat)
            if(iprocindex == 0) then
!for very large UMats, the intrinic zeroing can cause a crash. In that case do an explicit zeroing
                if(UMatInt <= 1000000000) then
                    UMat = 0.0_dp
                else
                    do ii = 1, UMatInt
                        UMat(ii) = 0.0_dp
                    end do
                end if
            end if
!nBasisMax(2,3) is iSpinSkip = 1 if UHF and 2 if RHF/ROHF
            CALL SetupTMAT(nBasis, iSpinSkip, TMATINT)
            IF(TBIN) THEN
                CALL READFCIINTBIN(UMAT, ECORE)
            ELSE
                CALL READFCIINT(UMAT, umat_win, NBASIS, ECORE)
            end if
            write(stdout, *) 'ECORE=', ECORE
            if(tCalcWithField) then
                call SetupFieldInts(nBasis,nFields_it)
                do i=1,nFields_it
                    call ReadPropInts(nBasis,FieldFiles_it(i),FieldCore(i),OneEFieldInts(:,:,i))
                    WRITE(stdout,*) 'FieldCORE(', trim(FieldFiles_it(i)), ')=',FieldCore(i)
                    call MPIBCast(FieldCore(i),1)
                    IntSize = nBasis*nBasis
                    call MPIBCast(OneEFieldInts(:,:,i),IntSize)
                end do
                call UpdateOneEInts(nBasis,nFields_it)
            endif

            IF(tCalcPropEst) THEN
                call SetupPropInts(nBasis)
                do i = 1, iNumPropToEst
                    call ReadPropInts(nBasis, EstPropFile(i), PropCore(i), OneEPropInts(:, :, i))
                    WRITE(stdout,*) 'PropCORE(', trim(EstPropFile(i)), ')=',PropCore(i)
                    call MPIBCast(PropCore(i), 1)
                    IntSize = nBasis * nBasis
                    call MPIBCast(OneEPropInts(:, :, i), IntSize)
                end do
            end if
        ELSE
            ISPINSKIP = NBASISMAX(2, 3)
            IF(NBASISMAX(1, 3) >= 0) THEN
                IF(TUEG .OR. THUB) THEN
                    IF(THUB .AND. TREAL) THEN
                        !!C.. Real space hubbard
                        !!C.. we pre-compute the 2-e integrals
                        if((t_new_real_space_hubbard .and. .not. t_trans_corr_hop) &
                               .or. t_tJ_model .or. t_heisenberg_model) then
                            write(stdout, *) "Not precomputing HUBBARD 2-e integrals"
                            UMatInt = 1_int64
                            call shared_allocate_mpi(umat_win, umat,(/UMatInt/))
                            !allocate(UMat(1), stat=ierr)
                            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
                            UMAT(1) = UHUB
                        else if(t_new_real_space_hubbard .and. t_trans_corr_hop) then

                            call init_hopping_transcorr()
                            call init_umat_rs_hub_transcorr()
                        else
                            write(stdout, *) "Generating 2e integrals"
                            !!C.. Generate the 2e integrals (UMAT)
                            CALL GetUMatSize(nBasis, UMATINT)
                            call shared_allocate_mpi(umat_win, umat,(/UMatInt/))
                            !allocate(UMat(UMatInt), stat=ierr)
                            LogAlloc(ierr, 'UMat', int(UMatInt), HElement_t_SizeB, tagUMat)
                            UMat = 0.0_dp
                            write(stdout, *) "Size of UMat is: ", UMATINT
                            CALL CALCUMATHUBREAL(NBASIS, UHUB, UMAT)
                        end if

                    else if(THUB .AND. .NOT. TPBC) THEN
                        !!C.. we pre-compute the 2-e integrals
                        write(stdout, *) "Generating 2e integrals"
                        !!C.. Generate the 2e integrals (UMAT)
                        CALL GetUMatSize(nBasis, UMATINT)
                        call shared_allocate_mpi(umat_win, umat,(/UMatInt/))
                        !allocate(UMat(UMatInt), stat=ierr)
                        LogAlloc(ierr, 'UMat', int(UMatInt), HElement_t_SizeB, tagUMat)
                        UMat = 0.0_dp
                        !!C.. Non-periodic hubbard (mom space)
                        call gen_coul_hubnpbc
                    ELSE
                        !!C.. Most normal Hubbards
                        IF(.NOT. TUEG) THEN
                            ISPINSKIP = -1
                            NBASISMAX(2, 3) = -1
                            write(stdout, *) "Not precomputing HUBBARD 2-e integrals"
                            call shared_allocate_mpi(umat_win, umat,(/1_int64/))
                            !allocate(UMat(1), stat=ierr)
                            LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
                            UMAT(1) = UHUB / OMEGA
                            ! try to output the UMat anyway. to apply dmrg comparison
                            if(t_umat_output) call write_kspace_umat()

                        end if
                        !!C.. The UEG doesn't store coul integrals
                    end if
                ELSE
                    !!C.. We need to init the arrays regardless of whether we're storing H
                    !!C..Need to initialise the Fourier arrays
                    allocate(Fck(nMsh**3), stat=ierr)
                    LogAlloc(ierr, 'FCK', NMSH**3, 16, tagFCK)
                    allocate(ZIA(2 * (NMSH + 1) * NMAX * NMAX))
                    LogAlloc(ierr, 'ZIA', 2 * (NMSH + 1) * NMAX * NMAX, 16, tagZIA)
                    write(stdout, *) NMSH, NMAX
                    !!C..
!               CALL N_MEMORY_CHECK()
                    IF(NMAXZ == 0) THEN
                        !!C..  We're doing a 2D simulation
                        CALL INITFOU2D(NMSH, FCK, NMAX, ALAT, TALPHA, ALPHA, OMEGA, ZIA)
                    ELSE
                        CALL INITFOU(NMSH, FCK, NMAX, ALAT, TALPHA, ALPHA, OMEGA, ZIA)
                    end if
!               CALL N_MEMORY_CHECK()
                    !!C.. we pre-compute the 2-e integrals
                    write(stdout, *) "Generating 2e integrals"
                    !!C.. Generate the 2e integrals (UMAT)
                    CALL GetUMatSize(nBasis, UMATINT)
                    call shared_allocate_mpi(umat_win, umat,(/UMatInt/))
                    !allocate(UMat(UMatInt), stat=ierr)
                    LogAlloc(ierr, 'UMat', int(UMatInt), HElement_t_SizeB, tagUMat)
                    UMat = 0.0_dp
#ifndef CMPLX_
                    CALL GEN_COUL(NBASISMAX, nBasis, G1, NMSH, NMAX, FCK, UMAT, ZIA)
#else
                    ! From my (Oskar) limited understanding of GEN_COUL,
                    !  I don't think it makes sense for
                    !  complex code.
                    call stop_all(this_routine, 'Should not be here')
#endif
                    deallocate(ZIA)
                    LogDealloc(tagZIA)
                    LogDealloc(tagFCK)
                    Deallocate(FCK)
                end if
            ELSE
                write(stdout, *) "Not precomputing 2-e integrals"
                call shared_allocate_mpi(umat_win, umat,(/1_int64/))
                !allocate(UMat(1), stat=ierr)
                LogAlloc(ierr, 'UMat', 1, HElement_t_SizeB, tagUMat)
            end if
!         CALL N_MEMORY_CHECK()
            !!C.. we need to generate TMAT - Now setup in individual routines
            !CALL N_MEMORY(IP_TMAT,HElement_t_size*nBasis*nBasis,'TMAT')
            !TMAT=(0.0_dp)
            IF(THUB) THEN
                if(t_new_hubbard) then
                    if(t_k_space_hubbard) then
                        ! also change here to the new k-space implementation
                        call init_tmat_kspace(lat)

                    else if(t_new_real_space_hubbard) then
                        call init_tmat(lat)
                    end if
                else
                    CALL CALCTMATHUB(NBASIS, NBASISMAX, BHUB, TTILT, G1, TREAL, TPBC)
                end if
            ELSE
                !!C..Cube multiplier
                CST = PI * PI / (2.0_dp * ALAT(1) * ALAT(1))
                !!C.. the UEG has k=2pi n/L rather than pi n/L, so we need to multiply the
                !!C.. KE up by 4
                IF(NBASISMAX(1, 1) <= 0) CST = CST * 4
                CALL CALCTMATUEG(NBASIS, ALAT, G1, CST, NBASISMAX(1, 1) <= 0, OMEGA)
            end if
        end if

        if((tPickVirtUniform .and. tBuildSpinSepLists .and. tBuildOccVirtList) .and. .not. tReltvy) then
            call stop_all(this_routine, "pick-virt-uniform-mag option needs TREL=.TRUE. in FCIDUMP")
        end if

!      write(stdout,*) "ONE ELECTRON"
!      do i=1,nBasis
!          do j=1,nBasis
!              write(36,"(2I5,G25.10)") i,j,GetTMatEl(i,j)
!          end do
!      end do
!      write(stdout,*) "TWO ELECTRON"
!      do i=1,nBasis
!          do j=1,nBasis
!              do k=1,nBasis
!                  do l=1,nBasis
!                     IDi = GTID(i)
!                     IDj = GTID(j)
!                     IDk = GTID(k)
!                     IDl = GTID(l)
!                     Index1=UMatInd(idi,idj,idk,idl)
!                     write(37,"(9I5,G25.10)") i,j,k,l,
!idi,idj,idk,idl,Index1,GetUMatEl(NBasisMax,UMAT,ALAT,nBasis,ISpinSkip,G1,idi,idj,idk,idl)
!                 end do
!             end do
!         end do
!     end do

        ! Setup the umatel pointers as well
        call init_getumatel_fn_pointers()
    End Subroutine IntInit