anchor.diff:            bulk_temperature=set_temperature)
        v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
   
  -     call md_vel(solvent, tau)
   
        T = compute_temperature(solvent, solvent_cells, tz)
        call elem_v_com%append(v_com, i, i*tau)
  
positive.code:      call wall_mpcd_step(solvent, solvent_cells, state, &
            wall_temperature=wall_t, wall_v=wall_v, wall_n=[rho, rho], thermostat=thermostat, &
            bulk_temperature=set_temperature)
       call md_vel(solvent, tau)
    end do
  
    do i = 1, N_loop
  
negative.code:         write(*,*) i 
       end if
       call mpcd_stream_xforce_yzwall(solvent, solvent_cells, tau, gravity_field(1))
       call md_vel(solvent, tau)
       call apply_pbc(solvent, solvent_cells%edges)
       call random_number(solvent_cells% origin)
       solvent_cells% origin = solvent_cells% origin - 1
  