subroutine orthogonalise_replicas(iter_data)
! Apply a Gram Schmidt orthogonalisation to the different system
! replicas.
!
! |psi_2'> = |psi_2> - (|psi_1><psi_1|psi_2>)/(<psi_1|psi_1>)
type(fcimc_iter_data), intent(inout) :: iter_data
integer :: tgt_run, src_run, run, j, TotWalkersNew
integer :: tgt_sgn, src_sgn
real(dp) :: all_overlaps_real(inum_runs, inum_runs)
integer :: low_loop_bound, high_loop_bound, loop_step
real(dp) :: norms(inum_runs), overlaps_real(inum_runs, inum_runs)
real(dp) :: sgn(lenof_sign), sgn_orig, delta_real, r
#ifdef CMPLX_
! the min_part_type and max_part_type macros return the index of the sgns
! corresponding to the real and imag parts of a replica respectively
real(dp) :: overlaps_imag(inum_runs, inum_runs), all_overlaps_imag(inum_runs, inum_runs)
real(dp) :: delta_imag
#endif
logical :: tCoreDet
character(*), parameter :: this_routine = 'orthogonalise_replicas'
#ifndef PROG_NUMRUNS_
call stop_all(this_routine, "orthogonalise replicas requires mneci or kmneci")
#endif
! If we are using the symmetric orthogonaliser, then bypass this
! routine...
if (tOrthogonaliseSymmetric) then
call orthogonalise_replicas_lowdin(iter_data)
return
end if
norms = 0.0_dp
overlaps_real = 0.0_dp
#ifdef CMPLX_
overlaps_imag = 0.0_dp
#endif
do tgt_run = 1, inum_runs
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
! Loop over source runs, and subtract out components
sgn_orig = sgn(tgt_run)
! set loop parameters
if (tPairedReplicas) then
low_loop_bound = 2 - mod(tgt_run, 2)
high_loop_bound = tgt_run - 2
loop_step = 2
else
low_loop_bound = 1
high_loop_bound = tgt_run - 1
loop_step = 1
end if
do src_run = low_loop_bound, high_loop_bound, loop_step
#ifdef CMPLX_
! (a + ib)*(c + id) = ac - bd + i(ad +bc)
delta_real = -(sgn(min_part_type(src_run)) &
* all_overlaps_real(src_run, tgt_run) &
- sgn(max_part_type(src_run)) &
* all_overlaps_imag(src_run, tgt_run)) &
/ all_norms(src_run)
sgn(min_part_type(tgt_run)) = sgn(min_part_type(tgt_run)) + delta_real
delta_imag = -(sgn(min_part_type(src_run)) &
* all_overlaps_imag(src_run, tgt_run) &
+ sgn(max_part_type(src_run)) &
* all_overlaps_real(src_run, tgt_run)) &
/ all_norms(src_run)
sgn(max_part_type(tgt_run)) = sgn(max_part_type(tgt_run)) + delta_imag
#else
if (.not. near_zero(all_norms(src_run))) then
delta_real = -sgn(src_run) * all_overlaps_real(src_run, tgt_run) &
/ all_norms(src_run)
end if
! test if a small remaining overlap still gives the
! correct shift..
if (t_test_overlap .and. tgt_run == inum_runs) then
! draw a random sign and change by the
! specified epsilon
! or choose just a random number out of [-eps,+eps]
! or try a prefactor..
delta_real = overlap_eps * delta_real
if (n_stop_ortho > 0 .and. iter > n_stop_ortho) then
delta_real = 0.0_dp
end if
end if
sgn(tgt_run) = sgn(tgt_run) + delta_real
#endif
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...
if (sgn(tgt_run) >= 0 .neqv. sgn_orig >= 0) then
NoDied(tgt_run) = NoDied(tgt_run) + abs(sgn_orig)
iter_data%ndied(tgt_run) = iter_data%ndied(tgt_run) &
+ abs(sgn_orig)
NoBorn(tgt_run) = NoBorn(tgt_run) + abs(sgn(tgt_run))
iter_data%nborn(tgt_run) = iter_data%nborn(tgt_run) &
+ abs(sgn(tgt_run))
else if (abs(sgn(tgt_run)) >= abs(sgn_orig)) then
NoBorn(tgt_run) = NoBorn(tgt_run) &
+ abs(sgn(tgt_run) - sgn_orig)
iter_data%nborn(tgt_run) = iter_data%nborn(tgt_run) &
+ abs(sgn(tgt_run) - sgn_orig)
else
NoDied(tgt_run) = NoDied(tgt_run) &
+ abs(sgn(tgt_run) - sgn_orig)
iter_data%ndied(tgt_run) = iter_data%ndied(tgt_run) &
+ abs(sgn(tgt_run) - sgn_orig)
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
#ifdef CMPLX_
! (a+ib)*(a-ib) = a^2 + b^2
norms(tgt_run) = norms(tgt_run) &
+ sgn(min_part_type(tgt_run))**2 + sgn(max_part_type(tgt_run))**2
#else
norms(tgt_run) = norms(tgt_run) + sgn(tgt_run)**2
#endif
do run = tgt_run + 1, inum_runs
#ifdef CMPLX_
! we want the inner product of the vectors:
! ___
! <n| = \ ` [ a(n,r) - ib(n,r) ] <r|
! /__,
! r
! ___
! |m> = \ ` [ a(n,r) + ib(n,r) ] |r>
! /__,
! r
! ___
! <n|m> = \ ` a(n,r)*a(m,r) + b(n,r)*b(m,r)
! /__, + i [ a(n,r)*b(m,r) - b(n,r)*a(m,r) ]
! r
!
! Note that symmetry may no longer be assumed for the overlap matrix
! i.e. overlaps_real(n,m) = overlaps_real(m,n), but
! overlaps_imag(n,m) = -overlaps_imag(m,n)
! only storing overlaps for run(m) > tgt_run(n):
! <n|m> =
! _ _
! | |
! | x |
! | x x |
! | x x x |
! | x x x x |
! | x x x x x |
! |_x x x x x x _|
overlaps_real(tgt_run, run) = overlaps_real(tgt_run, run) &
+ sgn(min_part_type(run)) * sgn(min_part_type(tgt_run)) &
+ sgn(max_part_type(run)) * sgn(max_part_type(tgt_run))
overlaps_imag(tgt_run, run) = overlaps_imag(tgt_run, run) &
+ sgn(min_part_type(run)) * sgn(max_part_type(tgt_run)) &
+ sgn(max_part_type(run)) * sgn(min_part_type(tgt_run))
overlaps_real(run, tgt_run) = overlaps_real(tgt_run, run)
overlaps_imag(run, tgt_run) = -overlaps_imag(tgt_run, run)
#else
overlaps_real(tgt_run, run) = overlaps_real(tgt_run, run) &
+ sgn(run) * sgn(tgt_run)
overlaps_real(run, tgt_run) = overlaps_real(tgt_run, run)
#endif
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_real, all_overlaps_real)
#ifdef CMPLX_
call MPISumAll(overlaps_imag, all_overlaps_imag)
#endif
end do
! Store a normalised overlap matrix for each of the replicas.
do src_run = 1, inum_runs - 1
do tgt_run = src_run + 1, inum_runs
if (all_norms(src_run) * all_norms(tgt_run) > EPS) then
replica_overlaps_real(src_run, tgt_run) = &
all_overlaps_real(src_run, tgt_run) / &
sqrt(all_norms(src_run) * all_norms(tgt_run))
#ifdef CMPLX_
replica_overlaps_imag(src_run, tgt_run) = &
all_overlaps_imag(src_run, tgt_run) / &
sqrt(all_norms(src_run) * all_norms(tgt_run))
#endif
end if
end do
end do
! Store in global data for perturbation accumulation.
all_overlaps = all_overlaps_real
! 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
end subroutine orthogonalise_replicas