subroutine orthogonalise_replica_pairs(iter_data)
! complex walkers not supported here. This routine is never called: deprecated code?
! Apply a Gram Schmidt orthogonalisation to the different system
! replicas.
! This routine works in cases where there are two replicas representing
! each excited state. Whenever an overlap has to be calculated between
! different excited states, both pairs are used and averaged. When a
! norm needs to be calculated, the overlap is taken between the two
! replicas for the same excited state.
! Underscores denote the excited state label and uperscores represent
! the replica label (i.e. 1 or 2):
! |psi_2^1'> = |psi_2^1> - (|psi_1^1><psi_1^2|psi_2^1> +
! |psi_1^2><psi_1^1|psi_2^1>)/(2<psi_1^1|psi_1^2>)
type(fcimc_iter_data), intent(inout) :: iter_data
integer :: tgt_state, src_state, run, j, irep, imod1, imod2, TotWalkersNew
real(dp) :: norms(inum_runs.div.2), overlaps(inum_runs, inum_runs)
real(dp) :: sgn(lenof_sign), sgn_orig(2), delta, r
logical :: tCoreDet
character(len=*), parameter :: this_routine = "orthogonalise_replica_pairs"
ASSERT(inum_runs == lenof_sign)
#ifndef PROG_NUMRUNS_
unused_var(iter_data)
call stop_all(this_routine, "orthogonalise replicas requires mneci.x")
#else
norms = 0.0_dp
overlaps = 0.0_dp
! Loop over all excited states (thus over all replica *pairs*).
do tgt_state = 1, inum_runs / 2
HolesInList = 0
do j = 1, int(TotWalkers)
! n.b. We are using a non-contiguous list (Hash algorithm).
call extract_sign(CurrentDets(:, j), sgn)
tCoreDet = check_determ_flag(CurrentDets(:, j))
if (IsUnoccDet(sgn) .and. (.not. tCoreDet)) then
HolesInList = HolesInList + 1
cycle
end if
if (tgt_state > 1) then
! The two replica signs for this excited state.
sgn_orig(1) = sgn(tgt_state * 2 - 1)
sgn_orig(2) = sgn(tgt_state * 2)
! Loop over both replicas for this excited state.
do irep = 1, 2
! These variables equal (1,2) and (2,1) when irep = 1 and
! 2, repsectively. They are used to access the correct
! indices of the sign arrays.
imod1 = mod(irep, 2)
imod2 = 1 - imod1
do src_state = 1, tgt_state - 1
delta = -sgn(src_state * 2 - imod1) * &
all_overlaps(src_state * 2 - imod2, tgt_state * 2 - imod1) / &
(2 * all_norms(src_state))
sgn(tgt_state * 2 - imod1) = sgn(tgt_state * 2 - imod1) + delta
delta = -sgn(src_state * 2 - imod2) * &
all_overlaps(src_state * 2 - imod1, tgt_state * 2 - imod1) / &
(2 * all_norms(src_state))
sgn(tgt_state * 2 - imod1) = sgn(tgt_state * 2 - imod1) + delta
end do
end do
call encode_sign(CurrentDets(:, j), sgn)
! Note that we shouldn't be able to kill all particles on a
! site, as we can only change a run if there are particles in
! a lower indexed run to change...
! The exception is when using semi-stochastic, where
! unoccupied determinants can be stored.
if (.not. tSemiStochastic) then
ASSERT(.not. IsUnoccDet(sgn))
end if
! Now we need to our accounting to make sure that NoBorn/
! Died/etc. counters remain reasonable.
!
! n.b. we don't worry about the delta=0 case, as adding 0 to
! the accumulators doesn't cause errors...
do irep = 1, 2
run = tgt_state * 2 - mod(irep, 2)
if (sgn(run) >= 0 .neqv. sgn_orig(irep) >= 0) then
NoDied(run) = NoDied(run) + abs(sgn_orig(irep))
iter_data%ndied(run) = iter_data%ndied(run) &
+ abs(sgn_orig(irep))
NoBorn(run) = NoBorn(run) + abs(sgn(run))
iter_data%nborn(run) = iter_data%nborn(run) &
+ abs(sgn(run))
else if (abs(sgn(run)) >= abs(sgn_orig(irep))) then
NoBorn(run) = NoBorn(run) &
+ abs(sgn(run) - sgn_orig(irep))
iter_data%nborn(run) = iter_data%nborn(run) &
+ abs(sgn(run) - sgn_orig(irep))
else
NoDied(run) = NoDied(run) &
+ abs(sgn(run) - sgn_orig(irep))
iter_data%ndied(run) = iter_data%ndied(run) &
+ abs(sgn(run) - sgn_orig(irep))
end if
end do
end if
! We should be accumulating the norm for this run, so it can
! be used in processing later runs. We should also be
! accumulating the overlap terms for doing the same
!
! n.b. These are the values _after_ orthogonalisation with
! the previous runs
norms(tgt_state) = norms(tgt_state) + sgn(2 * tgt_state - 1) * sgn(2 * tgt_state)
! Calculate overlaps for all 4 combinations of the 2 replicas
! of the 2 excited states.
do run = tgt_state + 1, inum_runs / 2
overlaps(tgt_state * 2 - 1, run * 2 - 1) = overlaps(tgt_state * 2 - 1, run * 2 - 1) &
+ sgn(tgt_state * 2 - 1) * sgn(run * 2 - 1)
overlaps(tgt_state * 2 - 1, run * 2) = overlaps(tgt_state * 2 - 1, run * 2) &
+ sgn(tgt_state * 2 - 1) * sgn(run * 2)
overlaps(tgt_state * 2, run * 2 - 1) = overlaps(tgt_state * 2, run * 2 - 1) &
+ sgn(tgt_state * 2) * sgn(run * 2 - 1)
overlaps(tgt_state * 2, run * 2) = overlaps(tgt_state * 2, run * 2) &
+ sgn(tgt_state * 2) * sgn(run * 2)
end do
end do
! And ensure that the norm/overlap data is accumulated onto all
! of the processors.
call MPISumAll(norms, all_norms)
call MPISumAll(overlaps, all_overlaps)
end do
! We now need to aggregate statistics here, rather than at the end
! of annihilation, as they have been modified by this routine...
!
! n.b. TotWalkersNew is integer, TotWalkers is int64 to ensure that it
! can be collected using AllTotWalkers (which may end up with far
! too many walkers).
TotWalkersNew = int(TotWalkers)
call CalcHashTableStats(TotWalkersNew, iter_data)
TotWalkers = TotWalkersNew
#endif
end subroutine orthogonalise_replica_pairs