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
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | iCacheFlag |
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