We investigate the dependence of the shear viscosity of suspensions of spherical colloids as a function of the volume fraction of the suspension, the colloid–colloid interactions and the shear rate. We couple molecular dynamics to describe the motion of the colloids with stochastic rotation dynamics (MD–SRD) for the fluid environment by means of stochastic collisions, in order to incorporate hydrodynamics effects leading to non-newtonian responses. The shear viscosity is computed using non-equilibrium simulations by imposing explicit velocity gradients. The impact of the colloid–colloid interactions is examined by modelling the inter-colloid pair potential with a repulsive power law, that allows interpolating different degrees of colloidal softness. The general rheological behaviour of our suspensions can be described with a Krieger–Dougherty like equation, which must be corrected to account for the variations in the maximum packing fraction and non-equilibrium effects arising from the flux of momentum imposed to the suspension, which appear when varying the softness of the inter-colloidal interactions. We further show evidence for non-newtonian behaviour at high Péclet numbers, characterised both by shear thinning and shear thickening, and thus demonstrate these effects can be successfully captured using MD–SRD methods.