anchor.diff:              qs%estimators%proj_energy_old = get_sanitized_projected_energy(qs)
               if (all(qs%vary_shift) .and. not(old_vary)) then
                   avg_start = iter
  -                time_avg_psip_list%pops(:,:) = qs%psip_list%pops(:,:)
  -                time_avg_psip_list%states(:,:) = qs%psip_list%states(:,:)
  -                time_avg_psip_list%nparticles = sum(abs(qs%psip_list%pops(:,:)))/qs%psip_list%pop_real_factor
  -                time_avg_psip_list%nstates = qs%psip_list%nstates
  +                time_avg_psip_list_pops(:qs%psip_list%nstates) = real(qs%psip_list%pops(1,:qs%psip_list%nstates))/qs%psip_list%pop_real_factor
  +                time_avg_psip_list_states(:,:qs%psip_list%nstates) = qs%psip_list%states(:,:qs%psip_list%nstates)
                   time_avg_psip_list_sq(1,:) = qs%psip_list%states(1,:)
                   time_avg_psip_list_sq(2,:) = (real(qs%psip_list%pops(1,:))/qs%psip_list%pop_real_factor)**2
                   nstates_sq = qs%psip_list%nstates
  
positive.code:                 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(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_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
                                 ! 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_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
                             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)
                         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_sq(2,:nstates_sq) =  time_avg_psip_list_sq(2,:nstates_sq)/(-avg_start+iter+1)
                  end if
  
negative.code:                 end if
  
                  if(all(qs%vary_shift)) then
                      !time_avg_psip_list%nparticles = time_avg_psip_list%nparticles*(iter-avg_start)
                      time_avg_psip_list_pops(:nstates_sq) =  time_avg_psip_list_pops(:nstates_sq)*(iter - avg_start)
  
                      time_avg_psip_list_sq(2,:nstates_sq) =  time_avg_psip_list_sq(2,:nstates_sq)*(iter-avg_start)
                  end if
  