  To understand massive star formation requires study of its initial
conditions. Two massive starless core candidates, C1-N & C1-S, have been
detected in IRDC G028.37+00.07 in $\rm N_2D^+$(3-2) with $ALMA$. From their
line widths, either the cores are subvirial and are thus young structures on
the verge of near free-fall collapse, or they are threaded by $\sim1$ mG
$B$-fields that help support them in near virial equilibrium and potentially
have older ages. We modeled the deuteration rate of $\rm N_2H^+$ to constrain
collapse rates of the cores. First, to measure their current deuterium
fraction, $D_{\rm
  frac}^{\rm N_2H^+}$ $\equiv [\rm N_2D^+]/[N_2H^+]$, we observed multiple
transitions of $\rm N_2H^+$ and $\rm N_2D^+$ with $CARMA$, $SMA$, $JCMT$,
$NRO~45m$ and $IRAM~30m$, to complement the $ALMA$ data. For both cores we
derived $D_{\rm
  frac}^{\rm N_2H^+}\sim0.3$, several orders of magnitude above the cosmic
[D]/[H] ratio. We then carried out chemodynamical modeling, exploring how
collapse rate relative to free-fall, $\alpha_{\rm ff}$, affects the level of
$D_{\rm frac}^{\rm N_2H^+}$ that is achieved from a given initial condition. To
reach the observed $D_{\rm frac}^{\rm
  N_2H^+}$, most models require slow collapse with $\alpha_{\rm
  ff}\sim0.1$, i.e., $\sim1/10$th of free-fall. This makes it more likely that
the cores have been able to reach a near virial equilibrium state and we
predict that strong $B$-fields will eventually be detected. The methods
developed here will be useful for measurement of the pre-stellar core mass
function.
