Subroutine CalcInit()
use constants, only: dp
use SystemData, only: G1, Alat, Beta, BRR, ECore, LMS, nBasis, nBasisMax, STot, nMsh, nEl, tSmallBasisForThreeBody
use SystemData, only: tUEG, nOccAlpha, nOccBeta, ElecPairs, tExactSizeSpace, tMCSizeSpace, MaxABPairs, tMCSizeTruncSpace
use SystemData, only: tContact
use IntegralsData, only: FCK, CST, nMax, UMat
use IntegralsData, only: HFEDelta, HFMix, NHFIt, tHFCalc
Use DetCalc, only: DetInv, nDet, tRead
Use DetCalcData, only: ICILevel
use hilbert_space_size, only: FindSymSizeofSpace, FindSymSizeofTruncSpace
use hilbert_space_size, only: FindSymMCSizeofSpace, FindSymMCSizeExcitLevel
use global_utilities
use sltcnd_mod, only: initSltCndPtr, sltcnd_0_base, sltcnd_0_tc
use excitation_types, only: Excite_0_t
real(dp) CalcT2, GetRhoEps
INTEGER I, IC, J, norb
INTEGER nList
HElement_t(dp) HDiagTemp, h_2_temp, h_3_temp
type(Excite_0_t) :: NoExc
character(*), parameter :: this_routine = 'CalcInit'
!Checking whether we have large enoguh basis for ultracold atoms and
!three-body excitations
if(tContact .and. ((nBasis / 2) < (noccAlpha + 2) .or. (nBasis / 2) < (noccBeta + 2))) then
if(noccAlpha == 1 .or. noccBeta == 1) then
tSmallBasisForThreeBody = .false.
else
write(stdout, *) 'There is not enough unoccupied orbitals for a poper three-body ', &
'excitation! Some of the three-body excitations are possible', &
'some of or not. If you really would like to calculate this system, ', &
'you have to implement the handling of cases, which are not possible.'
stop
end if
else
tSmallBasisForThreeBody = .true.
end if
! initialize the slater condon rules
call initSltCndPtr()
allocate(MCDet(nEl))
call LogMemAlloc('MCDet', nEl, 4, this_routine, tagMCDet)
IF(NPATHS == -1) THEN
write(stdout, *) 'NPATHS=-1. SETTING NPATHS to NDET'
NPATHS = NDET
end if
IF(NDET > 0 .AND. ABS(DETINV) + NPATHS > NDET) THEN
write(stdout, *) 'DETINV+NPATHS=', ABS(DETINV) + NPATHS, '>NDET=', NDET
write(stdout, *) 'Setting DETINV and NPATHS to 0'
DETINV = 0
NPATHS = 0
end if
IF(THFCALC) THEN
write(stdout, *) "Calculating Hartree-Fock Basis"
write(stdout, *) "Max Iterations:", NHFIT
write(stdout, *) "FMIX,EDELTA", HFMIX, HFEDELTA
end if
IF(TMONTE) THEN
write(stdout, *) 'MC Determinant Symmetry:'
write(stdout, *)(MDK(I), I=1, 4)
end if
! Thus would appear to be obsolete
! IF(G_VMC_FAC.LE.0) THEN
! write(stdout,*) "G_VMC_FAC=",G_VMC_FAC
! call stop_all(this_routine, "G_VNC_FAC LE 0")
! end if
IF(.not. near_zero(BETAP)) THEN
I_P = NINT(BETA / BETAP)
IF(.not. tFCIMC) THEN
write(stdout, *) 'BETAP=', BETAP
write(stdout, *) 'RESETTING P '
IF(I_P > 100000) write(stdout, *) '*** WARNING I_P=', I_P
end if
end if
IF(.not. tFCIMC) write(stdout, *) 'BETA, P :', BETA, I_P
!C DBRAT=0.001
!C DBETA=DBRAT*BETA
! actually i have to initialize the matrix elements here
if(t_lattice_model) then
if(t_tJ_model) then
if(tGUGA) then
call init_get_helement_tj_guga()
else
call init_get_helement_tj()
end if
else if(t_heisenberg_model) then
if(tGUGA) then
call init_get_helement_heisenberg_guga()
else
call init_get_helement_heisenberg()
end if
else if(t_new_real_space_hubbard) then
call init_get_helement_hubbard()
else if(t_k_space_hubbard) then
call init_get_helement_k_space_hub
end if
end if
IF(.NOT. TREAD) THEN
IC = 0
HDiagTemp = get_helement(fDet, fDet, 0)
write(stdout, *) '<D0|H|D0>=', real(HDiagTemp, dp)
write(stdout, *) '<D0|T|D0>=', CALCT(FDET, NEL)
if (t_3_body_excits) then
if (t_k_space_hubbard) then
h_2_temp = get_2_body_diag_transcorr(fdet)
h_3_temp = get_3_body_diag_transcorr(fdet)
write(stdout, *) "<D0|U|D0>=", h_2_temp
write(stdout, *) "<D0|L|D0>=", h_3_temp
else
HDiagTemp = sltcnd_0_tc(fdet, NoExc)
h_2_temp = sltcnd_0_base(fdet, NoExc) - calct(fdet, nel)
h_3_temp = HDiagTemp - sltcnd_0_base(fdet, NoExc)
write(stdout, *) "<D0|U|D0>=", h_2_temp
write(stdout, *) "<D0|L|D0>=", h_3_temp
end if
else
write(stdout, *) "<D0|U|D0>=", real(HDiagTemp,dp) - calct(fdet, nel)
end if
IF(TUEG) THEN
! The actual KE rather than the one-electron part of the Hamiltonian
write(stdout, *) 'Kinetic=', CALCT2(FDET, NEL, G1, ALAT, CST)
end if
end if
! Find out the number of alpha and beta electrons. For restricted calculations, these should be the same.
! TODO: in GUGA this information is not quite correct, since there
! is no notion of alpha/beta orbitals only positively or negatively
! coupled orbitals, with respect to the total spin quantum number
! but for now, leave it in to not break the remaining code, which
! esp. in the excitation generator depends on those values!
! But change this in future and include a corresponding CalcInitGUGA()
if(tGUGA) then
write(stdout, *) " !! NOTE: running a GUGA simulation, so following info makes no sense!"
write(stdout, *) " but is kept for now to not break remaining code!"
end if
nOccAlpha = 0
nOccBeta = 0
do i = 1, NEl
j = fdet(i)
ic = 0
IF(G1(J)%Ms == 1) THEN
! Orbital is an alpha orbital
nOccAlpha = nOccAlpha + 1
ELSE
nOccBeta = nOccBeta + 1
end if
end do
write(stdout, "(A,I5,A,I5,A)") " FDet has ", nOccAlpha, " alpha electrons, and ", nOccBeta, " beta electrons."
ElecPairs = (NEl * (NEl - 1)) / 2
MaxABPairs = (nBasis * (nBasis - 1) / 2)
! And stats on the number of different types of electron pairs
! that can be found
AA_elec_pairs = nOccAlpha * (nOccAlpha - 1) / 2
BB_elec_pairs = nOccBeta * (nOccBeta - 1) / 2
par_elec_pairs = AA_elec_pairs + BB_elec_pairs
AB_elec_pairs = nOccAlpha * nOccBeta
if(AA_elec_pairs + BB_elec_pairs + AB_elec_pairs /= ElecPairs) &
call stop_all(this_routine, "Calculation of electron pairs failed")
write(stdout, *) ' ', AA_elec_pairs, &
' alpha-alpha occupied electron pairs'
write(stdout, *) ' ', BB_elec_pairs, &
' beta-beta occupied electron pairs'
write(stdout, *) ' ', AB_elec_pairs, &
' alpha-beta occupied electron pairs'
! Get some stats about available numbers of holes, etc.
ASSERT(.not. btest(nbasis, 0))
norb = nbasis / 2
nholes = nbasis - nel
nholes_a = norb - nOccAlpha
nholes_b = norb - nOccBeta
! And count the available hole pairs!
hole_pairs = nholes * (nholes - 1) / 2
AA_hole_pairs = nholes_a * (nholes_a - 1) / 2
BB_hole_pairs = nholes_b * (nholes_b - 1) / 2
AB_hole_pairs = nholes_a * nholes_b
par_hole_pairs = AA_hole_pairs + BB_hole_pairs
if(par_hole_pairs + AB_hole_pairs /= hole_pairs) &
call stop_all(this_routine, "Calculation of hole pairs failed")
write(stdout, *) ' ', AA_hole_pairs, 'alpha-alpha (vacant) hole pairs'
write(stdout, *) ' ', BB_hole_pairs, 'beta-beta (vacant) hole pairs'
write(stdout, *) ' ', AB_hole_pairs, 'alpha-beta (vacant) hole pairs'
IF(tExactSizeSpace) THEN
IF(ICILevel == 0) THEN
CALL FindSymSizeofSpace(6)
ELSE
CALL FindSymSizeofTruncSpace(6)
end if
end if
IF(tMCSizeSpace) THEN
CALL FindSymMCSizeofSpace(6)
end if
if(tMCSizeTruncSpace) then
CALL FindSymMCSizeExcitLevel(6)
end if
IF(TMCDET) THEN
!C.. Generate the determinant from which we start the MC
NLIST = 1
CALL GENSYMDETSS(get_basisfn(MDK), NEL, G1, BRR, NBASIS, MCDET, NLIST, NBASISMAX)
IF(NLIST == 0) THEN
!C.. we couldn't find a det of that symmetry
call stop_all(this_routine, 'Cannot find MC start determinant of correct symmetry')
end if
ELSE
DO I = 1, NEL
MCDET(I) = FDET(I)
end do
end if
IF(TMONTE) THEN
write(stdout, "(A)", advance='no') 'MC Start Det: '
call write_det(stdout, mcDet, .true.)
end if
!C.. we need to calculate a value for RHOEPS, so we approximate that
!C.. RHO_II~=exp(-BETA*H_II/p). RHOEPS is a %ge of this
!C.. we have put TMAT instead of ZIA
IF(I_HMAX /= -20) THEN
!C.. If we're using rhos,
RHOEPS = GETRHOEPS(RHOEPSILON, BETA, NEL, BRR, I_P)
! write(stdout,*) "RHOEPS:",RHOEPS
ELSE
!C.. we're acutally diagonalizing H's, so we just leave RHOEPS as RHOEPSILON
RHOEPS = RHOEPSILON
end if
End Subroutine CalcInit