anchor.diff:                  if (all(qs%vary_shift)) then
                   do i = 1, qs%psip_list%nstates
                       state = qs%psip_list%states(:,i) 
  -                    call binary_search_i0_list_trot(time_avg_psip_list%states, state, 1, time_avg_psip_list%nstates, hit, pos)
  +                    call binary_search_i0_list_trot(time_avg_psip_list_states, state, 1, nstates_sq, hit, pos)
                       !call binary_search(time_avg_psip_list%states, state, 1, time_avg_psip_list%nstates, hit, pos)
                       if (hit) then
  -                          time_avg_psip_list%nparticles = time_avg_psip_list%nparticles - abs(real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor)
  -                          time_avg_psip_list%pops(:,pos) = time_avg_psip_list%pops(:,pos) + qs%psip_list%pops(:,i) 
  -                          time_avg_psip_list%nparticles = time_avg_psip_list%nparticles + abs(real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor)
  +                          !time_avg_psip_list%nparticles = time_avg_psip_list%nparticles - abs(real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor)
  +                          time_avg_psip_list_pops(pos) = time_avg_psip_list_pops(pos) + (real(qs%psip_list%pops(1,i))/qs%psip_list%pop_real_factor)
  +                          !time_avg_psip_list%nparticles = time_avg_psip_list%nparticles + abs(real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor)
                             time_avg_psip_list_sq(2,pos) = time_avg_psip_list_sq(2,pos) + (real(qs%psip_list%pops(1,i))/qs%psip_list%pop_real_factor)**2 
                          else
  -                           do j = time_avg_psip_list%nstates, pos, -1
  +                           do j = nstates_sq, pos, -1
                                  ! i is the number of determinants that will be inserted below j.
                                  k = j + 1 
  -                               time_avg_psip_list%states(:,k) = time_avg_psip_list%states(:,j)
  -                               time_avg_psip_list%pops(:,k) = time_avg_psip_list%pops(:,j)
  -                               time_avg_psip_list%dat(:,k) = time_avg_psip_list%dat(:,j)
  +                               time_avg_psip_list_states(:,k) = time_avg_psip_list_states(:,j)
  +                               time_avg_psip_list_pops(k) = time_avg_psip_list_pops(j)
  +                               !time_avg_psip_list%dat(:,k) = time_avg_psip_list%dat(:,j)
                                  time_avg_psip_list_sq(1,k) = time_avg_psip_list_sq(1,j)
                                  time_avg_psip_list_sq(2,k) = time_avg_psip_list_sq(2,j)
                              end do
  -                           call insert_new_walker(sys, time_avg_psip_list, annihilation_flags, pos, &
  -                               state, qs%psip_list%pops(:,i), qs%ref)
  -                           time_avg_psip_list%nstates = time_avg_psip_list%nstates+1
  +                           !call insert_new_walker(sys, time_avg_psip_list, annihilation_flags, pos, &
  +                           !    state, qs%psip_list%pops(:,i), qs%ref)
  +                           time_avg_psip_list_states(:,pos) = qs%psip_list%states(:,i)
  +                           time_avg_psip_list_pops(pos) = (real(qs%psip_list%pops(1,i))/qs%psip_list%pop_real_factor)
  +                           !time_avg_psip_list%nstates = time_avg_psip_list%nstates+1
                              time_avg_psip_list_sq(1,pos) = qs%psip_list%states(1,i)
                              time_avg_psip_list_sq(2,pos) = (real(qs%psip_list%pops(1,i))/qs%psip_list%pop_real_factor)**2
                              nstates_sq = nstates_sq + 1
                              ! Extract the real sign from the encoded sign.
  -                           real_population = real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor
  -                           time_avg_psip_list%nparticles = time_avg_psip_list%nparticles + abs(real_population)
  +                           !real_population = real(time_avg_psip_list%pops(:,pos))/time_avg_psip_list%pop_real_factor
  +                           !time_avg_psip_list%nparticles = time_avg_psip_list%nparticles + abs(real_population)
                          end if
                   end do
                   ! Update time average.
  -                time_avg_psip_list%nparticles = time_avg_psip_list%nparticles/(iter-avg_start+1)
  -                time_avg_psip_list%pops(:,:time_avg_psip_list%nstates) =  time_avg_psip_list%pops(:,:time_avg_psip_list%nstates)/(-avg_start+iter+1)
  +                !time_avg_psip_list%nparticles = time_avg_psip_list%nparticles/(iter-avg_start+1)
  +                time_avg_psip_list_pops(:nstates_sq) =  time_avg_psip_list_pops(:nstates_sq)/(-avg_start+iter+1)
   
                   time_avg_psip_list_sq(2,:nstates_sq) =  time_avg_psip_list_sq(2,:nstates_sq)/(-avg_start+iter+1)
                   end if
  
positive.code: 
          if (parent) then
              write (io_unit, '(1X, "Time-averaged cluster populations",/)')
              do i = 1, time_avg_psip_list%nstates
                  call write_qmc_var(io_unit, time_avg_psip_list%states(1,i))
                  call write_qmc_var(io_unit, time_avg_psip_list%pops(1,i))
                  write (io_unit,'()')
              end do
              write (io_unit, '(1X, "Time-averaged cluster populations squared",/)')
              do i = 1, time_avg_psip_list%nstates
                  call write_qmc_var(io_unit, int(time_avg_psip_list_sq(1,i)))
                  call write_qmc_var(io_unit, time_avg_psip_list_sq(2,i))
                  write (io_unit,'()')
  
negative.code:         real(p), allocatable :: real_population(:), var_energy(:)
          logical :: old_vary
          integer :: avg_start
          real(p) :: p_ref
  
          old_vary = .false.
          avg_start = 0
  