anchor.diff:          ! 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,max_level-1 and p(max_level)
  -        ! is such that \sum_{n_s=0}^{max_level} p(n_s) = 1, and hence
  -        ! p(max_level) = 1/2^max_level.
  +        ! 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_sum = 0.0_p
  +        psize = 0.0_p
           cluster%nexcitors = -1
  -        do i = 0, max_cluster_size -1
  -            psize = 1.0_p/2**(i+1)
  -            psize_sum = psize_sum + psize
  -            if (rand < psize_sum) then
  +        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 reached end of loop, then must be the largest possible cluster.
  -        if (i == max_cluster_size) then
  -            cluster%nexcitors = max_cluster_size
  -            psize = 1.0_p/2**max_cluster_size
  +        ! 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
  -        cluster%pselect = cluster%pselect*psize
   
           ! 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
  
positive.code: 
          logical :: soft_exit
  
          integer :: D0_pos, max_cluster_size
          logical :: hit
  
          real :: t1, t2
  
negative.code: 
                  ! Normalisation factor for cluster%amplitudes...
                  cluster%amplitude = real(cluster_population,p)/(D0_normalisation**(cluster%nexcitors-1))
                  ! Factorial due to the series expansion of exp(\hat{T}).
                  cluster%amplitude = cluster%amplitude!/factorial(cluster%nexcitors)
              else
                  ! Simply set excitation level to a too high (fake) level to avoid
                  ! this cluster being used.
  