The mechanical properties of clay minerals are largely dependent upon the chemical compositions and the mesoscale fabrics of the constituent particles. This paper describes results of a series of mesoscale molecular dynamics simulations of the hydrostatic compression and shear strain behavior for initially randomly oriented assemblies of 103 illite primary particles. The particles are simulated as rigid-body ellipsoids that interact through the single-site, Gay–Berne potential function. This corresponds to a coarse-grained model based on prior atomistic scale computation of the potential of mean force for water-mediated interactions between pairs of particles through free energy perturbation method. We investigate the mesoscale fabrics of the NPT-equilibrated assemblies for confining pressures ranging from 1.0 to 125 atm, including path dependence associated with unloading and reloading. We analyze and quantify the geometric arrangement including particle orientation, specific surface area, properties of particle stacks/aggregates, and interstack pair correlation functions. The compression of each particle assembly is associated with large irrecoverable changes in void ratio, while unloading and reloading involves much smaller, largely recoverable volumetric strains. The results are qualitatively similar to macroscopic compression behavior reported in laboratory tests. We simulate the uniaxial and shear behavior at each of the equilibrated pressure states through a series of strain-controlled steps, allowing full relaxation of the virial stresses computed at each step. The simulations investigate directional and path dependence of the shear behavior for strain deviations up to 0.2%. The results show the onset on nonlinear stiffness properties at strain levels \(\sim\)0.01% and hysteretic behavior upon unloading and reloading. Small-strain stiffness properties of the particle assemblies are qualitatively in good agreement with quasi-static, elastic stiffness properties reported for illitic clays.