Subroutine SysInit
Use global_utilities
use SymData, only: tAbelian, TwoCycleSymGens, nSymLabels
use constants, only: Pi, Pi2, THIRD
use read_fci
use sym_mod
use SymExcitDataMod, only: kPointToBasisFn
character(*), parameter :: this_routine = 'SysInit'
integer ierr
CHARACTER(len=1) CPAR(3)
CHARACTER(len=3) CPARITY
! For init of mom
TYPE(BasisFN) G
! For the UEG
real(dp) FKF, Rs
! General variables
INTEGER i, j, k, l, iG, k2_max_2d
INTEGER len
type(timer), save :: proc_timer
real(dp) SUM
! Called functions
TYPE(BasisFN) FrzSym
real(dp) dUnscaledE
real(dp), allocatable :: arr_tmp(:, :)
integer, allocatable :: brr_tmp(:)
! UEG2
integer :: AllocateStatus
real(dp), parameter :: EulersConst = 0.5772156649015328606065120900824024_dp
integer, allocatable :: temp_sym_vecs(:, :)
ECORE = 0.0_dp
! //AJWT TBR
! IFDET=0
! TRHOIJND=.false.
proc_timer%timer_name = 'SysInit '
call set_timer(proc_timer)
!C ==-------------------------------------------------------------------==
!C..Input parameters
write(stdout, '(A)') '======== SYSTEM =========='
write(stdout, '(A,I5)') ' NUMBER OF ELECTRONS : ', NEL
! IF(TSPN) THEN
! write(stdout,*) ' Restricting the spin state of the system, TSPN : ' , TSPN
! ELSE
! write(stdout,*) ' No restriction on the spin state of the system, TSPN : ' , TSPN
! end if
if (TSPN) then
write(stdout, *) ' Restricting the S_z spin-projection of the system, TSPN : ', TSPN
write(stdout, '(A,I5)') ' S_z quantum number : ', LMS
else
write(stdout, *) ' No restriction on the S_z spin-projection of the system, TSPN : ', TSPN
end if
NBASISMAX(1:5, 1:7) = 0
TSPINPOLAR = .FALSE.
DO I = 1, 3
SymRestrict%k(I) = IPARITY(I)
end do
SymRestrict%Ms = IPARITY(4)
SymRestrict%Sym%s = IPARITY(5)
IF (TSPN) THEN
!C.. If we're doing monte carlo without generating a list of
!C.. determinants, we cannot as yet force spin or parity, except by
!C.. restricting the basis set. This will only work for Ms=NEL/2
!C IF((.NOT.TBLOCK).AND.(.NOT.TCALCHMAT.OR.NTAY.LE.0)) THEN
!C write(stdout,*) 'TSPN set to TRUE. Determinant list not being',
!C & ' used for MC. Forcing MS=Nel/2'
!C IF(MOD(NEL,2).EQ.0) THEN
!C LMS=NEL/2
!C ELSE
!C LMS=NEL
!C end if
!C TSPINPOLAR=.TRUE.
!C end if
IF (MOD(LMS + NEL * 2, 2) /= MOD(NEL, 2)) THEN
write(stdout, *) 'LMS=', LMS, ' not achievable with', NEL, ' electrons'
write(stdout, *) 'Resetting LMS'
LMS = MOD(NEL, 2)
end if
LMS2 = LMS
ELSE
if (tUEG2) then
TSPN = .TRUE.
if (dimen == 1) then
! 1D calculations are always spin polarised
LMS = NEL
!TSPINPOLAR = .TRUE.
else
LMS = 0
end if
LMS2 = LMS
end if
end if
! global ms should only be outputted if we restrict the spin system
! otherwise not defined ...
if (TSPN) write(stdout, *) ' GLOBAL MS : ', LMS
IF ((NBASISMAX(2, 3) == 1) .or. tROHF) THEN
!If we are dealing with an open shell system, the calculation of symreps will sometimes fail.
!This will have consequences for the rest of the program, so in a slightly hacky way, we can simply
!not reorder the orbitals by energy, so that they remain in symmetries.
!The reason it fails it that it looks for a complete set of orbitals which are degenerate
!and ignores those in the symmetry classification. Unfortunately in ROHF and UHF there aren't such things.
write(stdout, '(A)') " Open shell system - SYMIGNOREENERGIES set. "
! tHFNOORDER=.true.
tSymIgnoreEnergies = .true.
end if
if (tUseBrillouin) THEN
write(stdout, '(A)') " Using Brillouin's Theorem to ignore single excitations"
end if
if (tStoreAsExcitations) THEN
write(stdout, '(A)') " Storing determinants as excitations from the HF determinant. WARNING this may not work!"
IF (nEL < 8) call stop_all(this_routine, ' tStoreAsExcitations requires nEl>=8.')
end if
TwoCycleSymGens = .false.
IF (TCPMD) THEN
write(stdout, '(A)') ' *** GENERIC SYSTEM USING KOHN-SHAM ORBITALS *** '
CALL CPMDSYSTEMINIT(LEN)
IF (TPARITY) THEN
write(stdout, "(A)", advance='no') ' SYMMETRIES : '
CALL WRITEALLSYM(5, SymRestrict, .true.)
end if
IF (THFORDER) write(stdout, '(A)') " Ordering according to 1-electron energies. "
else if (tVASP) THEN
write(stdout, '(A)') ' *** GENERIC SYSTEM USING HF ORBITALS PRODUCED BY VASP *** '
CALL VASPSystemInit(LEN)
else if (TREADINT) THEN
!C.. we read in the integrals from FCIDUMP and ignore most config
!C..
write(stdout, '(A)') ' *** GENERIC SYSTEM *** '
IF (THUB) THEN
THUB = .FALSE.
write(stdout, '(A)') " Setting THUB=.FALSE. "
end if
TwoCycleSymGens = .true.
IF (TDFREAD) THEN
write(stdout, '(A)') " Reading Density fitted integrals. "
LMSBASIS = LMS
CALL InitDFBasis(nBasisMax, Len)
else if (tRIIntegrals) THEN
LMSBASIS = LMS
IF (tRIIntegrals) THEN
write(stdout, '(A)') " Reading RI integrals. "
CALL InitRIBasis(nBasisMax, Len)
LMSBASIS = LMS
ELSE
write(stdout, '(A)') " Reading in all integrals from FCIDUMP file, but storing them in a cache... "
tAbelian = .true.
end if
CALL INITFROMFCID(NEL, NBASISMAX, LEN, LMSBASIS, TBIN)
nBasisMax(3, 3) = 1 !No momentum conservation
nBasisMax(4, 1) = -1 !God knows...
nBasisMax(4, 2) = 1 !Ditto
!.. Correspond to ISS=0
! NBASISMAX(2,3)=-1 !indicate we generate on the fly
iSpinSkip = -1
!NBASISMAX(2,3) !indicate we generate on the fly
!C.. say we're a UHF det so all singles are 0
! IF(LMS.EQ.0) THEN
! NBASISMAX(4,5)=1
! ELSE
! NBASISMAX(4,5)=2
! end if
IF (NBASISMAX(2, 3) == 1) then
write(stdout, '(A)') " Unrestricted calculation. Cave Arthropodia. "
else if (tROHF) then
write(stdout, '(A)') " High-spin restricted calculation. Seriously Cave Arthropodia. "
end if
ELSE
LMSBASIS = LMS
! write(stdout,*) "TBIN:",tBin
CALL INITFROMFCID(NEL, NBASISMAX, LEN, LMSBASIS, TBIN)
!C.. say we're a UHF det so all singles are 0
! IF(LMS.EQ.0) THEN
! NBASISMAX(4,5)=1
! ELSE
! NBASISMAX(4,5)=2
! end if
IF (tStoreSpinOrbs) then
write(stdout, '(A)') " Unrestricted calculation. Cave Arthropodia. "
else if (tROHF) then
write(stdout, '(A)') " High-spin restricted calculation. Seriously Cave Arthropodia. "
end if
end if
ELSE
! ======================================================
if (tUEG2) then
!determines the recip lattice, the real volume, R_s and the Fermi vector
call LatticeInit(RS, FKF)
! engergy cutoff not given -> set to 2* Fermi vector
if (.not. torbEcutoff) then
orbEcutoff = (2.0_dp * FKF / k_lattice_constant)**2.0_dp
torbEcutoff = .true.
end if
! if cell size not given
if (NMAXX == 0 .and. NMAXY == 0 .and. NMAXZ == 0) then
call CalcCell
end if
TTILT = .FALSE.
ALAT = 0.0_dp !shouldn't be used in the UEG2 part...
NMAX = MAX(NMAXX, NMAXY, NMAXZ)
NNR = NMSH * NMSH * NMSH
write(stdout, '(A)') ' *** In UEG2 *** '
write(stdout, '(A)') ' *** UNIFORM ELECTRON GAS CALCULATION *** '
write(stdout, '(A,F20.16)') ' Electron Gas Rs set to ', FUEGRS
IF (TPARITY) THEN
write(stdout, *) ' MOMENTUM : ', (IPARITY(I), I=1, 3)
end if
write(stdout, '(A,I5)') ' Dimension : ', Dimen
write(stdout, *) ' Reciprocal lattice constant : ', k_lattice_constant
write(stdout, '(A,I5)') ' NMAXX : ', NMAXX
write(stdout, '(A,I5)') ' NMAXY : ', NMAXY
write(stdout, '(A,I5)') ' NMAXZ : ', NMAXZ
write(stdout, '(A,I5)') ' NMSH : ', NMSH
write(stdout, *) " Wigner-Seitz radius Rs=", RS
write(stdout, *) " Fermi vector kF^2=", FKF**2
write(stdout, *) " Fermi Energy EF=", FKF * FKF / 2
write(stdout, *) " Unscaled fermi vector kF=", FKF / k_lattice_constant
write(stdout, *) " Unscaled Fermi Energy nmax**2=", (FKF * FKF / 2) / (0.5 * (2 * PI / ALAT(5))**2)
IF (.not. (OrbECutoff.isclose.1e-20_dp)) write(stdout, *) " Orbital Energy Cutoff:", OrbECutoff
write(stdout, '(1X,A,F19.5)') ' VOLUME : ', OMEGA
write(stdout, *) ' TALPHA : ', TALPHA
write(stdout, '(1X,A,F19.5)') ' ALPHA : ', ALPHA
ALPHA = (OMEGA)**THIRD * ALPHA
write(stdout, '(1X,A,F19.5)') ' SCALED ALPHA : ', ALPHA
write(stdout, *) 'Madelung constant: ', Madelung
!C..
!C..Calculate number of basis functions
!C.. UEG allows from -NMAX->NMAX
IF (TSPINPOLAR) THEN
NBASISMAX(4, 1) = 1
!C.. spinskip
NBASISMAX(2, 3) = 1
ELSE
NBASISMAX(4, 1) = -1
!C.. spinskip
! If spinskip is unset
IF (NBASISMAX(2, 3) == 0) NBASISMAX(2, 3) = 2
end if
NBASISMAX(4, 2) = 1
NBASISMAX(1, 1) = -NMAXX
NBASISMAX(1, 2) = NMAXX
NBASISMAX(2, 1) = -NMAXY
NBASISMAX(2, 2) = NMAXY
NBASISMAX(3, 1) = -NMAXZ
NBASISMAX(3, 2) = NMAXZ
NBASISMAX(1, 3) = -1
LEN = (2 * NMAXX + 1) * (2 * NMAXY + 1) * (2 * NMAXZ + 1) * ((NBASISMAX(4, 2) - NBASISMAX(4, 1)) / 2 + 1)
!C.. UEG
NBASISMAX(3, 3) = -1
!C.. we actually store twice as much in arr as we need.
!C.. the ARR(1:LEN,1) are the energies of the orbitals ordered according to
!C.. BRR. ARR(1:LEN,2) are the energies of the orbitals with default
!C.. ordering.
!C.. ARR is reallocated in IntFreezeBasis if orbitals are frozen so that it
!C.. has the correct size and shape to contain the eigenvalues of the active
!C.. basis.
write(stdout, '(A,I5)') " NUMBER OF SPIN ORBITALS IN BASIS : ", Len
allocate(Arr(LEN, 2), STAT=ierr)
LogAlloc(ierr, 'Arr', 2 * LEN, 8, tagArr)
! // TBR
! IP_ARRSTORE=IP_ARR
ARR = 0.0_dp
allocate(Brr(LEN), STAT=ierr)
LogAlloc(ierr, 'Brr', LEN, 4, tagBrr)
BRR(1:LEN) = 0
allocate(G1(Len), STAT=ierr)
LogAlloc(ierr, 'G1', LEN, BasisFNSizeB, tagG1)
G1(1:LEN) = NullBasisFn
IF (TCPMD) THEN
write(stdout, '(A)') '*** INITIALIZING BASIS FNs FROM CPMD ***'
CALL CPMDBASISINIT(NBASISMAX, ARR, BRR, G1, LEN)
NBASIS = LEN
iSpinSkip = NBasisMax(2, 3)
else if (tVASP) THEN
write(stdout, '(A)') '*** INITIALIZING BASIS FNs FROM VASP ***'
CALL VASPBasisInit(ARR, BRR, G1, LEN) ! This also modifies nBasisMax
NBASIS = LEN
iSpinSkip = NBasisMax(2, 3)
else if (TREADINT .AND. TDFREAD) THEN
write(stdout, '(A)') '*** Creating Basis Fns from Dalton output ***'
call InitDaltonBasis(Arr, Brr, G1, Len)
nBasis = Len
call GenMolpSymTable(1, G1, nBasis)
else if (TREADINT) THEN
!This is also called for tRiIntegrals and tCacheFCIDUMPInts
write(stdout, '(A)') '*** CREATING BASIS FNs FROM FCIDUMP ***'
CALL GETFCIBASIS(NBASISMAX, ARR, BRR, G1, LEN, TBIN)
NBASIS = LEN
!C.. we're reading in integrals and have a molpro symmetry table
IF (lNoSymmetry) THEN
write(stdout, *) "Turning Symmetry off"
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
CALL GENMOLPSYMTABLE(1, G1, NBASIS)
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
ELSE
CALL GENMOLPSYMTABLE(NBASISMAX(5, 2) + 1, G1, NBASIS)
end if
ELSE
!C.. Create plane wave basis functions
write(stdout, *) "Creating plane wave basis."
IG = 0
DO I = NBASISMAX(1, 1), NBASISMAX(1, 2)
DO J = NBASISMAX(2, 1), NBASISMAX(2, 2)
DO K = NBASISMAX(3, 1), NBASISMAX(3, 2)
DO L = NBASISMAX(4, 1), NBASISMAX(4, 2), 2
!if (dimen ==1 .AND. L ==1) cycle
G%k(1) = I
G%k(2) = J
G%k(3) = K
G%Ms = L
IF (KALLOWED(G, NBASISMAX)) THEN
CALL GetUEGKE(I, J, K, ALAT, tUEGTrueEnergies, tUEGOffset, k_offset, SUM, dUnscaledE)
IF (dUnscaledE > OrbECutoff) CYCLE
IG = IG + 1
ARR(IG, 1) = SUM
ARR(IG, 2) = SUM
BRR(IG) = IG
!C..These are the quantum numbers: n,l,m and sigma
G1(IG)%K(1) = I
G1(IG)%K(2) = J
G1(IG)%K(3) = K
G1(IG)%MS = L
G1(IG)%Sym = TotSymRep()
end if
end do
end do
end do
end do
if (real_lattice_type == 'sc' .AND. maxval(G1%K(1)) >= NMAXX) write(stdout, *) 'ERROR IN CELL CALCULATION! '
!C..Check to see if all's well
write(stdout, *) ' NUMBER OF BASIS FUNCTIONS : ', IG
NBASIS = IG
! calculate k-vectors in cartesian coordinates
allocate(kvec(NBASIS, 3), STAT=AllocateStatus)
IG = 0
DO I = 1, NBASIS
IG = IG + 1
kvec(IG, 1) = k_lattice_vectors(1, 1) * G1(IG)%K(1) &
+ k_lattice_vectors(2, 1) * G1(IG)%K(2) &
+ k_lattice_vectors(3, 1) * G1(IG)%K(3)
kvec(IG, 2) = k_lattice_vectors(1, 2) * G1(IG)%K(1) &
+ k_lattice_vectors(2, 2) * G1(IG)%K(2) &
+ k_lattice_vectors(3, 2) * G1(IG)%K(3)
kvec(IG, 3) = k_lattice_vectors(1, 3) * G1(IG)%K(1) &
+ k_lattice_vectors(2, 3) * G1(IG)%K(2) &
+ k_lattice_vectors(3, 3) * G1(IG)%K(3)
end do
IF (LEN /= IG) THEN
if (OrbECutoff > -1e20_dp) then
write(stdout, *) " Have removed ", LEN - IG, " high energy orbitals "
! Resize arr and brr.
allocate(arr_tmp(nbasis, 2), brr_tmp(nbasis), stat=ierr)
arr_tmp = arr(1:nbasis, :)
brr_tmp = brr(1:nbasis)
deallocate(arr, brr, stat=ierr)
LogDealloc(tagarr)
LogDealloc(tagbrr)
allocate(arr(nbasis, 2), brr(nbasis), stat=ierr)
LogAlloc(ierr, 'Arr', 2 * nbasis, 8, tagArr)
LogAlloc(ierr, 'Brr', nbasis, 4, tagBrr)
arr = arr_tmp
brr = brr_tmp
deallocate(arr_tmp, brr_tmp, stat=ierr)
else
write(stdout, *) " LEN=", LEN, "IG=", IG
call stop_all(this_routine, ' LEN NE IG ')
end if
end if
if (.not. tHub) CALL GENMOLPSYMTABLE(1, G1, NBASIS)
end if
IF (tFixLz) THEN
write(stdout, '(A)') "****** USING Lz SYMMETRY *******"
write(stdout, '(A,I5)') "Pure spherical harmonics with complex orbitals used to constrain Lz to: ", LzTot
write(stdout, *) "Due to the breaking of the Ml degeneracy, the fock energies are slightly wrong, "&
&//"on order of 1.0e-4_dp - do not use for MP2!"
if (nsymlabels > 4) then
call stop_all(this_routine, "D2h point group detected. Incompatable with Lz symmetry conserving "&
&//"orbitals. Have you transformed these orbitals into spherical harmonics correctly?!")
end if
end if
!C.. (.NOT.TREADINT)
!C.. Set the initial symmetry to be totally symmetric
FrzSym = NullBasisFn
FrzSym%Sym = TotSymRep()
CALL SetupFreezeSym(FrzSym)
!C..Now we sort them using SORT2 and then SORT
!C.. This sorts ARR and BRR into order of ARR [AJWT]
IF (.NOT. THFNOORDER) THEN
CALL ORDERBASIS(NBASIS, ARR, BRR, ORBORDER, NBASISMAX, G1)
ELSE
!.. copy the default ordered energies.
CALL DCOPY(NBASIS, ARR(1, 1), 1, ARR(1, 2), 1)
end if
! write(stdout,*) THFNOORDER, " THFNOORDER"
if (.not. tMolpro) then
!If we are calling from molpro, we write the basis later (after reordering)
CALL writebasis(stdout, G1, nBasis, ARR, BRR)
end if
IF (NEL > NBASIS) call stop_all(this_routine, 'MORE ELECTRONS THAN BASIS FUNCTIONS')
CALL neci_flush(stdout)
NOCC = NEl / 2
IF (TREADINT) THEN
!C.. we're reading in integrals and have a molpro symmetry table
IF (lNoSymmetry) THEN
write(stdout, *) "Turning Symmetry off"
CALL GENMOLPSYMREPS()
ELSE
CALL GENMOLPSYMREPS()
end if
else if (TCPMD) THEN
!C.. If TCPMD, then we've generated the symmetry table earlier,
!C.. but we still need the sym reps table.
call stop_all(this_routine, 'CPMD interface deprecated')
! CALL GENCPMDSYMREPS(G1,NBASIS,ARR,1.e-5_dp)
else if (tVASP) THEN
!C.. If VASP-based calculation, then we've generated the symmetry table earlier,
!C.. but we still need the sym reps table. DEGENTOL=1.0e-6_dp. CHECK w/AJWT.
CALL GENSYMREPS(G1, NBASIS, ARR, 1.e-6_dp)
else if (THUB .AND. .NOT. TREAL) THEN
CALL GenHubMomIrrepsSymTable(G1, nBasis, nBasisMax)
CALL GENHUBSYMREPS(NBASIS, ARR, BRR)
CALL writebasis(stdout, G1, nBasis, ARR, BRR)
ELSE
!C.. no symmetry, so a simple sym table
CALL GENMOLPSYMREPS()
end if
!// TBR
! write(stdout,*) ' ETRIAL : ',ETRIAL
! IF(FCOUL.NE.1.0_dp) write(stdout,*) "WARNING: FCOUL is not 1.0_dp. FCOUL=",FCOUL
IF (FCOULDAMPBETA > 0) write(stdout, *) "FCOUL Damping. Beta ", FCOULDAMPBETA, " Mu ", FCOULDAMPMU
call halt_timer(proc_timer)
!calculate tau if not given
if (TAU < 0.0_dp) then
call CalcTau
end if
if (tMadelung .AND. near_zero(Madelung) .AND. dimen == 3) then
Madelung = calc_madelung()
else if (tMadelung .AND. near_zero(Madelung) .AND. dimen /= 3) then
call stop_all(this_routine, "Calculation of Madelung constant works in 3D only!")
end if
return
end if !UEG2
! ======================================================
IF (TUEG) THEN
write(stdout, '(A)') ' *** UNIFORM ELECTRON GAS CALCULATION *** '
if (iPeriodicDampingType /= 0) then
! We are using either a screened or an attenuated Coulomb
! potential for calculating the exchange integrals.
! This means that we need to be able to distinguish between
! exchange integrals and normal Coulomb integrals and hence we
! should refer to spin-orbitals throughout.
nBasisMax(2, 3) = 1
tStoreSpinOrbs = .true.
end if
if (t_ueg_transcorr) write(stdout, *) 'Using Transcorrelated method on UEG'
IF (.not. near_zero(FUEGRS)) THEN
write(stdout, '(A,F20.16)') ' Electron Gas Rs set to ', FUEGRS
write(stdout, '(A,I4)') ' DIMENSION set to ', dimen
if (dimen == 3) then
OMEGA = BOX * BOX * BOX * BOA * COA
!C.. required density is (3/(4 pi rs^3))
!C.. need omega to be (NEL* 4 pi rs^3 / 3)
!C.. need box to be (NEL*4 pi/(3 BOA COA))^(1/3) rs
BOX = (NEL * 4.0_dp * PI / (3.0_dp * BOA * COA))**(1.0_dp / 3.0_dp)
BOX = BOX * FUEGRS
write(stdout, '(A, F20.16)') " Resetting box size to ", BOX
else if (dimen == 2) then
OMEGA = BOX * BOX * BOA
!C.. required density is (1/( pi rs^2))
!C.. need omega to be (NEL* pi rs^2)
!C.. need box to be (NEL* pi/ BOA)^(1/2) rs
BOX = (NEL * PI / BOA)**(1.0_dp / 2.0_dp)
BOX = BOX * FUEGRS
write(stdout, '(A, F20.16)') " Resetting box size to ", BOX
else
write(stdout, '(A, I4)') " Dimension problem ", dimen
stop
end if
end if
end if
IF (THUB) write(stdout, '(A)') ' *** HUBBARD MODEL *** '
!C..
IF (.NOT. THUB .AND. .NOT. TUEG) THEN
! Just have even/odd symmetry so it's a two cycle symmetry
! generation issue.
TwoCycleSymGens = .true.
write(stdout, '(A)') " Electron in cubic box. "
IF (TPARITY) THEN
write(stdout, '(A)') ' ******************************* '
write(stdout, *) ' PARITY IS ON '
DO I = 1, 3
IF (IPARITY(I) == 1) THEN
CPAR(I) = 'G'
else if (IPARITY(I) == -1) THEN
CPAR(I) = 'U'
ELSE
call stop_all(this_routine, ' !!! PROBLEM WITH PARITY !!! ')
end if
end do
CPARITY = CPAR(1)//CPAR(2)//CPAR(3)
write(stdout, *) ' PARITY : ', CPARITY
ELSE
write(stdout, *) ' PARITY IS OFF '
end if
write(stdout, *) ' ******************************* '
! //TBR
! IF((.NOT.TBLOCK).AND.(.NOT.TCALCHMAT.OR.NTAY.LT.0)) STOP 'CANNOT USE PARITY WITHOUT LIST OF DETS'
ELSE
IF (TPARITY) THEN
write(stdout, *) ' MOMENTUM : ', (IPARITY(I), I=1, 3)
end if
end if
!C..
! W.D: are those variable ever used actually?
NMAX = MAX(NMAXX, NMAXY, NMAXZ)
NNR = NMSH * NMSH * NMSH
write(stdout, '(A,I5)') ' NMAXX : ', NMAXX
write(stdout, '(A,I5)') ' NMAXY : ', NMAXY
write(stdout, '(A,I5)') ' NMAXZ : ', NMAXZ
write(stdout, '(A,I5)') ' NMSH : ', NMSH
!C.. 2D check
IF (NMAXZ == 0) THEN
write(stdout, '(A)') ' NMAXZ=0. 2D calculation using C/A=1/A '
COA = 1 / BOX
end if
!C..
IF (THUB) THEN
write(stdout, '(1X,A,F19.5)') ' HUBBARD T : ', BHUB
write(stdout, '(1X,A,F19.5)') ' HUBBARD U : ', UHUB
if (abs(nn_bhub) > EPS) then
write(stdout, '(1X,A,F19.5)') ' HUBBARD T* : ', nn_bhub
print *, "Also next-nearest neighbor hopping!"
end if
IF (TTILT) write(stdout, *) ' TILTED LATTICE: ', ITILTX, ",", ITILTY
IF (TTILT .AND. ITILTX > ITILTY) call stop_all(this_routine, 'ERROR: ITILTX>ITILTY')
if (t_new_hubbard) then
if (iprocindex == root) then
print *, "New Hubbard Implementation! "
print *, "lattice used: "
call lat%print_lat()
end if
end if
ELSE
write(stdout, '(1X,A,F19.5)') ' BOX LENGTH : ', BOX
write(stdout, '(1X,A,F19.5)') ' B/A : ', BOA
write(stdout, '(1X,A,F19.5)') ' C/A : ', COA
TTILT = .FALSE.
end if
ALAT(1) = BOX
ALAT(2) = BOX * BOA
ALAT(3) = BOX * COA
IF (near_zero(fRc) .AND. iPeriodicDampingType /= 0) THEN
ALAT(4) = BOX * ((BOA * COA) / (4 * PI / 3))**THIRD
ELSE
ALAT(4) = fRc
end if
IF (THUB) THEN
if (t_new_hubbard) then
omega = real(lat%get_nsites(), dp)
if (iprocindex == root) then
print *, " periodic boundary conditions: ", lat%is_periodic()
print *, "Real space basis: ", t_new_real_space_hubbard
end if
else if (t_heisenberg_model .or. t_tJ_model) then
root_print "New tJ/Heisenberg Implementation! "
root_print "lattice used: "
if (iprocindex == root) then
call lat%print_lat()
end if
omega = real(lat%get_nsites(), dp)
root_print " periodic boundary conditions: ", lat%is_periodic()
rs = 1.0_dp
else
write(stdout, *) ' X-LENGTH OF HUBBARD CHAIN:', NMAXX
write(stdout, *) ' Y-LENGTH OF HUBBARD CHAIN:', NMAXY
write(stdout, *) ' Z-LENGTH OF HUBBARD CHAIN:', NMAXZ
write(stdout, *) ' Periodic Boundary Conditions:', TPBC
write(stdout, *) ' Real space basis:', TREAL
IF (TTILT .AND. THUB) THEN
OMEGA = real(NMAXX, dp) * NMAXY * (ITILTX * ITILTX + ITILTY * ITILTY)
ELSE
OMEGA = real(NMAXX, dp) * (NMAXY) * (NMAXZ)
end if
end if
RS = 1.0_dp
ELSE
if (dimen == 3) then
OMEGA = ALAT(1) * ALAT(2) * ALAT(3)
RS = (3.0_dp * OMEGA / (4.0_dp * PI * NEL))**THIRD
ALAT(5) = RS
IF (iPeriodicDampingType /= 0) THEN
IF (iPeriodicDampingType == 1) THEN
write(stdout, *) " Using attenuated Coulomb potential for exchange interactions."
else if (iPeriodicDampingType == 2) THEN
write(stdout, *) " Using cut-off Coulomb potential for exchange interactions."
end if
write(stdout, *) " Rc cutoff: ", ALAT(4)
end if
write(stdout, *) " Wigner-Seitz radius Rs=", RS
FKF = (9 * PI / 4)**THIRD / RS
write(stdout, *) " Fermi vector kF=", FKF
write(stdout, *) " Fermi Energy EF=", FKF * FKF / 2
write(stdout, *) " Unscaled Fermi Energy nmax**2=", (FKF * FKF / 2) / (0.5 * (2 * PI / ALAT(5))**2)
else if (dimen == 2) then
OMEGA = ALAT(1) * ALAT(2)
RS = dsqrt(OMEGA / (PI * NEL))
ALAT(5) = RS
IF (iPeriodicDampingType /= 0) THEN
IF (iPeriodicDampingType == 1) THEN
write(stdout, *) " Using attenuated Coulomb potential for exchange interactions."
else if (iPeriodicDampingType == 2) THEN
write(stdout, *) " Using cut-off Coulomb potential for exchange interactions."
end if
write(stdout, *) " Rc cutoff: ", ALAT(4)
end if
write(stdout, *) " Wigner-Seitz radius Rs=", RS
FKF = dsqrt(2.0_dp) / RS
write(stdout, *) " Fermi vector kF=", FKF
write(stdout, *) " Fermi Energy EF=", FKF * FKF / 2
write(stdout, *) " Unscaled Fermi Energy nmax**2=", (FKF * FKF / 2) / (0.5 * (2 * PI / ALAT(5))**2) !??????????
else if (dimen == 1) then
OMEGA = ALAT(1)
RS = OMEGA / (2 * NEL)
ALAT(5) = RS
write(stdout, *) 'Be cautios, the 1D rs and kF values have not been checked thorougHly!'
IF (iPeriodicDampingType /= 0) THEN
IF (iPeriodicDampingType == 1) THEN
write(stdout, *) " Using attenuated Coulomb potential for exchange interactions."
else if (iPeriodicDampingType == 2) THEN
write(stdout, *) " Using cut-off Coulomb potential for exchange interactions."
end if
write(stdout, *) " Rc cutoff: ", ALAT(4)
end if
write(stdout, *) " Wigner-Seitz radius Rs=", RS
FKF = PI * RS
write(stdout, *) " Fermi vector kF=", FKF
write(stdout, *) " Fermi Energy EF=", FKF * FKF / 2
write(stdout, *) " Unscaled Fermi Energy nmax**2=", (FKF * FKF / 2) / (0.5 * (2 * PI / ALAT(5))**2)
else
write(stdout, *) 'Dimension problem'
stop
end if
end if
IF (.not. (OrbECutoff.isclose.1e-20_dp)) write(stdout, *) " Orbital Energy Cutoff:", OrbECutoff
write(stdout, '(1X,A,F19.5)') ' VOLUME : ', OMEGA
write(stdout, *) ' TALPHA : ', TALPHA
write(stdout, '(1X,A,F19.5)') ' ALPHA : ', ALPHA
ALPHA = MIN(ALAT(1), ALAT(2), ALAT(3)) * ALPHA
write(stdout, '(1X,A,F19.5)') ' SCALED ALPHA : ', ALPHA
!C..
!C..Calculate number of basis functions
!C.. UEG allows from -NMAX->NMAX
IF (TSPINPOLAR) THEN
NBASISMAX(4, 1) = 1
!C.. spinskip
NBASISMAX(2, 3) = 1
ELSE
NBASISMAX(4, 1) = -1
!C.. spinskip
! If spinskip is unset
IF (NBASISMAX(2, 3) == 0) NBASISMAX(2, 3) = 2
end if
NBASISMAX(4, 2) = 1
IF (THUB) THEN
if (t_new_hubbard) then
! i need a new setup routine for this for the new generic
! hubbard setup! essentialy this just sets up NBASISMAX ..
! what do i need from this legacy variable??
len = 2 * lat%get_nsites()
if (t_k_space_hubbard) then
! this indicates pbc and k-space.
! especially for the addelecsym function!
NBASISMAX(1, 3) = 0
else if (t_new_real_space_hubbard) then
NBASISMAX(1, 3) = 4
NBASISMAX(3, 3) = 0
end if
else if (t_heisenberg_model .or. t_tJ_model) then
len = 2 * lat%get_nsites()
nBasisMax(1, 3) = 4
nBasisMax(3, 3) = 0
else
IF (TTILT) THEN
CALL SETBASISLIM_HUBTILT(NBASISMAX, NMAXX, NMAXY, NMAXZ, LEN, TPBC, ITILTX, ITILTY)
! is supported now!
ELSE
CALL SETBASISLIM_HUB(NBASISMAX, NMAXX, NMAXY, NMAXZ, LEN, TPBC, TREAL)
end if
end if
else if (TUEG) THEN
NBASISMAX(1, 1) = -NMAXX
NBASISMAX(1, 2) = NMAXX
NBASISMAX(2, 1) = -NMAXY
NBASISMAX(2, 2) = NMAXY
NBASISMAX(3, 1) = -NMAXZ
NBASISMAX(3, 2) = NMAXZ
NBASISMAX(1, 3) = -1
if (dimen == 3) then
LEN = (2 * NMAXX + 1) * (2 * NMAXY + 1) * (2 * NMAXZ + 1) * ((NBASISMAX(4, 2) - NBASISMAX(4, 1)) / 2 + 1)
else if (dimen == 2) then
LEN = (2 * NMAXX + 1) * (2 * NMAXY + 1) * ((NBASISMAX(4, 2) - NBASISMAX(4, 1)) / 2 + 1)
else if (dimen == 1) then
LEN = (2 * NMAXX + 1) * ((NBASISMAX(4, 2) - NBASISMAX(4, 1)) / 2 + 1)
else
write(stdout, '(A, I4)') " Dimension problem ", dimen
stop
end if
!C.. UEG
NBASISMAX(3, 3) = -1
ELSE
NBASISMAX(1, 1) = 1
NBASISMAX(1, 2) = NMAXX
NBASISMAX(2, 1) = 1
NBASISMAX(2, 2) = NMAXY
NBASISMAX(3, 1) = 1
NBASISMAX(3, 2) = NMAXZ
NBASISMAX(1, 3) = 0
LEN = NMAXX * NMAXY * NMAXZ * ((NBASISMAX(4, 2) - NBASISMAX(4, 1)) / 2 + 1)
NBASISMAX(1, 3) = 0
!C.. particle in box
NBASISMAX(3, 3) = -2
tAbelian = .true.
end if
end if
!C.. (.NOT.TREADINT)
!C.. we actually store twice as much in arr as we need.
!C.. the ARR(1:LEN,1) are the energies of the orbitals ordered according to
!C.. BRR. ARR(1:LEN,2) are the energies of the orbitals with default
!C.. ordering.
!C.. ARR is reallocated in IntFreezeBasis if orbitals are frozen so that it
!C.. has the correct size and shape to contain the eigenvalues of the active
!C.. basis.
write(stdout, '(A,I5)') " NUMBER OF SPIN ORBITALS IN BASIS : ", Len
allocate(Arr(LEN, 2), STAT=ierr)
LogAlloc(ierr, 'Arr', 2 * LEN, 8, tagArr)
! // TBR
! IP_ARRSTORE=IP_ARR
ARR = 0.0_dp
allocate(Brr(LEN), STAT=ierr)
LogAlloc(ierr, 'Brr', LEN, 4, tagBrr)
BRR(1:LEN) = 0
allocate(G1(Len), STAT=ierr)
LogAlloc(ierr, 'G1', LEN, BasisFNSizeB, tagG1)
G1(1:LEN) = NullBasisFn
IF (TCPMD) THEN
write(stdout, '(A)') '*** INITIALIZING BASIS FNs FROM CPMD ***'
CALL CPMDBASISINIT(NBASISMAX, ARR, BRR, G1, LEN)
NBASIS = LEN
iSpinSkip = NBasisMax(2, 3)
else if (tVASP) THEN
write(stdout, '(A)') '*** INITIALIZING BASIS FNs FROM VASP ***'
CALL VASPBasisInit(ARR, BRR, G1, LEN) ! This also modifies nBasisMax
NBASIS = LEN
iSpinSkip = NBasisMax(2, 3)
else if (TREADINT .AND. TDFREAD) THEN
write(stdout, '(A)') '*** Creating Basis Fns from Dalton output ***'
call InitDaltonBasis(Arr, Brr, G1, Len)
nBasis = Len
call GenMolpSymTable(1, G1, nBasis)
else if (TREADINT) THEN
!This is also called for tRiIntegrals and tCacheFCIDUMPInts
write(stdout, '(A)') '*** CREATING BASIS FNs FROM FCIDUMP ***'
CALL GETFCIBASIS(NBASISMAX, ARR, BRR, G1, LEN, TBIN)
NBASIS = LEN
!C.. we're reading in integrals and have a molpro symmetry table
IF (lNoSymmetry) THEN
write(stdout, *) "Turning Symmetry off"
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
CALL GENMOLPSYMTABLE(1, G1, NBASIS)
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
ELSE
CALL GENMOLPSYMTABLE(NBASISMAX(5, 2) + 1, G1, NBASIS)
end if
ELSE
!C.. Create plane wave basis functions
if (treal) then
write(stdout, *) "Creating real-space basis."
else
write(stdout, *) "Creating plane wave basis."
end if
if (t_new_hubbard) then
BRR = [(i, i=1, 2 * lat%get_nsites())]
IG = 2 * lat%get_nsites()
if (t_new_real_space_hubbard) then
! i have to do everything what is done below with my new
! lattice class:
! something like: (loop over spin-orbital!)
ARR = 0.0_dp
else
ARR(:, 1) = [(bhub * lat%dispersion_rel_spin_orb(i), i=1, IG)]
ARR(:, 2) = [(bhub * lat%dispersion_rel_spin_orb(i), i=1, IG)]
end if
do i = 1, lat%get_nsites()
! todo: not sure if i really should set up the k-vectors
! for the real-space! maybe the convention actually is to
! set them all to 0! check that!
! do i still want to save the "real-space positions" in the
! k-vectors? i could.. but do i need it?
G1(2 * i - 1)%k = lat%get_k_vec(i)
G1(2 * i - 1)%ms = -1
G1(2 * i - 1)%Sym = TotSymRep()
G1(2 * i)%k = lat%get_k_vec(i)
G1(2 * i)%ms = 1
G1(2 * i)%Sym = TotSymRep()
! and in the k-space i still need to
end do
else if (t_heisenberg_model .or. t_tJ_model) then
BRR = [(i, i=1, 2 * lat%get_nsites())]
IG = 2 * lat%get_nsites()
ARR = 0.0_dp
do i = 1, lat%get_nsites()
! todo: not sure if i really should set up the k-vectors
! for the real-space! maybe the convention actually is to
! set them all to 0! check that!
! do i still want to save the "real-space positions" in the
! k-vectors? i could.. but do i need it?
G1(2 * i - 1)%k = lat%get_k_vec(i)
G1(2 * i - 1)%ms = -1
G1(2 * i - 1)%Sym = TotSymRep()
G1(2 * i)%k = lat%get_k_vec(i)
G1(2 * i)%ms = 1
G1(2 * i)%Sym = TotSymRep()
! and in the k-space i still need to
end do
else if (dimen == 3) then
IG = 0
DO I = NBASISMAX(1, 1), NBASISMAX(1, 2)
DO J = NBASISMAX(2, 1), NBASISMAX(2, 2)
DO K = NBASISMAX(3, 1), NBASISMAX(3, 2)
DO L = NBASISMAX(4, 1), NBASISMAX(4, 2), 2
G%k(1) = I
G%k(2) = J
G%k(3) = K
G%Ms = L
! i have to change this condition to accomodate both real
! and k-space hubbard models for the tilted and cubic
! lattice
! aperiodic does not work anyway and is not really needed..
! if tilted i want to check if k is allowed otherwise
if ((treal .and. .not. ttilt) .or. kAllowed(G, NBASISMAX)) then
IF (THUB) THEN
!C..Note for the Hubbard model, the t is defined by ALAT(1)!
call setupMomIndexTable()
call setupBreathingCont(2 * btHub / OMEGA)
IF (TPBC) THEN
CALL HUBKIN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
ELSE
CALL HUBKINN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
end if
else if (TUEG) THEN
CALL GetUEGKE(I, J, K, ALAT, tUEGTrueEnergies, tUEGOffset, k_offset, SUM, dUnscaledE)
! Bug for non-trueEnergy
! IF(dUnscaledE.gt.OrbECutoff) CYCLE
IF (tUEGTrueEnergies) then
IF (dUnscaledE > OrbECutoff) CYCLE
ELSE
IF (SUM > OrbECutoff) CYCLE
END IF
ELSE
SUM = (BOX**2) * ((I * I / ALAT(1)**2) + (J * J / ALAT(2)**2) + (K * K / ALAT(3)**2))
end if
IF (.NOT. TUEG .AND. SUM > OrbECutoff) CYCLE
IG = IG + 1
ARR(IG, 1) = SUM
ARR(IG, 2) = SUM
BRR(IG) = IG
!C..These are the quantum numbers: n,l,m and sigma
G1(IG)%K(1) = I
G1(IG)%K(2) = J
G1(IG)%K(3) = K
G1(IG)%MS = L
G1(IG)%Sym = TotSymRep()
end if
end do
end do
end do
end do
else if (dimen == 2) then
IG = 0
k2_max_2d = int(OrbECutoff) + 1
do k = 0, k2_max_2d
DO I = NBASISMAX(1, 1), NBASISMAX(1, 2)
DO J = NBASISMAX(2, 1), NBASISMAX(2, 2)
DO L = NBASISMAX(4, 1), NBASISMAX(4, 2), 2
G%k(1) = I
G%k(2) = J
G%k(3) = 0
G%Ms = L
! change to implement the tilted real-space
if ((treal .and. .not. ttilt) .or. KALLOWED(G, nBasisMax)) then
! IF((THUB.AND.(TREAL.OR..NOT.TPBC)).OR.KALLOWED(G,NBASISMAX)) THEN
IF (THUB) THEN
!C..Note for the Hubbard model, the t is defined by ALAT(1)!
IF (TPBC) THEN
CALL HUBKIN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
ELSE
CALL HUBKINN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
end if
else if (TUEG) THEN
CALL GetUEGKE(I, J, 0, ALAT, tUEGTrueEnergies, tUEGOffset, k_offset, SUM, dUnscaledE)
! Bug for non-trueEnergy
! IF(dUnscaledE.gt.OrbECutoff) CYCLE
IF (tUEGTrueEnergies) then
IF (dUnscaledE > OrbECutoff) CYCLE
ELSE
if ((sum > (k * 1.d0 + 1.d-20)) .or. (sum < (k * 1.d0 - 1.d-20))) CYCLE
IF (SUM > OrbECutoff) CYCLE
END IF
ELSE
SUM = (BOX**2) * ((I * I / ALAT(1)**2) + (J * J / ALAT(2)**2))
end if
IF (.NOT. TUEG .AND. SUM > OrbECutoff) CYCLE
IG = IG + 1
ARR(IG, 1) = SUM
ARR(IG, 2) = SUM
BRR(IG) = IG
!C..These are the quantum numbers: n,l,m and sigma
G1(IG)%K(1) = I
G1(IG)%K(2) = J
G1(IG)%K(3) = 0
G1(IG)%MS = L
G1(IG)%Sym = TotSymRep()
end if
end do
end do
end do
end do
!C..Check to see if all's well
else if (dimen == 1) then
IG = 0
DO I = NBASISMAX(1, 1), NBASISMAX(1, 2)
DO J = NBASISMAX(2, 1), NBASISMAX(2, 2)
DO K = NBASISMAX(3, 1), NBASISMAX(3, 2)
DO L = NBASISMAX(4, 1), NBASISMAX(4, 2), 2
G%k(1) = I
G%k(2) = J
G%k(3) = K
G%Ms = L
! change to implement the tilted real-space
if ((treal .and. .not. ttilt) .or. KALLOWED(G, nBasisMax)) then
! IF((THUB.AND.(TREAL.OR..NOT.TPBC)).OR.KALLOWED(G,NBASISMAX))
! THEN
IF (THUB) THEN
!C..Note for the Hubbard model, the t is defined by ALAT(1)!
call setupMomIndexTable()
call setupBreathingCont(2 * btHub / OMEGA)
IF (TPBC) THEN
CALL HUBKIN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
ELSE
CALL HUBKINN(I, J, K, NBASISMAX, BHUB, TTILT, SUM, TREAL)
end if
else if (TUEG) THEN
CALL GetUEGKE(I, J, K, ALAT, tUEGTrueEnergies, tUEGOffset, k_offset, SUM, dUnscaledE)
! Bug for non-trueEnergy
! IF(dUnscaledE.gt.OrbECutoff) CYCLE
IF (tUEGTrueEnergies) then
IF (dUnscaledE > OrbECutoff) CYCLE
ELSE
IF (SUM > OrbECutoff) CYCLE
END IF
ELSE
SUM = (BOX**2) * ((I * I / ALAT(1)**2) + (J * J / ALAT(2)**2) + (K * K / ALAT(3)**2))
end if
IF (.NOT. TUEG .AND. SUM > OrbECutoff) CYCLE
IG = IG + 1
ARR(IG, 1) = SUM
ARR(IG, 2) = SUM
BRR(IG) = IG
!C..These are the quantum numbers: n,l,m and sigma
G1(IG)%K(1) = I
G1(IG)%K(2) = J
G1(IG)%K(3) = K
G1(IG)%MS = L
G1(IG)%Sym = TotSymRep()
end if
end do
end do
end do
end do
else
write(stdout, '(A, I4)') " Dimension problem ", dimen
stop
end if
!C..Check to see if all's well
write(stdout, *) ' NUMBER OF BASIS FUNCTIONS : ', IG
NBASIS = IG
! try to order all of the stuff in ascending orbital number..
! but this could mean a lot of changes in other parts of the code..
! first i would have to sort..
! what about reordering the basis functions and the associated
! k-points in ascending single particle energy order?
! for the guga implementation this would be beneficiary for
! the efficiency of the code..
! for now implement it only for the GUGA case to not mess to
! much with the rest of the code..
if (tGUGA .and. (.not. t_guga_noreorder)) then
! i have to sort the alpha and betas seperately, due the
! possible degeneracies
call sort(arr(:, 1), brr, nskip=2)
call sort(arr(2:nBasis, 1), brr(2:nBasis), nskip=2)
call sort(arr(:, 2))
! i have to make a copy of the G1 array i guess to temporarily
! save the symmetry information..
allocate(temp_sym_vecs(nBasis, 4))
do i = 1, nBasis
temp_sym_vecs(i, 1) = G1(i)%k(1)
temp_sym_vecs(i, 2) = G1(i)%k(2)
temp_sym_vecs(i, 3) = G1(i)%k(3)
temp_sym_vecs(i, 4) = G1(i)%ms
end do
! and now restore the information in the G1
do i = 1, nBasis
G1(i)%k(1) = temp_sym_vecs(brr(i), 1)
G1(i)%k(2) = temp_sym_vecs(brr(i), 2)
G1(i)%k(3) = temp_sym_vecs(brr(i), 3)
G1(i)%ms = temp_sym_vecs(brr(i), 4)
! now it is okay to sort brr ascending or?
brr(i) = i
end do
! hm.. and what about degeneracies and the connected
! k-vectors.. is it enough to just leave them as they are,
! since they are degenerate anyway... in the "normal"
! implementation it also does not seem like there is any
! change in the association with the k-vectors, although
! the brr and arr arrays are sorted again..
end if
IF (LEN /= IG) THEN
IF (OrbECutoff > -1e20_dp) then
write(stdout, *) " Have removed ", LEN - IG, " high energy orbitals "
! Resize arr and brr.
allocate(arr_tmp(nbasis, 2), brr_tmp(nbasis), stat=ierr)
arr_tmp = arr(1:nbasis, :)
brr_tmp = brr(1:nbasis)
deallocate(arr, brr, stat=ierr)
LogDealloc(tagarr)
LogDealloc(tagbrr)
allocate(arr(nbasis, 2), brr(nbasis), stat=ierr)
LogAlloc(ierr, 'Arr', 2 * nbasis, 8, tagArr)
LogAlloc(ierr, 'Brr', nbasis, 4, tagBrr)
arr = arr_tmp
brr = brr_tmp
deallocate(arr_tmp, brr_tmp, stat=ierr)
else
write(stdout, *) " LEN=", LEN, "IG=", IG
call stop_all(this_routine, ' LEN NE IG ')
end if
end if
! also turn off symmetries for the real-space basis
if (thub .and. treal) then
write(stdout, *) "Turning Symmetry off for the real-space Hubbard"
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
call GENMOLPSYMTABLE(1, G1, NBASIS)
DO I = 1, nBasis
G1(I)%Sym%s = 0
end do
end if
if (.not. tHub) CALL GENMOLPSYMTABLE(1, G1, NBASIS)
end if
IF (tFixLz) THEN
write(stdout, '(A)') "****** USING Lz SYMMETRY *******"
write(stdout, '(A,I5)') "Pure spherical harmonics with complex orbitals used to constrain Lz to: ", LzTot
write(stdout, *) "Due to the breaking of the Ml degeneracy, the fock energies are slightly wrong, "&
&//"on order of 1.0e-4_dp - do not use for MP2!"
if (nsymlabels > 4) then
call stop_all(this_routine, "D2h point group detected. Incompatable with Lz symmetry conserving "&
&//"orbitals. Have you transformed these orbitals into spherical harmonics correctly?!")
end if
end if
!C.. (.NOT.TREADINT)
!C.. Set the initial symmetry to be totally symmetric
FrzSym = NullBasisFn
FrzSym%Sym = TotSymRep()
CALL SetupFreezeSym(FrzSym)
!C..Now we sort them using SORT2 and then SORT
!C.. This sorts ARR and BRR into order of ARR [AJWT]
IF (.NOT. THFNOORDER) THEN
CALL ORDERBASIS(NBASIS, ARR, BRR, ORBORDER, NBASISMAX, G1)
ELSE
!.. copy the default ordered energies.
CALL DCOPY(NBASIS, ARR(1, 1), 1, ARR(1, 2), 1)
end if
! write(stdout,*) THFNOORDER, " THFNOORDER"
if (.not. tMolpro) then
!If we are calling from molpro, we write the basis later (after reordering)
CALL writebasis(stdout, G1, nBasis, ARR, BRR)
end if
IF (NEL > NBASIS) &
call stop_all(this_routine, 'MORE ELECTRONS THAN BASIS FUNCTIONS')
CALL neci_flush(stdout)
IF (TREAL .AND. THUB) THEN
!C.. we need to allow integrals between different spins
NBASISMAX(2, 3) = 1
end if
NOCC = NEl / 2
IF (TREADINT) THEN
!C.. we're reading in integrals and have a molpro symmetry table
IF (lNoSymmetry) THEN
write(stdout, *) "Turning Symmetry off"
CALL GENMOLPSYMREPS()
ELSE
CALL GENMOLPSYMREPS()
end if
else if (TCPMD) THEN
!C.. If TCPMD, then we've generated the symmetry table earlier,
!C.. but we still need the sym reps table.
CALL GENCPMDSYMREPS(G1, NBASIS, ARR)
else if (tVASP) THEN
!C.. If VASP-based calculation, then we've generated the symmetry table earlier,
!C.. but we still need the sym reps table. DEGENTOL=1.0e-6_dp. CHECK w/AJWT.
CALL GENSYMREPS(G1, NBASIS, ARR, 1.e-6_dp)
else if (THUB .AND. .NOT. TREAL) THEN
! [W.D. 25.1.2018:]
! this also has to be changed for the new hubbard implementation
if (t_k_space_hubbard) then
! or i change the function below to account for the new
! implementation
call setup_symmetry_table()
! call gen_symreps()
else
CALL GenHubMomIrrepsSymTable(G1, nBasis, nBasisMax)
end if
! this function does not make sense..
CALL GENHUBSYMREPS(NBASIS, ARR, BRR)
CALL writebasis(stdout, G1, nBasis, ARR, BRR)
ELSE
!C.. no symmetry, so a simple sym table
CALL GENMOLPSYMREPS()
end if
!C..
!// TBR
! write(stdout,*) ' TREAD : ' , TREAD
!// TBR
! write(stdout,*) ' ETRIAL : ',ETRIAL
! IF(FCOUL.NE.1.0_dp) write(stdout,*) "WARNING: FCOUL is not 1.0_dp. FCOUL=",FCOUL
IF (FCOULDAMPBETA > 0) write(stdout, *) "FCOUL Damping. Beta ", FCOULDAMPBETA, " Mu ", FCOULDAMPMU
! ====================== GUGA implementation ===========================
! design decision to have as much guga related functionality stored in
! guga_*.F90 files and only call necessary routines.
! routine guga_init() is found in module guga_init.F90
! it prints out info and sets the nReps to be the number of orbitals
! have to call routine at this late stage in the initialisation,
! because it needs info of the number of orbitals from the integral
! input files..: lets hope FDET or something similar is not getting
! allocated before this point...
! CHANGE: switched init functions to guga_data
call halt_timer(proc_timer)
End Subroutine SysInit