  We construct the S-wave part of the electromagnetic vector annihilation
current to $O(\alpha_s v^2)$, where $v$ is the non-relativistic quark velocity,
for heavy quarks whose dynamics are described by the NRQCD action on the
lattice. The NRQCD vector current for $Q\bar{Q}$ annihilation is expressed as a
linear combination of lattice operators with quantum numbers L=0, $J^P=1^-$,
and the coefficients are determined by matching to the corresponding continuum
current in QCD to$O(v^2)$ at one-loop. The annihilation channel gives a complex
amplitude with Coulomb-exchange and infrared singularities, making a careful
choice for the contours of integration and infrared subtraction functions in
the numerical integration necessary. An automated vertex generation program
written in Python is employed, allowing us to use a realistic NRQCD action and
an improved gluon lattice action; a change in the definition of either action
is easily accommodated in this procedure. The final result is applicable to
simulations of electromagnetic decays of heavy quarkonia, notably the
$\Upsilon$ meson.
