subroutine real_time_read_input(file_reader)
class(FileReader_t), intent(inout) :: file_reader
type(TokenIterator_t) :: tokens
logical :: eof
character(100) :: w
character(*), parameter :: this_routine = "real_time_read_input"
integer, parameter :: lesser = -1, greater = 1
integer :: i, j
integer, allocatable :: buffer(:)
! set the flag that this is a real time calculation
t_real_time_fciqmc = .true.
! usually only real-valued FCIDUMPs
t_complex_ints = .false.
#ifndef CMPLX_
call stop_all(this_routine, "Real-time calculations require kneci or kmneci")
#endif
! and set default values for the real-time calculation
call set_real_time_defaults()
real_time: do while (file_reader%nextline(tokens, skip_empty=.true.))
w = to_upper(tokens%next())
select case (w)
! have to enter all the different input options here
case ("VERLET")
! using a verlet algorithm instead of the second order runge-kutta
tVerletScheme = .true.
if (tokens%remaining_items() > 0) iterInit = to_int(tokens%next())
if (stepsAlpha == 1) write(stdout, *) "Warning: STEPSALPHA is 1. Ignoring VERLET keyword"
case ("DAMPING")
! to reduce the explosive spread of walkers through the
! Hilbert space a small imaginery energy can be introduced in
! the Schroedinger equation id/dt y(t) = (H-E0-ie)y(t)
real_time_info%damping = to_realdp(tokens%next())
case ("ROTATE-TIME")
! If the time is to be rotated by some angle time_angle to increase
! stability, this can be set here
t_rotated_time = .true.
real_time_info%time_angle = to_realdp(tokens%next())
! use nicks perturbation & kp-fciqmc stuff here as much as
! possible too
case ("PROJECT-INITIAL-STATE")
! If we specify this, we do not create the overlap state specifically,
! but copy the currentdets.
! This can greatly help to overcome memory problems, as it basically
! halves the memory required for initial state preparation
tNewOverlap = .false.
! just compute the time-evolution of a singly excited state (with
! reference to the ground state. This gives the contribution of
! this state to the spectrum
case ("SINGLE")
alloc_popsfile_dets = .true.
! deprecated, replace by MULTI
tWritePopsNorm = .true.
! Now, overlap state and initial state are the same
tNewOverlap = .false.
allocate(pops_pert(1))
pops_pert%nannihilate = 1
pops_pert%ncreate = 1
allocate(pops_pert(1)%crtn_orbs(1))
allocate(pops_pert(1)%ann_orbs(1))
pops_pert(1)%ann_orbs(1) = to_int(tokens%next())
pops_pert(1)%crtn_orbs(1) = to_int(tokens%next())
call init_perturbation_annihilation(pops_pert(1))
call init_perturbation_creation(pops_pert(1))
case ("KSPACE")
! Apply the perturbations in kspace. This only does something for real
! space hubbard, else it is the default
t_kspace_operators = .true.
! Arbitrary perturbation on the initial state, always get the overlap
! with the initial state
case ("MULTI")
alloc_popsfile_dets = .true.
tWritePopsNorm = .true.
allocate(pops_pert(1))
allocate(overlap_pert(1))
allocate(buffer(nel))
j = 0
! first, read all orbitals to which particles shall be added
do
if (tokens%remaining_items() == 0) exit
i = to_int(tokens%next())
! -1 is the terminator for creation and indicates that all following
! orbitals are to be annihilated
if (i == -1) exit
j = j + 1
! as the size of pops_pert is unknown, use a buffer
buffer(j) = i
end do
! allocate creation operators
if (j > 0) then
pops_pert%ncreate = j
allocate(pops_pert(1)%crtn_orbs(j))
! we take the overlap with the initial state, so overlap_pert == pops_pert
overlap_pert%ncreate = j
allocate(overlap_pert(1)%crtn_orbs(j))
do i = 1, j
pops_pert(1)%crtn_orbs(i) = buffer(i)
overlap_pert(1)%crtn_orbs(i) = pops_pert(1)%crtn_orbs(i)
end do
end if
j = 0
! now, read in all orbitals from which particles shall be removed
do while (tokens%remaining_items() > 0)
j = j + 1
buffer(j) = to_int(tokens%next())
end do
! again, allocate annihilation operators
if (j > 0) then
pops_pert%nannihilate = j
overlap_pert%nannihilate = j
allocate(pops_pert(1)%ann_orbs(j))
allocate(overlap_pert(1)%ann_orbs(j))
do i = 1, j
pops_pert(1)%ann_orbs(i) = buffer(i)
overlap_pert(1)%ann_orbs(i) = pops_pert(1)%ann_orbs(i)
end do
end if
call init_perturbation_annihilation(pops_pert(1))
call init_perturbation_creation(pops_pert(1))
call init_perturbation_annihilation(overlap_pert(1))
call init_perturbation_creation(overlap_pert(1))
deallocate(buffer)
! the most important info is if it is the photoemmission(lesser GF)
! or photoabsorption (greater GF) and the orbital we want the
! corresponding operator apply on
! the type of GF considered also changes the sign of the FT exponent
! decision for now: input a specific GF matrix element and the type
! of the greensfunction to be calculated(lesser,greater) eg:
! lesser i j : <y(0)| a^+_i a_j |y(0)>
case ("LESSER")
alloc_popsfile_dets = .true.
! lesser GF -> photo emission: apply a annihilation operator
tOverlapPert = .true.
tWritePopsNorm = .true.
! i probably also can use the overlap-perturbed routines
! from nick
! but since applying <y(0)|a^+_i for all i is way cheaper
! and should be done for all possible and allowed i.
! and creating all those vectors should be done in the init
! step and stored, and then just calc. the overlap each time
! step
! store the information of the type of greensfunction
gf_type = lesser
! probably have to loop over spin-orbitals dont i? yes!
! if no specific orbital is specified-> loop over all j!
! but only do that later: input is a SPINORBITAL!
if (tokens%remaining_items() > 0) then
allocate(pops_pert(1))
pops_pert%nannihilate = 1
allocate(pops_pert(1)%ann_orbs(1))
pops_pert(1)%ann_orbs(1) = to_int(tokens%next())
call init_perturbation_annihilation(pops_pert(1))
else
call stop_all(this_routine, "Invalid input for Green's function")
end if
if (tokens%size() == 3) then
gf_count = 1
!allocate the perturbation object
! and also the lefthand perturbation object for overlap
allocate(overlap_pert(1))
overlap_pert%nannihilate = 1
allocate(overlap_pert(1)%ann_orbs(1))
! read left hand operator first
overlap_pert(1)%ann_orbs(1) = to_int(tokens%next())
call init_perturbation_annihilation(overlap_pert(1))
! If the created and annihilated orbital are the same, we
! do not need to explicitly construct the projection state,
! this might save a lot of memory
if (pops_pert(1)%ann_orbs(1) == overlap_pert(1)%ann_orbs(1)) &
tNewOverlap = .false.
else if (tokens%size() == 2) then
allGfs = 1
else
call stop_all(this_routine, "Invalid input for Green's function")
end if
case ("GREATER")
alloc_popsfile_dets = .true.
! greater GF -> photo absorption: apply a creation operator
tOverlapPert = .true.
tWritePopsNorm = .true.
! i probably also can use the overlap-perturbed routines
! from nick
! but since applying <y(0)|a_i for all i is way cheaper
! and should be done for all possible and allowed i.
! and creating all those vectors should be done in the init
! step and stored, and then just calc. the overlap each time
! step
! store type of greensfunction
gf_type = greater
! if no specific orbital is specified-> loop over all j!
! but only do that later
if (tokens%remaining_items() > 0) then
allocate(pops_pert(1))
pops_pert%ncreate = 1
allocate(pops_pert(1)%crtn_orbs(1))
pops_pert(1)%crtn_orbs(1) = to_int(tokens%next())
call init_perturbation_creation(pops_pert(1))
else
call stop_all(this_routine, "Invalid input for Green's function")
end if
if (tokens%size() == 3) then
! allocate the perturbation object
allocate(overlap_pert(1))
overlap_pert%ncreate = 1
allocate(overlap_pert(1)%crtn_orbs(1))
overlap_pert(1)%crtn_orbs(1) = to_int(tokens%next())
call init_perturbation_creation(overlap_pert(1))
! If the created and annihilated orbital are the same, we
! do not need to explicitly construct the projection state,
! this might save a lot of memory
if (pops_pert(1)%crtn_orbs(1) == overlap_pert(1)%crtn_orbs(1)) &
tNewOverlap = .false.
else if (tokens%size() == 2) then
allGfs = 2
else
call stop_all(this_routine, "Invalid input for Green's function")
end if
case ("SCALE-POPULATION")
tScalePopulation = .true.
case ("LOWER-THRESHOLD")
! indicates that the given rotation-threshold is not an upper
! but in fact a lower threshold, so the variation is switched
! on as soon as the walker number drops below
! this is a particularly useless thing in most cases, but for
! proving some stuff, it saves the day
tLowerThreshold = .true.
case ("FULLY-ROTATED")
! for testing purposes, it is useful to do pure imaginary
! time evolution with the rotated time algorithm -> this is
! enabled by this keyword
! in addition, this disables the usage of input POPSFILEs for
! more efficient ground state search (the real-time POPSFILE
! read-in settings are not useful for ground state search)
t_rotated_time = .true.
tWalkContGrow = .true.
real_time_info%time_angle = 2 * atan(1.0_dp)
case ("PRINT-POP")
! include the time-dependent population of targeted orbitals into
! the output. This requires them to be evaluated on the fly
numSnapShotOrbs = 0
allocate(buffer(tokens%size() + 1))
do
if (tokens%remaining_items() > 0) then
numSnapShotOrbs = numSnapShotOrbs + 1
! nBasis is not defined at this point, so we cannot check if
! there are too many items given - no serious input will contain
! more arguments than basis states anyway
buffer(numSnapShotOrbs) = to_int(tokens%next())
else
exit
end if
end do
allocate(snapShotOrbs(numSnapShotOrbs))
snapShotOrbs(1:numSnapShotOrbs) = buffer(1:numSnapShotOrbs)
deallocate(buffer)
case ("NOSHIFT")
! disabling the shift gives higher precision results as no
! renormalization of the norm by a dynamic factor is made
! note that the walker number will grow exponentially in this
! scenario, however
asymptoticShift = 0.0_dp
! might want to set DiagSft = 0.0_dp but there migth also be some cases
! in which this is unwanted
tStaticShift = .true.
case ("START-HF")
! do not read in an initial state from a POPSFILE and apply a perturbation
! but start right away in the HF as the initial state does not matter in
! principle for the spectrum
tReadPops = .false.
tStartSinglePart = .true.
case ("STABILIZE-WALKERS")
! enabling this activates the dynamic shift as soon as the walker number drops
! below 80% of the peak value
tStabilizerShift = .true.
if (tokens%remaining_items() > 0) then
asymptoticShift = to_realdp(tokens%next())
tStaticShift = .true.
end if
case ("UNCONSTRAINED-SHIFT")
! use an unconstrained shift mode that also allows
! negative shifts
tOnlyPositiveShift = .false.
write(stdout, *) &
"WARNING: Using an unconstrained shift can lead to instabilities"
case ("HF-OVERLAP")
! take the overlap not with the initial state but with the perturbed
! reference
tHFOverlap = .true.
case ("ENERGY-BENCHMARK")
! one can specify an energy which shall be added as a global shift
! to the hamiltonian. Useful for getting transition energies
benchmarkEnergy = to_realdp(tokens%next())
case ("DYNAMIC-CORE")
tDynamicCoreSpace = .true.
! if dynamic core is set, the core space for semistochastic treatment is
! updated every few hundred iterations according to the currently most
! occupied determinants
case ("COMPLEX-INTEGRALS")
! in the real-time implementation, since we need the complex
! functionality anyway, we have to additionally tell the
! program, that the FCIDUMP input is complex
! the default is that they are real!
t_complex_ints = .true.
case ("NSPAWNMAX")
! specify a maximum number of spawn attempts per determinant in
! regulation mode (i.e. for large number of spawns)
nspawnMax = to_int(tokens%next())
case ("COMPLEXWALKERS-COMPLEXINTS")
! if we really use complex integrals, we have to tell as the
! default is using real integrals with complex walkers
tComplexWalkers_RealInts = .false.
case ("RT-POPS")
! in addition to the 'normal' popsfile, a second one is supplied
! containing a time evolved state
tRealTimePopsfile = .true.
case ("OVERPOPULATE")
! enabling sets the options for time-dependent shift and rotation
! such that a positive shift will occur with a stable walker number
t_rotated_time = .true.
tDynamicAlpha = .true.
! this is done by pinning the shift to some positive value and
! then auto-adjusting the rotation
tStaticShift = .true.
! here, rotation and shift variation have to start at the same point
! (in principle, it is only required that the rotation does not start
! before shift variation) to prevent the rotation from converging on
! its own, circumventing the overpopulation via positive shift
tOverpopulate = .true.
! it is most efficient to turn on the shift after equilibration of the angle
! so this is done via the stabilize-walkers feature
tStabilizerShift = .true.
if (tokens%remaining_items() > 0) then
asymptoticShift = to_realdp(tokens%next())
else
asymptoticShift = 2.0_dp
end if
case ("DYNAMIC-ROTATION")
! this automatically adjusts the temporal rotation to find a minimal
! alpha guaranteeing a fixed walker number
tDynamicAlpha = .true.
t_rotated_time = .true.
if (tokens%remaining_items() > 0) alphaDamping = to_realdp(tokens%next())
case ("ROTATION-THRESHOLD")
! number of walkers at which the variation of rotation angle starts
! 0 by default
rotThresh = to_int(tokens%next())
case ("STEPSALPHA")
! length of the decay channel update cycle (in timesteps)
! i.e. angle of rotation and damping
stepsAlpha = to_int(tokens%next())
if (stepsAlpha == 1 .and. tVerletScheme) write(stdout, *) &
"Warning: STEPSALPHA is 1. Ignoring VERLET keyword"
case ("DYNAMIC-DAMPING")
! allow the damping to be time-dependent
! optional: damping parameter for the adjustment of eta
tDynamicDamping = .true.
if (tokens%remaining_items() > 0) etaDamping = to_realdp(tokens%next())
case ("LIMIT-SHIFT")
! limits the shift to some maximum value. On short times, the threshold
! can be exceeded.
tLimitShift = .true.
! optional argument: threshold value (absolute value!). Default is 3
if (tokens%remaining_items() > 0) shiftLimit = to_realdp(tokens%next())
case ("INFINITE-INIT")
! use the initiator adaptiation without any inititators - works well
! in some real-time applications
! this is not equivalent to switching on initiators without the
! addtoinitiator keyword as infinite-init will also remove all
! existing inititators
!
! Note that this option requires `core-inits OFF`
tInfInit = .true.
tAddtoInitiator = .true.
tTruncInitiator = .true.
case ("LOG-TRAJECTORY")
! This prints out the complex time trajectory in the form of alpha(iter)
! and tau(iter)
tLogTrajectory = .true.
case("QUAD-DAMP")
! Additional energy-dependent damping (quadratic in H)
if (tokens%remaining_items() > 0) then
real_time_info%quad_damp_fac = to_realdp(tokens%next())
else
real_time_info%quad_damp_fac = 0.5d0
end if
case ("GENERATE-CORESPACE")
! Now, we write out the most important determinants along the contour
! Also, the contour is logged
tGenerateCoreSpace = .true.
tLogTrajectory = .true.
! optionally, we can supply the number of states to log
ss_space_in%tpops = .true.
if (tokens%remaining_items() > 0) then
ss_space_in%npops = to_int(tokens%next())
else
ss_space_in%npops = 1000
end if
if (tSemiStochastic) call stop_all(this_routine, &
"GENERATE-CORESPACE NOT AVAILABLE IN SEMI-STOCHASTIC MODE")
case ("CORESPACE-THRESHOLD")
! Set the threshold from which on a determinant is in the corespace
wn_threshold = to_realdp(tokens%next())
case ("CORESPACE-LOG-INTERVAL")
! Set the number of iterations after which we get the new candidates for the
! corespace
corespace_log_interval = to_int(tokens%next())
case ("READ-TRAJECTORY")
! This reads in a trajectory and performs the time-evolution along
! it
tReadTrajectory = .true.
case ("LIVE-TRAJECTORY")
! Now we re-read the trajectory during runtime, this can be used to
! use a trajectory that is currently being determined
tReadTrajectory = .true.
tLiveTrajectory = .true.
case ("CORESPACE-OVERLAP")
! Get the Green's function for the corespace only. This performs the
! time-evolution only in the semistochastic space.
tGZero = .true.
! If the corespace-greensfunction is to be obtained, semi-stochastic
! has to be turned on
if (.not. tSemiStochastic) call stop_all(this_routine, &
"CORESPACE-OVERLAP ONLY AVAILABLE IN SEMI-STOCHASTIC MODE")
case ("ENDREALTIME")
exit real_time
case default
call stop_all(this_routine, "Keyword "//trim(w)//" not recognized in REALTIME block")
end select
end do real_time
end subroutine real_time_read_input