#include "macros.h" MODULE System use SystemData use CalcData, only: tTruncInitiator, InitiatorWalkNo, & occCASorbs, virtCASorbs, tPairedReplicas, & S2Init, tDynamicAvMcEx use FciMCData, only: tGenMatHEl, t_initialized_roi, t_adjoint_replicas use sort_mod use SymExcitDataMod, only: tBuildOccVirtList, tBuildSpinSepLists use constants use read_fci, only: FCIDUMP_name, load_orb_perm use util_mod, only: error_function, error_function_c, & near_zero, operator(.isclose.), get_free_unit, operator(.div.), & stop_all, neci_flush use lattice_mod, only: lattice, lat use k_space_hubbard, only: setup_symmetry_table use breathing_Hub, only: setupMomIndexTable, setupBreathingCont use tc_three_body_data, only: LMatEps, tSparseLMat use SD_spin_purification_mod, only: possible_purification_methods, & SD_spin_purification, spin_pure_J use MPI_wrapper, only: iprocindex, root use fcimcdata, only: pParallel use tc_three_body_data, only: LMatEps, tSparseLMat use guga_data, only: tGUGACore use fortran_strings, only: to_upper, to_lower, to_int, to_int64, to_realdp use tau_main, only: tau, assign_value_to_tau use gasci, only: GAS_specification, GAS_exc_gen, possible_GAS_exc_gen, & user_input_GAS_exc_gen, CumulGASSpec_t, LocalGASSpec_t, FlexibleGASSpec_t use gasci_util, only: t_output_GAS_sizes use gasci_pchb_main, only: GAS_PCHB_user_input, GAS_PCHB_user_input_vals use pchb_excitgen, only: FCI_PCHB_user_input, FCI_PCHB_user_input_vals use cpmdinit_mod, only: CPMDBASISINIT, GENCPMDSYMREPS, cpmdsysteminit use growing_buffers, only: buffer_int_1D_t use hubbard_mod, only: genhubmomirrepssymtable, genhubsymreps, & hubkin, hubkinn, setbasislim_hubtilt, setbasislim_hub use Determinants, only: getuegke, orderbasis, writebasis IMPLICIT NONE contains subroutine SetSysDefaults() != Set defaults for Calc data items. use default_sets USE SymData, only: tAbelianFastExcitGen USE SymData, only: tStoreStateList use OneEInts, only: tOneElecDiag t_calc_adjoint = .false. ! Default from SymExcitDataMod tBuildOccVirtList = .false. ! SYSTEM defaults - leave these as the default defaults ! Any further addition of defaults should change these after via ! specifying a new set of DEFAULTS. ! implementation of spin adapted GUGA approach t_twisted_bc = .false. t_consider_diff_bias = .false. tReltvy = .false. tComplexOrbs_RealInts = .false. tComplexWalkers_RealInts = .false. tReadFreeFormat = .true. tMolproMimic = .false. tNoSingExcits = .false. tOneElecDiag = .false. tMCSizeTruncSpace = .false. iMCCalcTruncLev = 0 tOddS_HPHF = .false. tRotatedOrbsReal = .false. !This is set if compiled in real, but reading in a complex FCIDUMP. tISKFuncs = .false. !This is for kpoint symmetry with inversion so that determinants can be combined. tKPntSym = .false. !This is for k-point symmetry with the symrandexcit2 excitation generators. tNoSinglesPossible = .false. tMCSizeSpace = .false. t_impurity_excitgen = .false. CalcDetPrint = 1000 CalcDetCycles = 10000 tFixLz = .false. tStoreSpinOrbs = .false. !by default we store/lookup integrals as spatial integrals tNoBrillouin = .true. tBrillouinsDefault = .true. tROHF = .false. tHPHFInts = .false. tHPHF = .false. tMaxHLGap = .false. UMatEps = 1.0e-8 LMatEps = 1.0e-10 tExactSizeSpace = .false. iRanLuxLev = 3 !This is the default level of quality for the random number generator. tNoSymGenRandExcits = .false. tNonUniRandExcits = .true. tCycleOrbs = .false. TSTARBIN = .false. TREADINT = .false. THFORDER = .false. TDFREAD = .false. tRIIntegrals = .false. TPBC = .false. TCPMD = .false. tVASP = .false. THUB = .false. TUEG = .false. tUEG2 = .false. tHeisenberg = .false. tLatticeGens = .false. tNoFailAb = .false. LMS = 0 TSPN = .false. STOT = 0 TPARITY = .false. IParity(:) = 0 dimen = 3 NMAXX = 0 NMAXY = 0 NMAXZ = 0 NMSH = 32 BOX = 1.0_dp BOA = 1.0_dp COA = 1.0_dp TUSEBRILLOUIN = .false. tUHF = .false. FUEGRS = 0.0_dp iPeriodicDampingType = 0 fRc = 0.0_dp TEXCH = .true. UHUB = 4 BHUB = -1 btHub = 0.0_dp TREAL = .false. tUEGTrueEnergies = .false. tUEGSpecifyMomentum = .false. tUEGOffset = .false. TTILT = .false. TALPHA = .false. ALPHA = 0.0_dp ISTATE = 1 OrbECutoff = 1e20_dp tOrbECutoff = .false. gCutoff = 1e20_dp ! This shouldn't be used tgCutoff = .false. tMadelung = .false. Madelung = 0.0_dp tUEGFreeze = .false. FreezeCutoff = 1e20_dp ! This shouldn't be used tMP2UEGRestrict = .false. tStoreAsExcitations = .false. TBIN = .false. tAbelianFastExcitGen = .true. TNoRenormRandExcits = .false. tStoreStateList = .false. !This will be turned to true by default if not in abelian symmetry tAssumeSizeExcitgen = .false. tLagrange = .false. tShake = .false. lNoSymmetry = .false. tRotateOrbs = .false. TimeStep = 0.01_dp ConvergedForce = 0.001_dp ShakeConverged = 0.001_dp tShakeApprox = .false. tShakeDelay = .false. ShakeStart = 1 tROIteration = .false. ROIterMax = 5000 tShakeIter = .false. ShakeIterMax = 5 tSeparateOccVirt = .false. tOffDiagSqrdMin = .false. tOffDiagSqrdMax = .false. tOffDiagMin = .false. tOffDiagMax = .false. tDoubExcMin = .false. tHijSqrdMin = .false. tOneElIntMax = .false. tOnePartOrbEnMax = .false. tHFSingDoubExcMax = .false. tERLocalization = .false. tVirtCoulombMax = .false. tVirtExchangeMin = .false. tRotateOccOnly = .false. tRotateVirtOnly = .false. tSpinOrbs = .false. tReadInCoeff = .false. tUseMP2VarDenMat = .false. tUseHFOrbs = .false. tFindCINatOrbs = .false. DiagWeight = 1.0_dp OffDiagWeight = 1.0_dp OneElWeight = 1.0_dp OrbEnMaxAlpha = 1.0_dp MaxMinFac = 1 DiagMaxMinFac = -1 OneElMaxMinFac = 1 tDiagonalizehij = .false. tHFNoOrder = .false. tSymIgnoreEnergies = .false. tPickVirtUniform = .false. tStoquastize = .false. tAllSymSectors = .false. tGenHelWeighted = .false. tGen_4ind_weighted = .false. tGen_4ind_part_exact = .false. tGen_4ind_lin_exact = .false. tGen_4ind_reverse = .false. tUEGNewGenerator = .false. tGen_4ind_2 = .false. tGen_4ind_2_symmetric = .false. tmodHub = .false. t_uniform_excits = .false. t_mol_3_body = .false. t_ueg_3_body = .false. t_ueg_transcorr = .false. t_trcorr_gausscutoff = .false. t_ueg_dump = .false. t_exclude_3_body_excits = .false. t_pcpp_excitgen = .false. t_fci_pchb_excitgen = .false. ! use weighted singles for the pchb excitgen? t_guga_pchb_weighted_singles = .false. tMultiReplicas = .false. t_adjoint_replicas = .false. ! GAS options tGAS = .false. #ifdef PROG_NUMRUNS_ ! by default, excitation generation already creates matrix elements tGenMatHEl = .true. tInfSumTCCalc = .false. tInfSumTCPrint = .false. tInfSumTCRead = .false. PotentialStrength = 1.0_dp TranscorrCutoff = 0 TranscorrIntCutoff = 0 TranscorrGaussCutoff = 1.d0 TContact = .false. TUnitary = .false. Tperiodicinmom = .false. t12FoldSym = .false. t_initialized_roi = .false. inum_runs = 1 #ifdef CMPLX_ lenof_sign = 2 #else lenof_sign = 1 #endif #endif !Feb08 defaults: IF (Feb08) THEN !...add defaults... end if ! Coulomb damping function currently removed. FCOULDAMPBETA = -1.0_dp COULDAMPORB = 0 k_offset = 0.0_dp end subroutine SetSysDefaults SUBROUTINE SysReadInput(file_reader, incompletely_parsed_tokens) use input_parser_mod, only: FileReader_t, TokenIterator_t USE SymData, only: tAbelianFastExcitGen USE SymData, only: tStoreStateList use OneEInts, only: tOneElecDiag class(FileReader_t), intent(inout) :: file_reader type(TokenIterator_t), intent(inout) :: incompletely_parsed_tokens LOGICAL eof CHARACTER(LEN=100) w INTEGER I, Odd_EvenHPHF, Odd_EvenMI integer :: ras_size_1, ras_size_2, ras_size_3, ras_min_1, ras_max_3, itmp type(TokenIterator_t) :: tokens character(*), parameter :: t_r = 'SysReadInput' character(*), parameter :: this_routine = 'SysReadInput' integer :: temp_n_orbs, buf(1000), n_orb ! The system block is specified with at least one keyword on the same ! line, giving the system type being used. w = to_upper(incompletely_parsed_tokens%next()) select case (w) case ("DFREAD") ! Instead, specify DensityFitted within the system block. TREADINT = .true. TDFREAD = .true. w = to_upper(incompletely_parsed_tokens%next()) select case (w) case ("ORDER") THFORDER = .true. case ("NOORDER") THFNOORDER = .true. end select case ("BINREAD") ! Instead, specify Binary within the system block. TREADINT = .true. TBIN = .true. w = to_upper(incompletely_parsed_tokens%next()) select case (w) case ("ORDER") THFORDER = .true. case ("NOORDER") THFNOORDER = .true. end select case ("READ", "GENERIC") TREADINT = .true. w = to_upper(incompletely_parsed_tokens%next(if_exhausted='')) select case (w) case ("ORDER") THFORDER = .true. case ("NOORDER") THFNOORDER = .true. end select case ("HUBBARD") THUB = .true. TPBC = .true. if (incompletely_parsed_tokens%remaining_items() > 0) then ! this indicates the new hubbard implementation ! for consistency turn off the old hubbard indication ! and for now this is only done for the real-space hubbard ! not for the k-space implementation todo ! and do i need to turn of tpbc also? try ! use the already provided setup routine and just modify the ! necessary stuff, like excitation generators! t_new_hubbard = .true. w = to_lower(incompletely_parsed_tokens%next()) select case (w) case ('real-space', 'real') treal = .true. t_new_real_space_hubbard = .true. t_lattice_model = .true. ! if no further input is given a provided fcidump is ! assumed! but this still needs to be implemented ! this fcidump gives the lattice structure! if (incompletely_parsed_tokens%remaining_items() > 0) then lattice_type = to_lower(incompletely_parsed_tokens%next()) else lattice_type = 'read' end if case ('momentum-space', 'k-space', 'momentum') ! reuse most of the old initialisation for the k-space ! hubbard. one has to be really careful to initialize all ! the correct stuff especially for the matrix element ! calculation with the HPHF option turned on! t_k_space_hubbard = .true. t_lattice_model = .true. tKPntSym = .true. case default print *, w call Stop_All(t_r, "not recognised keyword!") end select end if case ('TJ', 'TJ-MODEL') t_tJ_model = .true. t_lattice_model = .true. ! misuse the hubbard initialisation tHub = .true. tpbc = .true. treal = .true. case ('HEISENBERG') ! should i misuse the already provided setup for the hubbard ! model again? .. maybe.. ! maybe i should use a general flag like t_lattice_model ! especially for the matrix element evaluation and stuff t_heisenberg_model = .true. t_lattice_model = .true. thub = .true. tpbc = .true. treal = .true. case ("RIINTEGRALS") tRIIntegrals = .true. tReadInt = .true. case ("UEG") TUEG = .true. tOneElecDiag = .true. !One electron integrals diagonal case ("VASP") tVASP = .true. case ("CPMD") TCPMD = .true. w = to_upper(incompletely_parsed_tokens%next()) select case (w) case ("ORDER") THFORDER = .true. end select case ("BOX") tOneElecDiag = .true. !One electron integrals diagonal case default call stop_all(this_routine, "System type "//trim(w)//" not valid") end select ! Now parse the rest of the system block. system: do while (file_reader%nextline(tokens, skip_empty=.true.)) w = to_upper(tokens%next()) select case (w) ! Options for molecular (READ) systems: control how the integral file ! is read in. case ("BINARY") tBin = .true. case ("DENSITYFITTED") tDFRead = .true. ! General options. case ("RIINTEGRALS") tRIIntegrals = .true. case ("ELECTRONS", "NEL") NEL = to_int(tokens%next()) case ("SPIN-RESTRICT") if (tokens%remaining_items() > 0) then user_input_m_s = to_int(tokens%next()) end if TSPN = .true. ! ==================== GUGA Implementation ==================== ! activate total spin preserving graphical unitary group approach and ! default total spin operator eigenvalue to 0, or else give as integer ! CONVENTION: give S in units of h/2, so S directly relates to the ! number of unpaired electrons case ("GUGA") if (tokens%remaining_items() > 0) then STOT = to_int(tokens%next()) else STOT = 0 end if tGUGA = .true. if (t_new_hubbard) then t_guga_noreorder = .true. end if tGUGACore = .true. ! also set LMS value to the inputted STOT to misuse the reference ! determinant creation for a fixed LMS also for the GUGA approach... LMS = STOT ! ============================================================= case ("TEST-DOUBLE") t_test_double = .true. test_i = to_int(tokens%next()) test_j = to_int(tokens%next()) test_k = to_int(tokens%next()) test_l = to_int(tokens%next()) case ("TEST-SINGLE") t_test_single = .true. test_i = to_int(tokens%next()) test_j = to_int(tokens%next()) case ("GUGA-NOREORDER") ! do not reorder the orbitals in the hubbard + guga implementation t_guga_noreorder = .true. case ("SYMIGNOREENERGIES") tSymIgnoreEnergies = .true. case ("NOSYMMETRY") lNoSymmetry = .true. IF (tHub) THEN CALL Stop_All("SysReadInput", "Cannot turn off symmetry with the hubbard model.") end if case ("FREEFORMAT") ! Relax the formatting requirements for reading in FCIDUMP files. ! ! For historical reasons, QChem uses a very fixed format for ! outputting FCIDUMP files. As a result the columns of orbital ! indices will merge whenever there are more than 99 spatial ! orbitals. To correctly read these files a FIXED format is ! required for reading. Obviously, this is non-ideal when reading ! FCIDUMP formats from elsewhere. ! ! The QChem behaviour used to be default, but this has been ! deprecated. To obtain the fixed behaviour use ! "FREEFORMAT OFF" or "FREEFORMAT FALSE" tReadFreeFormat = .true. if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("OFF", "FALSE") tReadFreeFormat = .false. case default end select end if case ("MIXED-HUBBARD") t_mixed_hubbard = .true. tNoBrillouin = .true. tBrillouinsDefault = .false. pParallel = 0.0_dp case ("OLLE-HUBBARD") t_olle_hubbard = .true. tNoBrillouin = .true. tBrillouinsDefault = .false. pParallel = 0.0_dp case ("SYM") TPARITY = .true. do I = 1, 4 IPARITY(I) = to_int(tokens%next()) end do ! the last number is the symmetry specification - and is placed in position 5 IPARITY(5) = IPARITY(4) IPARITY(4) = 0 tSymSet = .true. case ("USEBRILLOUINTHEOREM") TUSEBRILLOUIN = .TRUE. tNoBrillouin = .false. tBrillouinsDefault = .false. case ("NOBRILLOUINTHEOREM") tNoBrillouin = .true. tBrillouinsDefault = .false. case ("UHF") ! indicates UHF type FCIDUMP ! Note, this keyword is required if we are doing an open shell calculation ! but do not want to include singles in the energy calculations ! (e.g. by nobrillouintheorem) tUHF = .true. case ("RS") FUEGRS = to_realdp(tokens%next()) case ("EXCHANGE-CUTOFF") iPeriodicDampingType = 2 if (tokens%remaining_items() > 0) then fRc = to_realdp(tokens%next()) end if case ("EXCHANGE-ATTENUATE") iPeriodicDampingType = 1 if (tokens%remaining_items() > 0) then fRc = to_realdp(tokens%next()) end if case ("EXCHANGE") w = to_upper(tokens%next()) select case (w) case ("ON") TEXCH = .TRUE. case ("OFF") TEXCH = .FALSE. case default call stop_all(this_routine, "EXCHANGE "//trim(w)//" not valid") end select case ("COULOMB") call stop_all(this_routine, "Coulomb feature removed") case ("COULOMB-DAMPING") call stop_all(this_routine, "Coulomb damping feature removed") case ("ENERGY-CUTOFF") tOrbECutoff = .true. OrbECutoff = to_realdp(tokens%next()) case ("G-CUTOFF") tgCutoff = .true. gCutoff = to_realdp(tokens%next()) case ("FREEZE-CUTOFF") tUEGFreeze = .true. FreezeCutoff = to_realdp(tokens%next()) case ("MADELUNG") tMadelung = .true. Madelung = to_realdp(tokens%next()) case ("UEG2") tUEG2 = .true. case ("STORE-AS-EXCITATIONS") tStoreAsExcitations = .true. case ("MP2-UEG-RESTRICT") tMP2UEGRestrict = .true. kiRestrict(1) = to_int(tokens%next()) kiRestrict(2) = to_int(tokens%next()) kiRestrict(3) = to_int(tokens%next()) kiMsRestrict = to_int(tokens%next()) kjRestrict(1) = to_int(tokens%next()) kjRestrict(2) = to_int(tokens%next()) kjRestrict(3) = to_int(tokens%next()) kjMsRestrict = to_int(tokens%next()) ! Options for model systems (electrons in a box/Hubbard). case ("CELL") NMAXX = to_int(tokens%next()) NMAXY = to_int(tokens%next()) NMAXZ = to_int(tokens%next()) case ('SPIN-TRANSCORR') ! make a spin-dependent transcorrelation factor t_spin_dependent_transcorr = .true. if (tokens%remaining_items() > 0) then trans_corr_param = to_realdp(tokens%next()) else trans_corr_param = 0.1_dp end if t_non_hermitian_2_body = .true. case ("NONHERMITIAN") ! just use a non-hermitian Hamiltonian, no additional tweaks ! note transcorrelation has only nonhermitian 2-body integrals ! so if you are doing a TCMF calculation, do ! `nonhermitian 2-body` tNoBrillouin = .true. tBrillouinsDefault = .false. if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("1-BODY") write(stdout, '(A)') 'Treating 1-body integrals as non-Hermitian.' t_non_hermitian_1_body = .true. case ("2-BODY") write(stdout, '(A)') 'Treating 2-body integrals as non-Hermitian.' t_non_hermitian_2_body = .true. case default write(stdout, '(A)') 'Treating all integrals as non-Hermitian.' ! by default, do both t_non_hermitian_1_body = .true. t_non_hermitian_2_body = .true. end select end if case("ADJOINT-CALCULATION") ! calculate the adjoint of H instead of H t_calc_adjoint = .true. write(stdout, *) "Calculating H^\dagger instead of H." if (.not. t_non_hermitian_1_body .and. .not. t_non_hermitian_2_body) then write(stdout, *) "WARNING: Adjoint matrix calculation for a & &Hermitian matrix. Assuming this is intended." end if case ('MOLECULAR-TRANSCORR') tNoBrillouin = .true. tBrillouinsDefault = .false. t_non_hermitian_2_body = .true. ! optionally supply the three-body integrals of the TC Hamiltonian t_3_body_excits = .true. if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("3-BODY") t_mol_3_body = .true. max_ex_level = 3 ! this uses a uniform excitation generator, switch off matrix ! element computation for HPHF tGenMatHEl = .false. case ("UEG") t_mol_3_body = .true. tGenMatHEl = .false. t12FoldSym = .true. tNoSinglesPossible = .true. case default t_mol_3_body = .true. tGenMatHEl = .false. end select end if if (t_mol_3_body) max_ex_level = 3 case ('UEG-TRANSCORR') t_ueg_transcorr = .true. t_non_hermitian_2_body = .true. do while (tokens%remaining_items() > 0) w = to_upper(tokens%next()) select case (w) case ("3-BODY") tTrcorrExgen = .false. tTrCorrRandExgen = .true. t_ueg_3_body = .true. tGenMatHEl = .false. max_ex_level = 3 tRPA_tc = .false. case ("TRCORR-EXCITGEN") tTrcorrExgen = .true. tTrCorrRandExgen = .false. case ("RAND-EXCITGEN") tTrCorrRandExgen = .true. end select end do case ('UEG-DUMP') t_ueg_dump = .true. case ('EXCLUDE-3-BODY-EX') ! Do not generate 3-body excitations, even in the molecular-transcorr mode t_exclude_3_body_excits = .true. case ('EXCLUDE-3-BODY-PARALLEL', 'EXCLUDE-3-BODY-PAR') ! exclude fully spin-parallel 3 body excitation t_exclude_pure_parallel = .true. case ('TRANSCORRELATED', 'TRANSCORR', 'TRANS-CORR') ! activate the transcorrelated Hamiltonian idea from hongjun for ! the real-space hubbard model t_trans_corr = .true. t_non_hermitian_2_body = .true. if (tokens%remaining_items() > 0) then trans_corr_param = to_realdp(tokens%next()) else ! defaul value 1 for now, since i have no clue how this behaves trans_corr_param = 1.0_dp end if case ("TRANSCORR-NEW") t_trans_corr = .true. t_trans_corr_new = .true. t_non_hermitian_2_body = .true. if (tokens%remaining_items() > 0) then trans_corr_param = to_realdp(tokens%next()) else ! defaul value 1 for now, since i have no clue how this behaves trans_corr_param = 1.0_dp end if case ('2-BODY-TRANSCORR', '2-BODY-TRANS-CORR', '2-BODY-TRANSCORRELATED', 'TRANSCORR-2BODY') ! for the tJ model there are 2 choices of the transcorrelation ! indicate that here! t_trans_corr_2body = .true. t_non_hermitian_2_body = .true. if (tokens%remaining_items() > 0) then trans_corr_param_2body = to_realdp(tokens%next()) else trans_corr_param_2body = 0.25_dp end if ! if it is the k-space hubbard also activate 3-body excitations here if (t_k_space_hubbard) then t_3_body_excits = .true. max_ex_level = 3 end if case ('NEIGHBOR-TRANSCORR', 'TRANSCORR-NEIGHBOR', 'N-TRANSCORR') t_trans_corr_2body = .true. t_non_hermitian_2_body = .true. if (tokens%remaining_items() > 0) then trans_corr_param_2body = to_realdp(tokens%next()) else trans_corr_param_2body = 0.25_dp end if ! if it is the k-space hubbard also activate 3-body excitations here if (t_k_space_hubbard) then t_3_body_excits = .true. max_ex_level = 3 end if case ("TRANSCORR-HOP", "HOP-TRANSCORR") t_trans_corr_hop = .true. t_non_hermitian_2_body = .true. if (tokens%remaining_items() > 0) then trans_corr_param = to_realdp(tokens%next()) else trans_corr_param = 0.5_dp end if case ("HOLE-FOCUS") t_hole_focus_excits = .true. if (tokens%remaining_items() > 0) then pholefocus = to_realdp(tokens%next()) else pholefocus = 0.5_dp end if case ("PRECOND-HUB") t_precond_hub = .true. case ("NO_REF_SHIFT") t_no_ref_shift = .true. ! Options for the type of the reciprocal lattice (eg sc, fcc, bcc) case ("REAL_LATTICE_TYPE") real_lattice_type = to_lower(tokens%next()) ! Options for the dimension (1, 2, or 3) case ("DIMENSION") dimen = to_int(tokens%next()) ! Options for transcorrelated method (only: UEG 2D 3D, Homogeneous 1D 3D ! gas with contact interaction) case ("TRANSCORRCUTOFF") if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("GAUSS") if (dimen /= 1) stop 'Gauss cutoff is developed only for 1D!' TranscorrGaussCutoff = to_realdp(tokens%next()) t_trcorr_gausscutoff = .true. case ("STEP") t_trcorr_gausscutoff = .false. TranscorrCutoff = to_int(tokens%next()) if (tContact .and. dimen == 3) then tInfSumTCCalc = .true. TranscorrIntCutoff = to_int(tokens%next()) end if end select end if !Options for turn off the RPA term(only: transcorrelated homogeneous 1D !gas with contact interaction), tRPA_tc is set to true for two particles, ! but it is turned off, if 3-body interactions are used case ("NORPATC") tRPA_tc = .false. case ("PERIODICINMOMSPACE") Tperiodicinmom = .true. ! Contact interaction for homogenous one dimensional Fermi gas is applied case ("CONTACTINTERACTION") tContact = .true. PotentialStrength = to_realdp(tokens%next()) if (dimen /= 1) & stop 'Contact interaction only for 1D!' ! Contact interaction forthe three dimensional Fermi gas in the unitary ! interaction case ("CONTACTINTERACTIONUNITARY") tContact = .true. tUnitary = .true. if (dimen /= 3) stop 'Unitary regime only for 3D!' ! Option for infinite summation at transcorrelated method for homogeneous ! 3D gas with contact interaction case ("TRANSCORRINFSUM") if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("CALC") tInfSumTCCalc = .true. case ("CALCANDPRINT") tInfSumTCCalc = .true. tInfSumTCPrint = .true. case ("READ") tInfSumTCCalc = .false. tInfSumTCRead = .true. end select end if ! This means that no a is generated when b would be made and rejected ! O(N^2) loop makes this a poor choice for larger systems. case ("UEGNOFAIL") tNoFailAb = .true. ! These are the new lattice excitation generators that conserve momentum ! during excitation generation for efficiency case ("LATTICE-EXCITGEN") tLatticeGens = .true. ! use the simplified random excitation generator for k-space hubbard that ! does not use a cumulative list (it is much more efficient) case ("UNIFORM-EXCITGEN") t_uniform_excits = .true. case ("MIXED-EXCITGEN") ! use a mix of weighted and uniform excitation generators ! for the transcorrelated k-space hubbard. ! triples are weighted and doubles will be done uniformly t_mixed_excits = .true. case ("MESH") NMSH = to_int(tokens%next()) case ("BOXSIZE") BOX = to_realdp(tokens%next()) if (tokens%remaining_items() > 0) then BOA = to_realdp(tokens%next()) COA = to_realdp(tokens%next()) else BOA = 1.0_dp COA = 1.0_dp end if case ("U") UHUB = to_realdp(tokens%next()) case ("B") BHUB = to_realdp(tokens%next()) case ("J") ! specify the tJ exchange here, the default is 1.0 ! this could also be used for the heisenberg model.. exchange_j = to_realdp(tokens%next()) case ("C") btHub = to_realdp(tokens%next()) tmodHub = .true. case ("NEXT-NEAREST-HOPPING") nn_bhub = to_realdp(tokens%next()) case ("TWISTED-BC") t_twisted_bc = .true. twisted_bc(1) = to_realdp(tokens%next()) if (tokens%remaining_items() > 0) then twisted_bc(2) = to_realdp(tokens%next()) end if case ("ANTI-PERIODIC-BC") if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("X") t_anti_periodic(1) = .true. case ("Y") t_anti_periodic(2) = .true. case ("XY", "YX") t_anti_periodic(1:2) = .true. end select else t_anti_periodic(1:2) = .true. end if case ("REAL") TREAL = .true. ! i think for the real-space i have to specifically turn off ! the symmetry and make other changes in the code to never try ! to set up the symmetry ! in case of the real-space lattice also turn off symmetries lNoSymmetry = .true. case ("APERIODIC") TPBC = .false. case ("OPEN-BC") ! open boundary implementation for the real-space hubbard ! model ! only applicable for 1 and 2 dimensional lattices! ! for the square lattice on can specify mixed boundary conditions ! (periodic + open) but not in the tilted where always full ! open boundary conditions are used if this keyword is present! if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("X") ! onl open in x-direction t_open_bc_x = .true. case ("Y") ! only open in y-direction t_open_bc_y = .true. case ("XY") ! open in both directions t_open_bc_x = .true. t_open_bc_y = .true. end select else t_open_bc_x = .true. t_open_bc_y = .true. end if case ("BIPARTITE", "BIPARTITE-ORDER") t_bipartite_order = .true. if (tokens%remaining_items() > 0) then t_input_order = .true. temp_n_orbs = to_int(tokens%next()) allocate(orbital_order(temp_n_orbs), source = 0) if (file_reader%nextline(tokens, skip_empty=.false.)) then do i = 1, temp_n_orbs orbital_order(i) = to_int(tokens%next()) end do else call stop_all(this_routine, 'Unexpected EOF reached.') end if end if case ("LATTICE") ! new hubbard implementation ! but maybe think of a better way to init that.. ! the input has to be like: ! lattice [type] [len_1] [*len_2] ! where length to is optional if it is necessary to input it. ! support t_lattice_model = .true. ! set some defaults: lattice_type = "read" length_x = -1 length_y = -1 if (tokens%remaining_items() > 0) then ! use only new hubbard flags in this case lattice_type = to_lower(tokens%next()) end if if (tokens%remaining_items() > 0) then length_x = to_int(tokens%next()) end if if (tokens%remaining_items() > 0) then length_y = to_int(tokens%next()) end if if (tokens%remaining_items() > 0) then length_z = to_int(tokens%next()) end if ! maybe i have to reuse the cell input functionality or set it ! here also, so that the setup is not messed up ! i should set those quantities here again.. nmaxx = length_x nmaxy = length_y nmaxz = 1 case ("UEG-OFFSET") tUEGOffset = .true. k_offset(1) = to_realdp(tokens%next()) k_offset(2) = to_realdp(tokens%next()) k_offset(3) = to_realdp(tokens%next()) case ("UEG-SCALED-ENERGIES") tUEGTrueEnergies = .true. case ("UEG-MOMENTUM") tUEGSpecifyMomentum = .true. k_momentum(1) = to_int(tokens%next()) k_momentum(2) = to_int(tokens%next()) k_momentum(3) = to_int(tokens%next()) case ("TILT") TTILT = .true. ITILTX = to_int(tokens%next()) ITILTY = to_int(tokens%next()) case ("ALPHA") TALPHA = .true. ALPHA = to_realdp(tokens%next()) case ("STATE") ISTATE = to_int(tokens%next()) if (ISTATE /= 1) then call stop_all(this_routine, "Require ISTATE to be left set as 1") end if case ("STOQUASTIZE") tStoquastize = .true. case ("FAST-EXCITGEN") tAbelianFastExcitGen = .true. ! tAbelianFastExcitGen is a temporary flag. ! It is used to indicate that if an Abelian symmetry group is presents ! the excitation generators should use optimized routines ! to take this into account. Not all excitation generator functions ! currently work with this. USE WITH CARE if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("OFF") tAbelianFastExcitGen = .false. end select end if case ("NORENORMRANDEXCITS") ! Since we have already calculated the number of excitations possible for each symmetry type, there ! no need to renormalise all excitations with weight 1. As long as pairs of allowed occupied and ! virtual orbitals can be chosen without any bias, then we can generate random excitations in O[1] time. ! This is default off since it will change previous results, however it is strongly recommended to be ! on for virtually all unweighted MC calculations, since it should speed up generation, especially in ! low symmetry and/or large systems. tNoRenormRandExcits = .true. case ("STORESTATELIST") !This flag indicates that we want to store the full list of symmetry state pairs. !This is done by default in non-abelian symmetries, but in abelian symmetries, will !only be done if specified. This may be wanted, since it means that random excitations !will be created quicker currently, since NoRenormRandExcits can only work with it on in !abelian symmetry. tStoreStateList = .true. case ("ASSUMESIZEEXCITGEN") !This flag indicates that we want to setup the excitations faster, by returning a maximum size of the !generators instantly in the first setup, and then we do not need to run through the excitations twice. !However, certain bits of the generator are not stored, namely nAllowPPS and SymProds, as well as the !iterator information. This means that we can only randomly generate excitations, and any attempt to !use these assumed size excitation generators to generate the whole list of excitations, will result !in bad, bad times. tAssumeSizeExcitgen = .true. case ("MOMINVSYM") call stop_all(t_r, 'Deprecated function. Look in defunct_code folder if you want to see it') case ("HPHF") tHPHF = .true. if (tokens%remaining_items() > 0) then Odd_EvenHPHF = to_int(tokens%next()) if (Odd_EvenHPHF == 1) then !Want to converge onto an Odd S State tOddS_HPHF = .true. else if (Odd_EvenHPHF == 0) then !Want to converge onto an Even S State !tOddS_HPHF should be false by default. else call stop_all("SysReadInput", "Invalid variable given to HPHF option: 0 = Even S; 1 = Odd S") end if end if case ("ROTATEORBS") ! The ROTATEORBS calculation initiates a routine which takes the HF orbitals ! and finds the optimal set of transformation coefficients to fit a particular criteria specified below. ! This new set of orbitals can then used to produce a ROFCIDUMP file and perform the FCIMC calculation. tRotateOrbs = .true. if (tokens%remaining_items() > 0) then TimeStep = to_realdp(tokens%next()) ConvergedForce = to_realdp(tokens%next()) end if ! The SHAKE orthonormalisation algorithm is automatically turned on with a default of 5 iterations. tShake = .true. tShakeIter = .true. case ("DIAGONALIZEHIJ") ! Instead of doing an orbital rotation according to the P.E's below, this keyword sets the rotate orbs routine to ! take the input <i|h|j>, and diagonalise it to give the rotation coefficients. tDiagonalizehij = .true. case ("HFSINGDOUBEXCMAX") ! This maximises the square of the four index integrals corresponding to single and double excitations from ! the HF. tHFSingDoubExcMax = .true. MaxMinFac = -1 case ("OFFDIAGSQRDMIN") tOffDiagSqrdMin = .true. MaxMinFac = 1 if (tokens%remaining_items() > 0) then OffDiagWeight = to_realdp(tokens%next()) ELSE OffDiagWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which minimise the sum of the |<ij|kl>|^2 elements. ! The following real value sets the importance of the minimisation relative to the ER maximisation, providing ! the ERLOCALIZATION keyword is also present. This is the same for all OFFDIAG keywords below. case ("OFFDIAGSQRDMAX") tOffDiagSqrdMax = .true. MaxMinFac = -1 if (tokens%remaining_items() > 0) then OffDiagWeight = to_realdp(tokens%next()) ELSE OffDiagWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which maximise the sum of the |<ij|kl>|^2 elements. case ("OFFDIAGMIN") tOffDiagMin = .true. MaxMinFac = 1 if (tokens%remaining_items() > 0) then OffDiagWeight = to_realdp(tokens%next()) ELSE OffDiagWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which minimise the sum of the <ij|kl> elements. case ("OFFDIAGMAX") tOffDiagMax = .true. MaxMinFac = -1 if (tokens%remaining_items() > 0) then OffDiagWeight = to_realdp(tokens%next()) ELSE OffDiagWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which maximise the sum of the <ij|kl> elements. case ("DOUBEXCITEMIN") tDoubExcMin = .true. MaxMinFac = 1 if (tokens%remaining_items() > 0) then OffDiagWeight = to_realdp(tokens%next()) ELSE OffDiagWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which minimise the double excitation hamiltonian elements. case ("HIJSQRDMIN") tHijSqrdMin = .true. OneElMaxMinFac = 1 tRotateVirtOnly = .true. if (tokens%remaining_items() > 0) then OneElWeight = to_realdp(tokens%next()) ELSE OneElWeight = 1.0 end if ! This sets the orbital rotation to find the coefficients which minimis the square of the one electron integrals, <i|h|j>. ! i can be occupied or virtual, j only virtual, i=<j. case ("ONEELINTMAX") tOneElIntMax = .true. MaxMinFac = -1 tRotateVirtOnly = .true. ! This sets the orbital rotation to find the coefficients which maximise the one electron integrals, <i|h|i>. case ("ONEPARTORBENMAX") tOnePartOrbEnMax = .true. MaxMinFac = -1 tRotateVirtOnly = .true. if (tokens%remaining_items() > 0) then OrbEnMaxAlpha = to_realdp(tokens%next()) ELSE OrbEnMaxAlpha = 1.0_dp end if ! This sets the orbital rotation to find the coefficients which maximise the one particle orbital energies. ! i.e maximise the fock energies, epsilon_i = <i|h|i> + sum_j [<ij||ij>]. case ("ERLOCALIZATION") tERLocalization = .true. DiagMaxMinFac = -1 if (tokens%remaining_items() > 0) then DiagWeight = to_realdp(tokens%next()) ELSE DiagWeight = 1.0 end if ! This sets the orbital rotation to an Edmiston-Reudenberg localisation. This maximises the self repulsion energy, i.e ! maximises the sum of the <ii|ii> terms. case ("VIRTCOULOMBMAX") tVirtCoulombMax = .true. MaxMinFac = -1 tRotateVirtOnly = .true. ! This sets the orbital rotation to maximise the sum_i_j <ij|ij> terms where both i and j are virtual spatial orbitals. case ("VIRTEXCHANGEMIN") tVirtExchangeMin = .true. MaxMinFac = 1 tRotateVirtOnly = .true. ! This sets the orbital rotation to minimise the sum_i_j <ij|ji> terms where both i and j are virtual spatial orbitals. case ("SHAKE") ! This will use the shake algorithm to iteratively enforce orthonormalisation ! on the rotation coefficients calculated in the ROTATEORBS ! routine. It finds a force matrix which moves the coefficients at a tangent ! to the constraint surface, from here, a normalisation ! will require just a small adjustment to ensure complete orthonormalisation, ! but not majorly affecting the new coefficients. tShake = .true. if (tokens%remaining_items() > 0) then ShakeConverged = to_realdp(tokens%next()) end if case ("SHAKEAPPROX") ! This turns on the shake approximation algorithm. ! To be used if the matrix inversion required in the full shake algorithm cannot ! be performed. ! The approximation applies the iterative scheme to find lambda, ! to each constraint in succession, rather than simultaneously. tShake = .false. tShakeApprox = .true. case ("SHAKEITER") ! Much like 'ROIteration', this overwrites the convergence criteria for the iterations of the shake constraints ! and instead performs only the number of iterations specified on this line. tShakeIter = .true. ShakeIterMax = to_int(tokens%next()) case ("SHAKEDELAY") ! This option sets the shake orthonomalisation algorithm to only kick in after a certain number of rotatation iterations. This ! potentially allows a large shift in the coefficients away from their starting point, followed by orthonormalisation to a ! significantly different position. tShakeDelay = .true. ShakeStart = to_int(tokens%next()) case ("LAGRANGE") ! This will use a non-iterative lagrange multiplier for each component of each rotated vector ! in the rotateorbs routines in order to ! attempt to maintain orthogonality. This currently does not seem to work too well! tShake = .false. tLagrange = .true. case ("SEPARATEOCCVIRT") ! This option applies to the orbital rotation. ! If present, the virtual and occuppied orbitals are localised separately. ! This has the advantage of keeping the HF reference determinant the same. tSeparateOccVirt = .true. case ("ROTATEOCCONLY") ! This option applies to orbital rotation. It separates the orbitals into occupied and virtual, ! and rotates the occupied while keeping the virtual as the HF. tSeparateOccVirt = .true. tRotateOccOnly = .true. case ("ROTATEVIRTONLY") ! This option rotates the virtual orbitals while keeping the occupied as the HF. tSeparateOccVirt = .true. tRotateVirtOnly = .true. case ("ROITERATION") ! Specifying this keyword overwrites the convergence limit from the ROTATEORBS line, ! and instead runs the orbital rotation for as many ! iterations as chosen on this line. tROIteration = .true. ROIterMax = to_int(tokens%next()) case ("MAXHLGAP") tMaxHLGap = .true. case ("ROTATEDORBS") ! This is a flag which is to be included if the FCIDUMP being read in is one containing rotated orbitals. ! The only difference is that the ! H elements for single excitations are no longer 0 (as for HF), ! and walkers on singly excited determinants must be included in the energy ! calculations. tRotatedOrbs = .true. tNoBrillouin = .true. tBrillouinsDefault = .false. case ("SPINORBS") ! This flag simply uses spin orbitals to perform the rotation rather than spatial orbitals. ! IF UHF=.true. is present in the FCIDUMP file, this will happen automatically, otherwise this keyword is required. tSpinOrbs = .true. case ("READINTRANSMAT") ! This sets the rotation routine to read in a transformation matrix (coefft1) given by the file TRANSFORMMAT, ! and use the transformation ! routine and print rofcidump routines to transform the orbitals and print out a new dump file. tReadInCoeff = .true. case ("USEMP2VDM") ! Once in the rotation routine, the MP2 variational density matrix is calculated, ! and this is used to transform the orbitals and print and ! new FCIDUMP file. tUseMP2VarDenMat = .true. tSeparateOccVirt = .true. tRotateVirtOnly = .true. tShake = .false. case ("USECINATORBS") ! This rotation option is slightly different, it first requires a spawning calculation ! from which the amplitudes of the wavefunction are ! obtained. From here, the one electron reduced density matrix is calculated, and the eigenvectors ! of this are used to rotate the HF orbitals. ! A new ROFCIDUMP file is then produced in the new natural orbitals. tFindCINatOrbs = .true. tShake = .false. case ("USEHFORBS") tUseHFOrbs = .true. tShake = .false. tSeparateOccVirt = .true. case ("RANLUXLEV") !This is the level of quality for the random number generator. Values go from 1 -> 4. 3 is default. iRanLuxLev = to_int(tokens%next()) case ("CALCEXACTSIZESPACE") !This option will calculate the exact size of the symmetry allowed space of determinants. Will scale badly. tExactSizeSpace = .true. case ("CALCMCSIZETRUNCSPACE") !This option will approximate the exact size of the symmetry allowed truncated space of determinants by MC. !The variance on the value will decrease as 1/N_steps tMCSizeTruncSpace = .true. iMCCalcTruncLev = to_int(tokens%next()) CalcDetCycles = to_int64(tokens%next()) CalcDetPrint = to_int64(tokens%next()) case ("CALCMCSIZESPACE") !This option will approximate the exact size of the symmetry allowed space of determinants by MC. ! The variance on the value will decrease as 1/N_steps tMCSizeSpace = .true. CalcDetCycles = to_int64(tokens%next()) CalcDetPrint = to_int64(tokens%next()) case ("NONUNIFORMRANDEXCITS") !This indicates that the new, non-uniform O[N] random excitation generators are to be used. !CYCLETHRUORBS can be useful if we have small basis sets or v high restrictive symmetry and will eliminate !large numbers of unsuccessful random draws of orbitals by calculating how many allowed orbitals there are !and cycling through them until the allowed one is drawn, rather than randomly drawing and redrawing until !an allowed orbital is found. For large basis sets, the chance of drawing a forbidden orbital is small !enough that this should be an unneccesary expense. tNonUniRandExcits = .true. do while (tokens%remaining_items() > 0) w = to_upper(tokens%next()) non_uniform_rand_excits: select case (w) case ("NOSYM_GUGA") call Stop_All(this_routine, "'nosym-guga' option deprecated!") case ("NOSYM_GUGA_DIFF") call Stop_All(this_routine, "'nosym-guga-diff' option deprecated!") case ("UEG_GUGA", "UEG-GUGA") tGen_sym_guga_ueg = .true. if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("MIXED") tgen_guga_mixed = .true. if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("SEMI") t_guga_mixed_init = .false. t_guga_mixed_semi = .true. end select end if case ("CRUDE") tgen_guga_crude = .true. end select end if case ("MOL_GUGA") tGen_sym_guga_mol = .true. case ("MOL_GUGA_WEIGHTED", "MOL-GUGA-WEIGHTED") tGen_sym_guga_mol = .true. tgen_guga_weighted = .true. case ("GUGA-CRUDE") ! try a crude excitation approximation, where no ! spin-flips in the excitation range are allowed tgen_guga_crude = .true. if (t_k_space_hubbard) then tGen_sym_guga_ueg = .true. else tgen_guga_weighted = .true. tGen_sym_guga_mol = .true. end if case ('GUGA-MIXED') ! try a mix of the crude and full implementation: ! initiators spawn fully, while non-initiators ! only spawn in the crude approximation tgen_guga_mixed = .true. tgen_guga_crude = .true. if (t_k_space_hubbard) then tGen_sym_guga_ueg = .true. else tgen_guga_weighted = .true. tGen_sym_guga_mol = .true. end if if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ('SEMI') t_guga_mixed_init = .false. t_guga_mixed_semi = .true. end select end if case ('GUGA-APPROX-EXCHANGE') ! in this case I force an exchange at the ! chosen indices of exchange type and 3-ind ! excitations to reduce the complexity of the ! algorithm t_approx_exchange = .true. if (t_k_space_hubbard) then tGen_sym_guga_ueg = .true. else tgen_guga_weighted = .true. tGen_sym_guga_mol = .true. end if if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ('NON-INITS') ! only do the approx. for noninits t_approx_exchange = .false. t_approx_exchange_noninits = .true. end select end if case ('GUGA-CRUDE-EXCHANGE') ! calculate exchange type excitation like ! determinants.. for all states, initiators and ! non-inits (also for 3-index excitation of ! mixed type) t_crude_exchange = .true. if (t_k_space_hubbard) then tGen_sym_guga_ueg = .true. else tgen_guga_weighted = .true. tGen_sym_guga_mol = .true. end if if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ('NON-INITS') ! only to the approx. for non-inits t_crude_exchange = .false. t_crude_exchange_noninits = .true. end select end if case ("CYCLETHRUORBS") tCycleOrbs = .true. case ("NOSYMGEN") !Do not generate symmetry-allowed excitations only, but all excitations. Spin-symmetry is still taken into account. tNoSymGenRandExcits = .true. case ("IMPORTANCESAMPLE") !Importance sample the excitations for FCIMCPar CALL Stop_All("ReadSysInp", "IMPORTANCESAMPLE option deprecated") ! tImportanceSample=.true. case ("PICK-VIRT-UNIFORM") ! Pick virtual orbitals randomly and uniformly in the ! 3rd generation of random excitation generators ! (symrandexcit3.F90) tPickVirtUniform = .true. tBuildOccVirtList = .true. case ("PICK-VIRT-UNIFORM-MAG") ! Pick virtual orbitals randomly and uniformly in the ! 3rd generation of random excitation generators ! (symrandexcit3.F90) tPickVirtUniform = .true. tBuildOccVirtList = .true. tBuildSpinSepLists = .true. case ("HEL-WEIGHTED-SLOW") ! Pick excitations from any site with a generation ! probability proportional to the connectiong HElement ! --> Lots of enumeration ! --> Very slow ! --> Maximal possible values of tau. tGenHelWeighted = .true. case ("4IND-WEIGHTED") ! Weight excitations based on the magnitude of the ! 4-index integrals (ij|ij) tGen_4ind_weighted = .true. case ("4IND-REVERSE") ! Weight excitations based on the magnitude of the ! actual hamiltonian matrix elements (at least for ! doubles). This is effectively the "reverse" of ! 4IND-WEIGHTED as above. tGen_4ind_reverse = .true. case ("4IND-WEIGHTED-PART-EXACT") ! Weight excitations as in 4IND-WEIGHTED, except for ! double excitations with the same spin which are ! weighted according to: ! sqrt(((ii|aa) + (jj|aa))(<ij|ab>-<ij|ba>)) tGen_4ind_weighted = .true. tGen_4ind_part_exact = .true. case ("4IND-WEIGHTED-LIN-EXACT") ! Weight excitations as in 4IND-WEIGHTED, except for ! double excitations with the same spin which are ! weighted according to: ! (1/M)(<ij|ab> - <ij|ba>) ! (The second half of this only affecting the choice ! of electron b) tGen_4ind_weighted = .true. tGen_4ind_lin_exact = .true. case ("4IND-WEIGHTED-2") ! Second attempt at 4ind weighted generator ! ! This version disables symmetric selection of the ! orbitals. I.e. for non-parallel electron choice, the ! orbitals are chosen one spin first, then the other ! ! --> Very slightly better Hij/pgen ratio ! --> Much faster runtime tGen_4ind_2 = .true. tGen_4ind_part_exact = .true. tGen_4ind_2_symmetric = .false. case ("4IND-WEIGHTED-UNBOUND") ! [Werner Dobrautz 26.4.2017:] ! new tests for optimizations for the latest ! excitation generator implementation, without using ! the artificial lower bounds for the energy tGen_4ind_2 = .true. tGen_4ind_part_exact = .true. tGen_4ind_2_symmetric = .false. tGen_4ind_unbound = .true. ! make a few small tests for the frequency histograms if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case (w) case ("IIAA") ! weight with <ii|aa> instead of <ia|ia> t_iiaa = .true. case ("RATIO") ! weigh with the ratio <ia|ia>/<ja|ja> t_ratio = .true. case ("IIAA-RATIO", "RATIO-IIAA") t_iiaa = .true. t_ratio = .true. end select end if case ("4IND-WEIGHTED-2-SYMMETRIC") ! The other version of this generator. This permits ! selecting orbitals in both directions tGen_4ind_2 = .true. tGen_4ind_part_exact = .true. tGen_4ind_2_symmetric = .true. case ("PCPP") ! the precomputed power-pitzer excitation generator t_pcpp_excitgen = .true. case ("PCHB") t_fci_pchb_excitgen = .true. if (allocated(FCI_PCHB_user_input)) deallocate(FCI_PCHB_user_input) allocate(FCI_PCHB_user_input) if (tokens%remaining_items() == 0) then call stop_all(t_r, "No item specified. Please specify one of {LOCALISED, DELOCALISED, MANUAL}.") end if select case (to_upper(tokens%next())) case ("LOCALISED") FCI_PCHB_user_input%option_selection = FCI_PCHB_user_input_vals%option_selection%LOCALISED case ("DELOCALISED") FCI_PCHB_user_input%option_selection = FCI_PCHB_user_input_vals%option_selection%DELOCALISED case ("MANUAL") FCI_PCHB_user_input%option_selection = FCI_PCHB_user_input_vals%option_selection%MANUAL allocate(FCI_PCHB_user_input%options) FCI_PCHB_user_input%options%singles = FCI_PCHB_user_input_vals%options%singles%from_str(tokens%next()) FCI_PCHB_user_input%options%doubles = FCI_PCHB_user_input_vals%options%doubles%from_str(tokens%next()) if (tokens%remaining_items() == 1) then select case (to_upper(tokens%next())) case("SPIN-ORB-RESOLVED") FCI_PCHB_user_input%options%doubles%spin_orb_resolved = .true. case("NO-SPIN-ORB-RESOLVED") FCI_PCHB_user_input%options%doubles%spin_orb_resolved = .false. case default call stop_all(t_r, "Should not be here.") end select else if (tokens%remaining_items() > 1) then call stop_all(t_r, "Should not be here.") end if case default write(stderr, *) "Please specify one of {LOCALISED, DELOCALISED, MANUAL}." write(stderr, *) "With localised and delocalised the fastest PCHB sampler is chosen automatically." write(stderr, *) "Specify delocalised for Hartree-Fock like orbitals and" write(stderr, *) "localised for localised orbitals." write(stderr, *) "With manual one can select the PCHB sampler themselves." call stop_all(t_r, "Wrong input to PCHB.") end select case ("GAS-CI") w = to_upper(tokens%next()) select case (w) case ('ON-THE-FLY-HEAT-BATH') user_input_GAS_exc_gen = possible_GAS_exc_gen%ON_FLY_HEAT_BATH case ('DISCONNECTED') user_input_GAS_exc_gen = possible_GAS_exc_gen%DISCONNECTED case ('DISCARDING') user_input_GAS_exc_gen = possible_GAS_exc_gen%DISCARDING case ('PCHB') user_input_GAS_exc_gen = possible_GAS_exc_gen%PCHB if (allocated(GAS_PCHB_user_input)) deallocate(GAS_PCHB_user_input) allocate(GAS_PCHB_user_input) if (tokens%remaining_items() == 0) then call stop_all(t_r, "No item specified. Please specify one of {LOCALISED, DELOCALISED, MANUAL}.") end if select case (to_upper(tokens%next())) case ("LOCALISED") GAS_PCHB_user_input%option_selection = GAS_PCHB_user_input_vals%option_selection%LOCALISED case ("DELOCALISED") GAS_PCHB_user_input%option_selection = GAS_PCHB_user_input_vals%option_selection%DELOCALISED case ("MANUAL") GAS_PCHB_user_input%option_selection = GAS_PCHB_user_input_vals%option_selection%MANUAL allocate(GAS_PCHB_user_input%options) GAS_PCHB_user_input%options%use_lookup = .true. GAS_PCHB_user_input%options%singles = GAS_PCHB_user_input_vals%options%singles%from_str(tokens%next()) GAS_PCHB_user_input%options%doubles = GAS_PCHB_user_input_vals%options%doubles%from_str(tokens%next()) if (tokens%remaining_items() == 1) then select case (to_upper(tokens%next())) case("SPIN-ORB-RESOLVED") GAS_PCHB_user_input%options%doubles%spin_orb_resolved = .true. case("NO-SPIN-ORB-RESOLVED") GAS_PCHB_user_input%options%doubles%spin_orb_resolved = .false. case default call stop_all(t_r, "Should not be here.") end select else if (tokens%remaining_items() > 1) then call stop_all(t_r, "Should not be here.") end if case default write(stderr, *) "Please specify one of {LOCALISED, DELOCALISED, MANUAL}." write(stderr, *) "With localised and delocalised the fastest PCHB sampler is chosen automatically." write(stderr, *) "Specify delocalised for Hartree-Fock like orbitals and" write(stderr, *) "localised for localised orbitals." write(stderr, *) "With manual one can select the PCHB sampler themselves." call stop_all(t_r, "Wrong input to GAS-CI PCHB.") end select case default call Stop_All(t_r, trim(w)//" not a valid keyword") end select case ("GUGA-PCHB") ! use an explicit guga-pchb keyword and flag t_guga_pchb = .true. case ("UEG") ! Use the new UEG excitation generator. ! TODO: This probably isn't the best way to do this tUEGNewGenerator = .true. tLatticeGens = .true. case default call Stop_All("ReadSysInp", trim(w)//" not a valid keyword") end select non_uniform_rand_excits end do case ("PCHB-WEIGHTED-SINGLES") ! Enable using weighted single excitations with the pchb excitation generator write(stdout, *) trim(w)//" is deprecated." write(stdout, *) "For SD basis please use `nonuniformrandexcits pchb & &singles on-fly-heat-bath` instead." write(stdout, *) "For GUGA please use `GUGA-PCHB-WEIGHTED-SINGLES`." call stop_all(this_routine, trim(w)//" is deprecated.") case ("GUGA-PCHB-WEIGHTED-SINGLES") t_guga_pchb_weighted_singles = .true. ! enable intermediately some pchb+guga testing case("ANALYZE-PCHB") t_analyze_pchb = .true. case("OLD-GUGA-PCHB") ! original, rather unoptimized guga-pchb implementation ! still testing planned for final verdict on best guga-pchb ! implementation! t_old_pchb = .true. case("EXCHANGE-GUGA-PCHB") ! additional guga-pchb implementation for exchange type ! contributions. testing for final 'best' guga-pchb ! version still needs to be done t_exchange_pchb = .true. case("SPAWNLISTDETS") !This option will mean that a file called SpawnOnlyDets will be read in, ! and only these determinants will be allowed to be spawned at. CALL Stop_All("ReadSysInp", "SPAWNLISTDETS option deprecated") case ("UMATEPSILON") ! For systems that read in from an FCIDUMP file, any two-electron ! integrals are screened against a threshold parameter. Below ! this they are ignored. ! ! By default, this parameter is 10-e8, but it can be changed here. UMatEps = to_realdp(tokens%next()) case ("LMATEPSILON") ! Six-index integrals are screened, too, with the default being 1e-10 LMatEps = to_realdp(tokens%next()) case ("NOSINGEXCITS") !This will mean that no single excitations are ever attempted to be generated. tNoSingExcits = .true. case ("DIAGONALTMAT") !Implies that the orbital basis are eigenfunctions of the KE operator, so TMAT can be stored as a diagonal matrix to save space. tOneElecDiag = .true. !One electron integrals diagonal case ("ROHF") !This is an option for open-shell systems to specify that the integrals are *restricted* open-shell integrals. !This will save memory (around a factor of 16) for the integral storage, ! but the FCIDUMP file should be the same as before (ie in UHF form). tROHF = .true. tNoBrillouin = .true. tBrillouinsDefault = .false. IF (tFindCINatOrbs) CALL Stop_All("ReadSysInp", "For orbital rotations of open shell "& &//"systems, UMAT must be stored in spin & & orbitals - cannot be compressed using ROHF.") case ("LZTOT") tFixLz = .true. LzTot = to_int(tokens%next()) case ("KPOINTS") tKPntSym = .true. case ("IMPURITY-EXCITGEN") ! Use the impurity model excitation generator (star geometry) t_impurity_excitgen = .true. case ("MOLPROMIMIC") !Mimic the run-time behaviour of molpros NECI implementation tMolpro = .true. tMolproMimic = .true. case ("READ_ROFCIDUMP") call stop_all(t_r, 'Deprecated function. Use FCIDUMP-NAME ROFCIDUMP instead.') case ("FCIDUMP-NAME") FCIDUMP_name = tokens%next() case ("COMPLEXORBS_REALINTS") !We have complex orbitals, but real integrals. This means that we only have 4x permutational symmetry, !so we need to check the (momentum) symmetry before we look up any integrals tComplexOrbs_RealInts = .true. case ("COMPLEXWALKERS-REALINTS") ! In case complex walkers shall be used but not complex basis functions, ! such that the integrals are real and have full symmetry tComplexWalkers_RealInts = .true. t_complex_ints = .false. case ("SYSTEM-REPLICAS") ! How many copies of the simulation do we want to run in parallel? ! This can only be done using mneci.x, where the size of the ! representation (i.e. lenof_sign) is permitted to vary at runtime #ifdef PROG_NUMRUNS_ inum_runs = to_int(tokens%next()) tMultiReplicas = .true. #ifdef CMPLX_ lenof_sign = 2 * inum_runs #else lenof_sign = inum_runs #endif if (inum_runs > inum_runs_max) then write(stdout, *) 'Maximum SYSTEM-REPLICAS: ', inum_runs_max call stop_all(t_r, 'SYSTEM-REPLICAS is greater than maximum & &permitted value') end if #else itmp = to_int(tokens%next()) #ifdef DOUBLERUN_ if (itmp /= 2) then #else if (itmp /= 1) then #endif call stop_all(t_r, "mneci.x build must be used for running & &with multiple simultaneous replicas") end if #endif case("ADJOINT-REPLICAS") ! some replicas will be evolved according to the adjoint H, ! useful to get the left eigenvector in ST-FCIQMC #if defined(PROG_NUMRUNS_) || defined(DOUBLERUN_) t_adjoint_replicas = .true. #else call stop_all(this_routine, & "mneci or dneci necessary for 'adjoint-replicas'") #endif case ("HEISENBERG") tHeisenberg = .true. case("PERMUTE-ORBS") ! Apply a permutation of the orbital indices to the ! ordering given in the FCIDUMP file - only has an ! effect when reading an FCIDUMP file, has no effect for ! hubbard/heisenberg/ueg systems etc n_orb = 0 do while (tokens%remaining_items() > 0) n_orb = n_orb + 1 buf(n_orb) = to_int(tokens%next()) end do call load_orb_perm(buf(1:n_orb)) case ("GAS-SPEC") tGAS = .true. block logical :: recoupling integer :: nGAS, iGAS, i_sg, n_sg integer :: i_orb ! n_orbs are the number of spatial orbitals per GAS space ! cn_min, cn_max are cumulated particle numbers per GAS space integer, allocatable :: n_orbs(:), cn_min(:), cn_max(:), & spat_GAS_orbs(:), beta_orbs(:), supergroups(:, :) w = to_upper(tokens%next()) if (w == 'LOCAL') then allocate(LocalGASSpec_t :: GAS_specification) else if (w == 'CUMULATIVE') then allocate(CumulGASSpec_t :: GAS_specification) else if (w == 'FLEXIBLE') then allocate(FlexibleGASSpec_t :: GAS_specification) else call stop_all(t_r, 'You may pass either LOCAL or CUMULATIVE constraints.') end if select type(GAS_specification) type is (FlexibleGASSpec_t) nGAS = to_int(tokens%next()) n_sg = to_int(tokens%next()) allocate(supergroups(nGAS, n_sg), source=0) do i_sg = 1, n_sg if (file_reader%nextline(tokens, skip_empty=.false.)) then do iGAS = 1, nGAS supergroups(iGAS, i_sg) = to_int(tokens%next()) end do else call stop_all(t_r, 'Error in GAS spec.') end if end do class default nGAS = to_int(tokens%next()) allocate(n_orbs(nGAS), cn_min(nGAS), cn_max(nGAS), source=0) do iGAS = 1, nGAS if (file_reader%nextline(tokens, skip_empty=.false.)) then n_orbs(iGAS) = to_int(tokens%next()) cn_min(iGAS) = to_int(tokens%next()) cn_max(iGAS) = to_int(tokens%next()) else call stop_all(t_r, 'Error in GAS spec.') end if end do end select if (file_reader%nextline(tokens, skip_empty=.false.)) then call read_spat_GAS_orbs(tokens, spat_GAS_orbs, recoupling) else call stop_all(t_r, 'Error in GAS spec.') end if select type(GAS_specification) type is(CumulGASSpec_t) GAS_specification = CumulGASSpec_t(cn_min, cn_max, spat_GAS_orbs, recoupling) type is(LocalGASSpec_t) GAS_specification = LocalGASSpec_t(cn_min, cn_max, spat_GAS_orbs, recoupling) type is(FlexibleGASSpec_t) GAS_specification = FlexibleGASSpec_t(supergroups, spat_GAS_orbs, recoupling) call GAS_specification%write_to(stdout) class default call stop_all(t_r, "Invalid type for GAS specification.") end select end block case("OUTPUT-GAS-HILBERT-SPACE-SIZE") t_output_GAS_sizes = .true. case ("SD-SPIN-PURIFICATION") allocate(SD_spin_purification, source=possible_purification_methods%FULL_S2) spin_pure_J = to_realdp(tokens%next()) if (tokens%remaining_items() > 0) then w = to_upper(tokens%next()) select case(w) case ("ONLY-LADDER-OPERATOR") SD_spin_purification = possible_purification_methods%ONLY_LADDER case ("TRUNCATE-LADDER-OPERATOR") SD_spin_purification = possible_purification_methods%TRUNCATED_LADDER case default call stop_all(t_r, "Wrong alternative purification method.") end select end if case ("ENDSYS") exit system case default call stop_all(this_routine, "Keyword " //trim(w)//" not recognized in SYSTEM block") end select end do system if (t_lattice_model) then if (t_k_space_hubbard) then lat => lattice(lattice_type, length_x, length_y, length_z, & .not. t_open_bc_x,.not. t_open_bc_y,.not. t_open_bc_z, 'k-space') else if (t_new_real_space_hubbard) then lat => lattice(lattice_type, length_x, length_y, length_z, & .not. t_open_bc_x,.not. t_open_bc_y, & .not. t_open_bc_z, 'real-space', t_bipartite_order = t_bipartite_order) else lat => lattice(lattice_type, length_x, length_y, length_z, & .not. t_open_bc_x,.not. t_open_bc_y,.not. t_open_bc_z, & t_bipartite_order = t_bipartite_order) end if end if if (NEL == 0) then call stop_all(this_routine, "Number of electrons cannot be zero.") end if if (.not. tUEG2) then if (THUB .OR. TUEG .OR. .NOT. (TREADINT .OR. TCPMD .or. tVASP)) then if (NMAXX == 0) then call stop_all(this_routine, "Must specify CELL - & &the number of basis functions in each dim.") end if if (.NOT. THUB .AND. near_zero(BOX)) then call stop_all(this_routine, "Must specify BOX size.") end if if (TTILT .AND. .NOT. THUB) then call stop_all(this_routine, "TILT can only be specified with HUBBARD.") end if end if end if END SUBROUTINE SysReadInput 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 Subroutine SysCleanup() use sym_mod, only: EndSym CALL ENDSYM() End Subroutine SysCleanup logical function AreSameSpatialOrb(i, j) ! Test whether spin orbitals i and j are from the same spatial orbital. ! Returns true if i and j are the *same* spin orbital or are the alpha ! and beta spin orbitals of the same spatial orbital. integer :: i, j integer :: a, b AreSameSpatialOrb = .false. if (i == j) then AreSameSpatialOrb = .true. else a = min(i, j) b = max(i, j) if (G1(a)%Ms == -1 .and. b - a == 1) then ! a is the alpha and b is the beta of the same spatial orbital. AreSameSpatialOrb = .true. end if end if end function AreSameSpatialOrb !================================================================ ! UEG2 Subroutines !================================================================ SUBROUTINE LatticeInit(RS, FKF) ! initiates the reciprocal lattice, real volume and the Fermi vector real(dp), intent(out) :: RS, FKF ! check dimension if (dimen == 3) then ! 3D OMEGA = 4.0_dp / 3.0_dp * PI * FUEGRS**3 * NEL RS = (3.0_dp * OMEGA / (4.0_dp * PI * NEL))**THIRD FKF = (9 * PI / 4)**THIRD / RS ! define lattice vectors and lattice constant in reciprocal space if (real_lattice_type == "sc") then k_lattice_constant = 2.0_dp * PI / OMEGA**THIRD k_lattice_vectors(1, 1:3) = (/1, 0, 0/) k_lattice_vectors(2, 1:3) = (/0, 1, 0/) k_lattice_vectors(3, 1:3) = (/0, 0, 1/) else if (real_lattice_type == "bcc") then k_lattice_constant = 2.0_dp * PI / (2.0_dp * OMEGA)**THIRD k_lattice_vectors(1, 1:3) = (/0, 1, 1/) k_lattice_vectors(2, 1:3) = (/1, 0, 1/) k_lattice_vectors(3, 1:3) = (/1, 1, 0/) else if (real_lattice_type == "fcc") then k_lattice_constant = 2.0_dp * PI / (4.0_dp * OMEGA)**THIRD k_lattice_vectors(1, 1:3) = (/-1, 1, 1/) k_lattice_vectors(2, 1:3) = (/1, -1, 1/) k_lattice_vectors(3, 1:3) = (/1, 1, -1/) else write(stdout, '(A)') 'lattice type not valid' end if else if (dimen == 2) then !2D write(stdout, '(A)') ' NMAXZ=0 : 2D calculation' OMEGA = PI * FUEGRS**2 * NEL RS = (OMEGA / (PI * NEL))**(1.0_dp / 2.0_dp) FKF = sqrt(2.0_dp) / RS ! define lattice vectors and lattice constant in reciprocal space k_lattice_constant = 2.0_dp * PI / OMEGA**(1.0_dp / 2.0_dp) k_lattice_vectors(1, 1:3) = (/1, 0, 0/) k_lattice_vectors(2, 1:3) = (/0, 1, 0/) k_lattice_vectors(3, 1:3) = (/0, 0, 0/) else if (dimen == 1) then !1D write(stdout, '(A)') ' NMAXZ=0, NMAXY=0 : 1D calculation' OMEGA = 2.0_dp * FUEGRS * NEL RS = OMEGA / (2.0_dp * NEL) FKF = (PI / 2.0_dp) / RS !for spin polarised simulation ! define lattice vectors and lattice constant in reciprocal space k_lattice_constant = 2.0_dp * PI / OMEGA k_lattice_vectors(1, 1:3) = (/1, 0, 0/) k_lattice_vectors(2, 1:3) = (/0, 0, 0/) k_lattice_vectors(3, 1:3) = (/0, 0, 0/) else write(stdout, '(A)') 'Problem with dimension! ' end if return END SUBROUTINE LatticeInit SUBROUTINE CalcCell !Detemines the cell size for a given cutoff and lattice type integer :: ii, jj, kk, EE logical :: under_cutoff if (real_lattice_type == "sc" .OR. dimen < 3) then NMAXX = int(sqrt(orbEcutoff)) + 1 if (dimen > 1) NMAXY = int(sqrt(orbEcutoff)) + 1 if (dimen > 2) NMAXZ = int(sqrt(orbEcutoff)) + 1 else if (real_lattice_type == "bcc" .or. real_lattice_type == "fcc") then ! calculate needed cell size ii = 0 ! ii is always positiv. jj varies from -ii to ii, kk from -|jj| to |jj| under_cutoff = .true. do while (ii <= int(orbEcutoff) .and. under_cutoff) !until no E < cutoff was found under_cutoff = .false. jj = -ii do while (abs(jj) <= abs(ii) .and. .not. under_cutoff) !until E < cutoff is found or jj =ii kk = -abs(jj) do while (abs(kk) <= abs(jj) .and. .not. under_cutoff)!until E < cutoff is found or kk=jj !calculate unscaled energy for ii, jj, kk EE = (k_lattice_vectors(1, 1) * ii + k_lattice_vectors(2, 1) * jj + k_lattice_vectors(3, 1) * kk)**2 EE = EE + (k_lattice_vectors(1, 2) * ii + k_lattice_vectors(2, 2) * jj + k_lattice_vectors(3, 2) * kk)**2 EE = EE + (k_lattice_vectors(1, 3) * ii + k_lattice_vectors(2, 3) * jj + k_lattice_vectors(3, 3) * kk)**2 if (EE <= orbEcutoff) under_cutoff = .true. kk = kk + 1 end do jj = jj + 1 end do ii = ii + 1 end do NMAXX = ii!-1 NMAXY = ii!-1 NMAXZ = ii!-1 end if ! lattice type return END SUBROUTINE CalcCell SUBROUTINE CalcTau !Detemines tau for a given lattice type if (dimen == 3) then ! 3D ! Hij_min**-1 call assign_value_to_tau((k_lattice_constant**2 * OMEGA) / (4.0_dp * PI), & 'Initialization in System_neci.') if (tTruncInitiator) call assign_value_to_tau(TAU * InitiatorWalkNo, 'Initialization in System_neci.') if (tHPHF) call assign_value_to_tau(TAU / sqrt(2.0_dp), 'Initialization in System_neci.') call assign_value_to_tau(0.9_dp * TAU * 4.0_dp / (NEL * (NEL - 1)) / (NBASIS - NEL), & 'Initialization in System_neci.') if (TAU > k_lattice_constant**(-2) / OrbEcutoff) then ! using Hii call assign_value_to_tau(1.0_dp / (k_lattice_constant**(2) * OrbEcutoff), & 'Initialization in System_neci.') write(stdout, *) '***************** Tau set by using Hii *******************************' !write(stdout,*) 1.0_dp/((2.0_dp*PI/Omega**third)**2*orbEcutoff) else write(stdout, *) 'Tau set by using Hji' end if else if (dimen == 2) then !2D ! Hij_min**-1 call assign_value_to_tau((k_lattice_constant * OMEGA) / (2.0_dp * PI), & 'Initialization in System_neci.') if (tTruncInitiator) call assign_value_to_tau(TAU * InitiatorWalkNo, & 'Initialization in System_neci.') if (tHPHF) call assign_value_to_tau(TAU / sqrt(2.0_dp), 'Initialization in System_neci.') call assign_value_to_tau(0.9_dp * TAU * 4.0_dp / (NEL * (NEL - 1)) / (NBASIS - NEL), & 'Initialization in System_neci.') if (TAU > k_lattice_constant**(-2) / OrbEcutoff) then !!!!!!!! NOT WORKING YET!!!!!!! !using Hii call assign_value_to_tau(1.0_dp / (k_lattice_constant**(2) * OrbEcutoff), & 'Initialization in System_neci.') write(stdout, *) '***************** Tau set by using Hii *******************************' else write(stdout, *) 'Tau set by using Hji' end if else if (dimen == 1) then !1D call assign_value_to_tau(OMEGA / (-2.0_dp * log(1.0_dp / (2.0_dp * sqrt(orbEcutoff)))), & 'Initialization in System_neci.') if (tTruncInitiator) call assign_value_to_tau(TAU * InitiatorWalkNo, & 'Initialization in System_neci.') if (tHPHF) call assign_value_to_tau(TAU / sqrt(2.0_dp), & 'Initialization in System_neci.') call assign_value_to_tau(0.9_dp * TAU * 4.0_dp / (NEL * (NEL - 1)) / (NBASIS - NEL), & 'Initialization in System_neci.') if (TAU > 0.9_dp * 1.0_dp / (0.5_dp * (k_lattice_constant)**2 * NEL * OrbEcutoff)) then ! using Hii call assign_value_to_tau(0.9_dp * 1.0_dp / (0.5_dp * (k_lattice_constant)**2 * NEL * OrbEcutoff), & 'Initialization in System_neci.') write(stdout, *) '***************** Tau set by using Hii *******************************' else write(stdout, *) 'Tau set by using Hji' end if end if !dimension write(stdout, *) 'Tau set to: ', TAU return END SUBROUTINE CalcTau function calc_madelung() result(res) real(dp) :: kappa integer :: i1, i2, i3, i4 integer :: n2 real(dp) :: k2, ek2, recipsum2 real(dp) :: t1, modr, er2, realsum2 real(dp) :: term2, term4 integer :: r_lattice_vectors(3, 3) real(dp) :: r_lattice_constant integer:: kvecX, kvecY, kvecZ integer :: tvecX, tvecY, tvecZ real(dp) :: inacc_madelung, temp_sum real(dp) :: k2max integer :: jj real(dp) :: xx, yy, zz real(dp) :: res inacc_madelung = 1.0d-15 ! Cannot be set lower than 1.0d-15 !------------------------------- 3D lattice ----------------------------------------- kappa = 2.8_dp / OMEGA**(1.0_dp / 3.0_dp) if (real_lattice_type == "sc") then r_lattice_constant = OMEGA**THIRD r_lattice_vectors(1, 1:3) = (/1, 0, 0/) r_lattice_vectors(2, 1:3) = (/0, 1, 0/) r_lattice_vectors(3, 1:3) = (/0, 0, 1/) else if (real_lattice_type == "fcc") then r_lattice_constant = 0.5_dp * (2.0_dp * OMEGA)**THIRD r_lattice_vectors(1, 1:3) = (/-1, 1, 1/) r_lattice_vectors(2, 1:3) = (/1, -1, 1/) r_lattice_vectors(3, 1:3) = (/1, 1, -1/) else if (real_lattice_type == "bcc") then r_lattice_constant = 0.5_dp * (4.0_dp * OMEGA)**THIRD r_lattice_vectors(1, 1:3) = (/0, 1, 1/) r_lattice_vectors(2, 1:3) = (/1, 0, 1/) r_lattice_vectors(3, 1:3) = (/1, 1, 0/) else write(stdout, '(A)') 'lattice type not valid' end if term2 = -pi / (kappa**2.0_dp * OMEGA) ! write(stdout,*) term2, "term2" term4 = -2.0_dp * kappa / sqrt(pi) ! write(stdout,*) term4, "term4" recipsum2 = 0.0_dp temp_sum = 0.0_dp k2max = 0.0_dp i4 = 1 do while (temp_sum * (1.0_dp + inacc_madelung) <= recipsum2) temp_sum = recipsum2 recipsum2 = 0.0_dp do i1 = -i4, i4 do i2 = -i4, i4 do i3 = -i4, i4 kvecX = k_lattice_vectors(1, 1) * i1 + k_lattice_vectors(2, 1) * i2 + k_lattice_vectors(3, 1) * i3 kvecY = k_lattice_vectors(1, 2) * i1 + k_lattice_vectors(2, 2) * i2 + k_lattice_vectors(3, 2) * i3 kvecZ = k_lattice_vectors(1, 3) * i1 + k_lattice_vectors(2, 3) * i2 + k_lattice_vectors(3, 3) * i3 n2 = kvecX**2 + kvecY**2 + kvecZ**2 k2 = (k_lattice_constant / 2.0_dp / PI)**2.0_dp * real(n2, dp) ek2 = (1.0_dp / OMEGA) * (1.0_dp / (pi * k2)) * exp(-pi**2.0_dp * k2 / kappa**2.0_dp) if (n2 > k2max) k2max = n2 if (n2 /= 0) then ! write(stdout,*) k2,ek2 ! for testing recipsum2 = recipsum2 + ek2 end if end do end do end do i4 = i4 + 1 ! write(stdout,*) "i4 kmax^2", i4-1, k2max*k_lattice_constant**2 end do realsum2 = 0.0_dp temp_sum = 0.0_dp i4 = 1 do while (temp_sum * (1.0_dp + inacc_madelung) <= realsum2) temp_sum = realsum2 realsum2 = 0.0_dp do i1 = -i4, i4 do i2 = -i4, i4 do i3 = -i4, i4 tvecX = r_lattice_vectors(1, 1) * i1 + r_lattice_vectors(2, 1) * i2 + r_lattice_vectors(3, 1) * i3 tvecY = r_lattice_vectors(1, 2) * i1 + r_lattice_vectors(2, 2) * i2 + r_lattice_vectors(3, 2) * i3 tvecZ = r_lattice_vectors(1, 3) * i1 + r_lattice_vectors(2, 3) * i2 + r_lattice_vectors(3, 3) * i3 n2 = tvecX**2 + tvecY**2 + tvecZ**2 t1 = r_lattice_constant * sqrt(real(n2, dp)) if (.not. near_zero(t1)) then er2 = error_function_c(kappa * t1) / t1 realsum2 = realsum2 + er2 end if end do end do end do ! write(stdout,*) 'i4, realsum2' , i4, realsum2 i4 = i4 + 1 end do ! write(stdout,*) "real space", realsum2 !------------------------------- output ----------------------------------------- res = realsum2 + recipsum2 + term2 + term4 write(stdout, *) "Calculating Madelung Constant - Fraser et al. PRB 53 4 1814" write(stdout, *) "kappa taken from CASINO manual to be", kappa write(stdout, *) "k2_max", k2max, k2max * k_lattice_constant**2 ! write(stdout,*) "omega", OMEGA ! write(stdout, *) "klatticeconstant", k_lattice_constant ! write(stdout, *) "rlatticeconstant", r_lattice_constant write(stdout, *) "Madelung constant", res return end function subroutine read_spat_GAS_orbs(tokens, spat_GAS_orbs, recoupling) use fortran_strings, only: operator(.in.), Token_t, split, can_be_int use input_parser_mod, only: TokenIterator_t type(TokenIterator_t), intent(inout) :: tokens integer, allocatable, intent(out) :: spat_GAS_orbs(:) logical, intent(out) :: recoupling integer :: times, iGAS, i type(Token_t), allocatable :: splitted(:) type(buffer_int_1D_t) :: buffer character(len=100) :: w routine_name("read_spat_GAS_orbs") call buffer%init() recoupling = .true. do while (tokens%remaining_items() > 0) w = to_upper(tokens%next()) if ('*' .in. w) then splitted = split(w, '*') times = to_int(splitted(1)%str) iGAS = to_int(splitted(2)%str) else if (can_be_int(w)) then times = 1 iGAS = to_int(w) else if (w == 'RECOUPLING') then recoupling = .true. exit else if (w == 'NO-RECOUPLING') then recoupling = .false. exit else call stop_all(this_routine, "Error in reading GAS orbitals.") end if do i = 1, times call buffer%push_back(iGAS) end do end do call buffer%dump_reset(spat_GAS_orbs) end subroutine END MODULE System