ABSTRACT We present a semi-analytic model for predicting kilonova light curves from the mergers of neutron stars with black holes (NSBH). The model is integrated into the mosfit platform, and can generate light curves from input binary properties and nuclear equation-of-state considerations, or incorporate measurements from gravitational wave (GW) detectors to perform multimessenger parameter estimation. The rapid framework enables the generation of NSBH kilonova distributions from binary populations, light curve predictions from GW data, and statistically meaningful comparisons with an equivalent binary neutron star (BNS) model in mosfit. We investigate a sample of kilonova candidates associated with cosmological short gamma-ray bursts, and demonstrate that they are broadly consistent with being driven by NSBH systems, though most have limited data. We also perform fits to the very well sampled GW170817, and show that the inability of an NSBH merger to produce lanthanide-poor ejecta results in a significant underestimate of the early (≲2 d) optical emission. Our model indicates that NSBH-driven kilonovae may peak up to a week after merger at optical wavelengths for some observer angles. This demonstrates the need for early coverage of emergent kilonovae in cases where the GW signal is either ambiguous or absent; they likely cannot be distinguished from BNS mergers by the light curves alone from ∼2 d after the merger. We also discuss the detectability of our model kilonovae with the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST).