Liquid alkanes in the molecular weight range of C20–C40 are the main constituents of lubricant basestocks, and their rheological properties are therefore of great concern in industrial lubricant applications. Using massively parallel supercomputers and an efficient parallel algorithm, we have carried out systematic studies of the rheological properties of a variety of model liquid alkanes ranging from linear to singly branched and multiply branched alkanes. We aim to elucidate the relationship between the molecular architecture and the viscous behavior. Nonequilibrium molecular dynamics simulations have been carried out for n-decane (C10H22), n-hexadecane (C16H34), n-tetracosane (C24H50), 10-n-hexylnonadecane (C25H52), and squalane (2, 6, 10, 15, 19, 23-hexamethyltetracosane, C30H62). At a high strain rate, the viscosity shows a power-law shear thinning behavior over several orders of magnitude in strain rate, with exponents ranging from −0.33 to −0.59. This power-law shear thinning is shown to be closely related to the ordering of the molecules. The molecular architecture is shown to have a significant influence on the power-law exponent. At a low strain rate, the viscosity behavior changes to a Newtonian plateau, whose accurate determination has been elusive in previous studies. The molecular order in this regime is essentially that of the equilibrium system, a signature of the linear response. The Newtonian plateau is verified by independent equilibrium molecular dynamics simulations using the Green–Kubo method. The reliable determination of the Newtonian viscosity from non-equilibrium molecular simulation permits us to calculate the viscosity index for squalane. The viscosity index is a widely used property to characterize the lubricant's temperature performance, and our studies represent the first approach towards its determination by molecular simulation.