Subroutine CalcInit() 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