anchor.diff:              end if
               ! Only one cluster of this size to choose => p_clust = 1
           case default
  -            ! Select cluster from the excitors on the current processor in
  -            ! a uniform manner.
  -            ! Not allowed to select the reference as it is not an excitor.
  -            ! First excitor 'seeds' the cluster:
  -            pos = int(genrand_real2()*(tot_walkers-1)) + 1 ! TODO: adjust for parallel
  -            if (pos >= D0_pos) pos = pos + 1
  -            cdet%f = walker_dets(:,pos)
  -            cluster_population = walker_population(1,pos)
  -            cluster%excitors(1)%f => walker_dets(:,pos)
  -            if (walker_population(1,pos) <= initiator_population) then
  -                cdet%initiator_flag = 1
  -            end if
  -            ! Now the rest...
  -            do i = 2, cluster%nexcitors
  -                ! Select a position between in the excitors list.
  -                pos = int(genrand_real2()*(tot_walkers-1)) + 1 ! TODO: adjust for parallel
  -                if (pos >= D0_pos) pos = pos + 1
  -                cluster%excitors(i)%f => walker_dets(:,pos)
  -                call collapse_cluster(walker_dets(:,pos), walker_population(1,pos), cdet%f, cluster_population)
  -                if (walker_population(1,pos) <= initiator_population) then
  -                    cdet%initiator_flag = 1
  +            ! Select cluster from the excitors on the current processor with
  +            ! probability for choosing an excitor proportional to the excip
  +            ! population on that excitor.
  +            do i = 1, cluster%nexcitors
  +                ! Select a position in the excitors list.
  +                pop = int(genrand_real2()*tot_excip_pop) + 1 ! TODO: adjust for parallel
  +                call binary_search(cumulative_excip_pop, pop, 1, tot_walkers, hit, pos)
  +                ! Not allowed to select the reference as it is not an excitor.
  +                ! Because we treat (for the purposes of the cumulative
  +                ! population) the reference to have 0 excips, then
  +                ! cumulative_excip_pop(D0_pos) = cumulative_excip_pop(D0_pos-1).
  +                ! The binary search algorithm assumes each value in the array
  +                ! being searched is unique, which is not true, so we can
  +                ! accidentally find D0_pos.  As we want to find pos such that
  +                ! cumulative_excip_pop(pos-1) < pop <= cumulative_excip_pop(pos),
  +                ! then this means we actually need the slot before D0_pos.
  +                ! Correcting for this accident is much easier than producing an
  +                ! array explicitly without D0...
  +                if (pos == D0_pos) pos = pos - 1
  +                if (i == 1) then
  +                    ! First excitor 'seeds' the cluster:
  +                    cdet%f = walker_dets(:,pos)
  +                    cluster_population = walker_population(1,pos)
  +                else
  +                    call collapse_cluster(walker_dets(:,pos), walker_population(1,pos), cdet%f, cluster_population)
                   end if
                   if (cluster_population == 0) exit
  +                if (walker_population(1,pos) <= initiator_population) cdet%initiator_flag = 1
  +                ! Probability of choosing this excitor = pop/tot_pop.
  +                cluster%pselect = (cluster%pselect*abs(walker_population(1,pos)))/tot_excip_pop
  +                cluster%excitors(i)%f => walker_dets(:,pos)
               end do
   
  -            ! We chose excitors uniformly.
  -            !  -> prob. of each excitor is 1./(number of excitors on processor).
  -            ! Overall probability is the product of this times the number of
  -            ! ways the cluster could have been selected, ie
  -            !   n_s!/n_excitors^n_s
  -            ! This is because we consider the cluster t_j t_i and t_i t_j to be
  -            ! the same cluster and hence must account for the number of ways the
  -            ! specific unique cluster can be generated.
  -            cluster%pselect = (cluster%pselect*factorial(cluster%nexcitors))/(real(tot_walkers-1,p)**cluster%nexcitors)
  +            ! We chose excitors with a probability proportional to their
  +            ! occupation.  However, because (for example) the cluster t_X t_Y
  +            ! and t_Y t_X collapse onto the same excitor (where X and Y each
  +            ! label an excitor), the probability of selecting a given cluster is
  +            ! proportional to the number of ways the cluster could have been
  +            ! formed.  (One can view this factorial contribution as the
  +            ! factorial prefactors in the series expansion of e^T.)
  +            cluster%pselect = cluster%pselect*factorial(cluster%nexcitors)
   
               if (cluster_population /= 0) cluster%excitation_level = get_excitation_level(f0, cdet%f)
               ! To contribute the cluster must be within a double excitation of
  
positive.code: 
                  ! Normalisation factor for cluster%amplitudes...
                  cluster%amplitude = real(cluster_population,p)/(D0_normalisation**(cluster%nexcitors-1))
              else
                  ! Simply set excitation level to a too high (fake) level to avoid
                  ! this cluster being used.
  
negative.code:         ! truncation_level+2 excitors.
          ! Following the process described by Thom in 'Initiator Stochastic
          ! Coupled Cluster Theory' (unpublished), each size, n_s, has probability
          ! p(n_s) = 1/2^(n_s+1), n_s=0,truncation_level and p(truncation_level+2)
          ! is such that \sum_{n_s=0}^{truncation_level+2} p(n_s) = 1.
          rand = genrand_real2()
          psize = 0.0_p
          cluster%nexcitors = -1
          do i = 0, max_size-1
              psize = psize + 1.0_p/2**(i+1)
              if (rand < psize) then
                  ! Found size!
                  cluster%nexcitors = i
                  cluster%pselect = cluster%pselect/2**(i+1)
                  exit
              end if
          end do
          ! If not set, then must be the largest possible cluster
          if (cluster%nexcitors == -1) then
              cluster%nexcitors = max_size
              cluster%pselect = cluster%pselect*(1.0_p - psize)
          end if
  
          ! Initiator approximation: no point using a CAS so just use the populations.
          ! This is sufficiently quick that we'll just do it in all cases, even
  