The addition of molecular dopants into organic semiconductors (OSCs) is a ubiquitous augmentation strategy to enhance the electrical conductivity of OSCs. Although the importance of optimizing OSC–dopant interactions is well-recognized, chemically generalizable structure–function relationships are difficult to extract due to the sensitivity and dependence of doping efficiency on chemistry, processing conditions, and morphology. Computational modeling for an integrated OSC–dopant design is an attractive approach to systematically isolate fundamental relationships, but requires the challenging simultaneous treatment of molecular reactivity and morphology evolution. We present the first computational study to couple molecular reactivity with morphology evolution in a molecularly doped OSC. Reactive Monte Carlo is employed to examine the evolution of OSC–dopant morphologies and doping efficiency with respect to dielectric, the thermodynamic driving for the doping reaction, and dopant aggregation. We observe that for well-mixed systems with experimentally relevant dielectric constants, doping efficiency is near unity with a very weak dependence on the ionization potential and electron affinity of OSC and dopant, respectively. At experimental dielectric constants, reaction-induced aggregation is observed, corresponding to the well-known insolubility of solution-doped materials. Simulations are qualitatively consistent with a number of experimental studies showing a decrease of doping efficiency with increasing dopant concentration. Finally, we observe that the aggregation of dopants lowers doping efficiency and thus presents a rational design strategy for maximizing doping efficiency in molecularly doped OSCs. This work represents an important first step toward the systematic integration of molecular reactivity and morphology evolution into the characterization of multi-scale structure–function relationships in molecularly doped OSCs.