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 10 3 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 – 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 inter-stack 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 volu-metric strains. The results are qualitatively similar to macroscopic compression 1 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 non-linear stiffness properties at strain levels ∼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.