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 .or. tZeroDelimitedFCIDUMP)) 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