© 2020 Author(s). Accurate calculation of the ion-ion recombination rate coefficient has been of long-standing interest as it controls the ion concentration in gas phase systems and in aerosols. We describe the development of a hybrid continuum-molecular dynamics (MD) approach to determine the ion-ion recombination rate coefficient. This approach is based on the limiting sphere method classically used for transition regime collision phenomena in aerosols. When ions are sufficiently far from one another, the ion-ion relative motion is described by diffusion equations, while within a critical distance, MD simulations are used to model ion-ion motion. MD simulations are parameterized using the Assisted Model Building with Energy Refinement force-field as well as by considering partial charges on atoms. Ion-neutral gas collisions are modeled in two mutually exclusive cubic domains composed of 103 gas atoms each, which remain centered on the recombining ions throughout calculations. Example calculations are reported for NH4+ recombination with NO2- in He, across a pressure range from 10 kPa to 10 000 kPa. Excellent agreement is found in comparison with calculations to literature values for the 100 kPa recombination rate coefficient (1.0 × 10-12 m3 s-1) in He. We also recover the experimentally observed increase in the recombination rate coefficient with pressure at sub-atmospheric pressures, and the observed decrease in the recombination rate coefficient in the high pressure continuum limit. We additionally find that non-dimensionalized forms of rate coefficients are consistent with recently developed equations for the dimensionless charged particle-ion collision rate coefficient based on Langevin dynamics simulations.