In this paper, an infinitesimal-strain based FFT formulation is extended to account for deformation twinning in hexagonal close-packed (HCP) materials. A model called the Complex Voxel (CV) model is developed that includes twinning as pseudo-slip and accounts for interaction between parent grains and corresponding twin variants by assuming a Taylor-type approximation at the voxel level. The macroscopic deformation behavior of Magnesium alloy AZ31, a representative HCP material, is simulated in three different loading directions. Detailed analysis of twinning dominated deformation reveals that: (a) Global Schmid factors indicate twin variant selection and the rate of twin growth; (b) Basal slip is the dominant active slip system during the twin nucleation and twin growth stages; (c) A comparison of numerical results of the present study with the published experimental and numerical studies indicate that the evolution of local stress states in parent grains and corresponding twin variants depend on their orientation with respect to the loading direction and/or the neighbouring grains. Further, the need of immediate experimental evidence to assist assumptions made during modelling deformation twinning is discussed. © 2018 Elsevier Ltd.