subroutine perform_real_time_fciqmc
! main real-time calculation routine
! do all the setup, read-in and calling of the "new" real-time MC loop
use real_time_data, only: gf_overlap
use FciMCData, only: TotImagTime
use real_time_procs, only: expand_corespace_buf
use fcimc_output, only: PrintHighPops
implicit none
character(*), parameter :: this_routine = "perform_real_time_fciqmc"
integer :: j, i, iterRK, err, all_err
real(dp) :: s_start, s_end, tstart(2), tend(2)
real(dp) :: totalTime
complex(dp), allocatable :: overlap_buf(:)
complex(dp), allocatable :: norm_buf(:)
character(255) :: rtPOPSFILE_name
logical :: t_comm_done
rtPOPSFILE_name = 'TIME_EVOLVED_POP'
write(stdout, *) " ========================================================== "
write(stdout, *) " ------------------ Real-time FCIQMC ---------------------- "
write(stdout, *) " ========================================================== "
! call the real-time setup routine and all the initialization
call init_real_time_calc_single()
! counts the number of iterations in Runge-Kutta when creating
! initial states for verlet
iterRK = 0
write(stdout, *) " Real-time FCIQMC initialized! "
! rewrite the major original neci core loop here and adapt it to
! the new necessary real-time stuff
! check nicks kp code, to have a guideline in how to go into that!
! as a first test if everything is set up correctly calculate the
! initial overlap for the same indiced of creation and annihilation
! operators <y(0)|a^+_i a_i |y(0)> = n_i
! and <y(0)|a_i a^+_i |y(0)> = 1 - n_i
! if normed correctly
! in the quantity perturbed_ground the left hand side and in
! CurrentDets the right hand side should be stored ..
! should write a routine which calulated the overlap of CurrentDets
! with the perturbed groundstate and stores it in a pre-allocated list
! shouldnt geourge have some of those routines..
! normsize and gf_count are initialized in init_real_time_calc_single
! -> cannot be used before
allocate(overlap_buf(gf_count), stat=i)
allocate(norm_buf(normsize), stat=i)
call update_gf_overlap()
write(stdout, *) "test on overlap at t = 0: "
if (.not. tRealTimePopsfile) then
if (gf_type == -1) then
write(stdout, *) " for lesser GF <y(0)| a^+_i a_j |y(0) >; i,j: ", &
overlap_pert(1)%ann_orbs(1), pops_pert(1)%ann_orbs(1)
else if (gf_type == 1) then
write(stdout, *) " for greater GF <y(0)| a_i a^+_j |y(0)> ; i,j: ", &
overlap_pert(1)%crtn_orbs(1), pops_pert(1)%crtn_orbs(1)
end if
end if
write(stdout, *) "Current GF:", gf_overlap(1, 1) / pert_norm(1, 1), pert_norm(1, 1), normsize
write(stdout, *) "Normalization", pert_norm(1, 1), dyn_norm_red(1, 1)
! enter the main real-time fciqmc loop here
fciqmc_loop: do while (.true.)
t_comm_done = .false.
! the timing stuff has to be done a bit differently in the
! real-time fciqmc, since there are 2 spawing and annihilation
! steps involved...
call set_timer(walker_time)
if (iProcIndex == root) s_start = neci_etime(tstart)
! the iter data is used in updating the shift, so it has to be reset before
! output
! this is a bad implementation : iter should be local
! if the trajectory is logged, print alpha and tau here
if (tLogTrajectory) call logTimeCurve()
if (StepsPrint > 0 .and. .not. tCoupleCycleOutput) then
if (mod(iter, StepsPrint) == 0) then
call iteration_output_wrapper(iter_data_fciqmc, TotParts, tPairedReplicas)
! mark that the communication has been done
t_comm_done = .true.
end if
end if
! update the overlap each time
! rmneci_setup: computation of instantaneous projected norm is shifted to here
if (mod(iter, StepsSft) == 0) then
! current overlap is now the one after iteration
! update the normalization due to the shift
norm_buf = calc_norm(CurrentDets, int(TotWalkers))
call MPIReduce(norm_buf, MPI_SUM, dyn_norm_psi)
call update_gf_overlap()
do j = 1, gf_count
do i = 1, normsize
current_overlap(i, j) = gf_overlap(i, j) / dyn_norm_red(i, j) * &
exp(shift_damping(((i - 1) / inum_runs + 1)))
end do
!normalize the greens function
! averaging should probably be done over the normalized
! GFs
overlap_buf(j) = sum(gf_overlap(:, j)) / sum(dyn_norm_red(:, j)) * &
sum(exp(shift_damping)) / inum_runs
overlap_real(j) = real(overlap_buf(j))
overlap_imag(j) = aimag(overlap_buf(j))
end do
call update_real_time_iteration(.not. t_comm_done)
end if
if (mod(iter, stepsAlpha) == 0 .and. (tDynamicAlpha .or. tDynamicDamping)) then
call adjust_decay_channels()
! the verlet scheme only works for constant time step, hence, upon adjusting
! the time step, we need to do another iteration of RK2
if (tVerletScheme) call end_verlet_sweep()
end if
if (tVerletScheme) call check_verlet_sweep(iterRK)
! if a threshold value is set, check it
if (tLimitShift) call trunc_shift()
if (tGenerateCoreSpace .and. (mod(iter, corespace_log_interval) == 0)) then
call expand_corespace_buf(core_space_buf, csbuf_size)
write(stdout, *) "New corespace buffer size", csbuf_size
end if
! perform the actual iteration(excitation generation etc.)
if (iterRK == 0) iter = iter + 1
if (tVerletScheme .and. .not. tVerletSweep) iterRK = iterRK + 1
if (tVerletSweep) then
call perform_verlet_iteration(err)
else
call perform_real_time_iteration(err)
end if
! exit the calculation if something failed
call MPISumAll(err, all_err)
if (all_err /= 0) then
! Print an error message, then exit
write(stdout, *) "Error occured during real-time iteration"
exit
end if
! check if somthing happpened to stop the iteration or something
call check_real_time_iteration()
! If required, set up the references
if (Iter == allDoubsInitsDelay + 1 .and. nRefs > 1 .and. tDelayGetRefs) &
call setup_reference_space(.true.)
if (iProcIndex == root) then
s_end = neci_etime(tend)
totalTime = real(s_end - s_global_start, dp)
end if
call MPIBcast(totalTime)
if (tTimeExit .and. totalTime > MaxTimeExit) then
! reset the maximum number of iterations to the index of the next one
nmcyc = Iter + StepsSft
! to prevent nmcyc from being updated
tTimeExit = .false.
end if
! load balancing
if (tLoadBalanceBlocks .and. mod(iter, 1000) == 0 .and. (.not. tSemiStochastic)) then
call adjust_load_balance(iter_data_fciqmc)
end if
if (mod(iter, 400) == 0 .and. tSemiStochastic .and. tDynamicCoreSpace) &
call refresh_semistochastic_space()
if (iProcIndex == root) then
s_end = neci_etime(tend)
IterTime = IterTime + (s_end - s_start)
end if
! update the normalization of the greensfunction according to damping (dynamic)
call update_shift_damping()
if (tSoftExitFound) exit fciqmc_loop
! we have to stop here since the gf_overlap array is now full
if (iter == nmcyc) exit fciqmc_loop
end do fciqmc_loop
call PrintHighPops()
if (tPopsFile) then
! as both the elapsed real-time and the elapsed imaginary time are required
! for dynamic alpha, both are stored. That is, the total elapsed real time
! is now stored instead of tau. If no tau is supplied upon read-in of the
! generated popsfile, it is estimated using the number of cycles, the current
! angle of rotation and the total elapsed real time
TotImagTime = elapsedImagTime
call assign_value_to_tau(elapsedRealTime, this_routine)
! THIS IS A HACK: We dont want to alter the POPSFILE functions themselves
! so we sneak in the shift_damping into some slot unimportant to rneci
AllSumNoatHF(1:inum_runs) = shift_damping
! Another hack: we also sneak the value of alpha somewhere here, the shift
! is not meaningful in real-time anyway
DiagSft = real_time_info%time_angle
call WriteToPopsfileParOneArr(CurrentDets, TotWalkers, rtPOPSFILE_name)
end if
! For testing purpose, report the value of the green's function at the end
if (gf_count > 0) &
write(stdout, *) "Final real part of Green function", overlap_real(1)
! We finish writing any extra output files like the tauContour and
! the CORESPACE file
if (tLogTrajectory) call closeTauContourFile()
if (tWriteCoreEnd) call write_most_pop_core_at_end(write_end_core_size)
! GENERATE-CORESPACE has precendence over WRITE-CORE-END
if (tGenerateCoreSpace) call get_corespace_from_buf(core_space_buf, csbuf_size)
deallocate(norm_buf, stat=i)
deallocate(overlap_buf, stat=i)
! rt_iter_adapt : avoid memleaks
call dealloc_real_time_memory()
end subroutine perform_real_time_fciqmc