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