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