  We extend current models of the halo occupation distribution (HOD) to include
a flexible, empirical framework for the forward modeling of the intrinsic
alignment (IA) of galaxies. A primary goal of this work is to produce mock
galaxy catalogs for the purpose of validating existing models and methods for
the mitigation of IA in weak lensing measurements. This technique can also be
used to produce new, simulation-based predictions for IA and galaxy clustering.
Our model is probabilistically formulated, and rests upon the assumption that
the orientations of galaxies exhibit a correlation with their host dark matter
(sub)halo orientation or with their position within the halo. We examine the
necessary components and phenomenology of such a model by considering the
alignments between (sub)halos in a cosmological dark matter only simulation. We
then validate this model for a realistic galaxy population in a set of
simulations in the IllustrisTNG suite. We create an HOD mock with TNG-like
correlations using our method, constraining the associated IA model parameters,
with the $\chi^2_{\rm dof}$ between our model's correlations and those of
Illustris matching as closely as 1.4 and 1.1 for orientation--position and
orientation--orientation correlation functions, respectively. By modeling the
misalignment between galaxies and their host halo, we show that the
3-dimensional two-point position and orientation correlation functions of
simulated (sub)halos and galaxies can be accurately reproduced from
quasi-linear scales down to $0.1~h^{-1}{\rm Mpc}$. We also find evidence for
environmental influence on IA within a halo. Our publicly-available software
provides a key component enabling efficient determination of Bayesian
posteriors on IA model parameters using observational measurements of
galaxy-orientation correlation functions in the highly nonlinear regime.
