Brownian dynamics (BD) simulation is used to calculate the viscoelasticity of model near-hard-sphere colloidal liquids using a continuous potential r–n interaction between the model colloidal particles. The exponent n was varied between 6 and 72. The real and reciprocal components of the complex shear viscosity, η′ and η″, were computed via time-correlation functions under no-shear conditions using a Green–Kubo formula. Also, oscillatory shear non-equilibrium BD was employed at finite strain amplitude in the linear-response regime. We find that the normalised stress autocorrelation function can be approximated very well by a two-parameter stretched exponential over the complete volume-fraction range. The parameters used to specify the stretched exponential and also the viscosities and long-time self-diffusion coefficients are quite sensitive to the value of n at a chosen volume fraction. The Newtonian viscosity decreases and the long-time self-diffusion coefficient increases with the softness of the interaction, in agreement with experiment. The value n= 36 gives best agreement with the experimental data, and therefore appears to be a good ‘effective’ interaction which we suggest includes the time-averaged effects of the many-body hydrodynamics to some extent. The state dependence of the derived spectrum of relaxation times is determined. As for experimental systems, the complex, viscosity scales with the dimensionless (‘longest’) relaxation time, Doτ1/a2, where a is the radius of the particle and Do is the self-diffusion coefficient in the zero-density limit. Also in the intermediate frequency regime 20 < a2ω/Do < 200 we find that both the real and imaginary parts of the complex shear viscosity decay as ca. ω–1/2 in agreement with experiment.