real_time_read_input.F90 Source File


Contents


Source Code

module real_time_read_input_module
  use real_time_data
  use real_time_init, only: set_real_time_defaults, benchmarkenergy
  use FciMCData, only: alloc_popsfile_dets, pops_pert
  use CalcData, only: tAddToInitiator, tTruncInitiator, tWalkContGrow, tStartSinglePart, &
       tWritePopsNorm, tReadPops, ss_space_in, tSemiStochastic
  use perturbations, only: init_perturbation_creation, init_perturbation_annihilation
  use kp_fciqmc_data_mod, only: tOverlapPert, overlap_pert, tScalePopulation
  use SystemData, only: nel, tComplexWalkers_RealInts, t_complex_ints
  use constants
  use input_parser_mod, only: FileReader_t, TokenIterator_t
  use fortran_strings, only: to_upper, to_lower, to_int, to_realdp
  use util_mod, only: stop_all

  implicit none

  contains

    ! need a real-time calc read_input routine to seperate that as much
    ! from the rest of the code as possible!
    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

end module real_time_read_input_module