anchor.diff:  
     end function compute_force_harmonic_trap
   
  +  subroutine rescale_at_walls
  +
  +    integer :: i
  +    integer :: cell_idx, n, start, wall_idx, cell(3)
  +    double precision :: local_v(3), local_k, local_T, factor
  +
  +    do cell_idx = 1, solvent_cells% N
  +       if (solvent_cells% cell_count(cell_idx) <= 1) cycle
  +
  +       start = solvent_cells% cell_start(cell_idx)
  +       n = solvent_cells% cell_count(cell_idx)
  +
  +       ! Find whether we are in a wall cell
  +       cell = compact_h_to_p(cell_idx - 1, solvent_cells% M) + 1
  +       if (cell(3) == 1) then
  +          wall_idx = 1
  +       else if (cell(3) == solvent_cells% L(3)) then
  +          wall_idx = 2
  +          local_T = 1.1
  +       else
  +          wall_idx = -1
  +       end if
  +       if (wall_idx==-1) cycle
  +       local_T = wall_t(wall_idx)
  +
  +       local_v = 0
  +       do i = start, start + n - 1
  +          local_v = local_v + solvent% vel(:, i)
  +       end do
  +       local_v = local_v / n
  +
  +       local_k = 0
  +       do i = start, start + n - 1
  +          local_k = local_k + sum((solvent% vel(:, i)-local_v)**2)/2
  +       end do
  +       local_k = local_k/(3*(n-1))
  +       factor = sqrt(local_T/(2*local_k))
  +
  +       do i = start, start + n - 1
  +          solvent% vel(:, i) = local_v + factor*(solvent% vel(:, i)-local_v)
  +       end do
  +
  +    end do
  +
  +  end subroutine rescale_at_walls
  +
   end program setup_sphere_thermo_trap
  
positive.code:   N_MD_steps = PTread_i(config, 'N_MD')
    dt = tau / N_MD_steps
    N_loop = PTread_i(config, 'N_loop')
  
    wall_t = PTread_dvec(config, 'wall_T', 2)
    T = PTread_d(config, 'T')
  
negative.code:      solvent%pos_old = solvent% pos
       colloids%pos_old = colloids% pos
  
       call rescale_at_walls
       call wall_mpcd_step(solvent, solvent_cells, state, &
            wall_temperature=wall_t, wall_v=wall_v, wall_n=[rho, rho])
  
  