  We propose a new method for calculating the Chern-Simons orbital
magnetoelectric coupling, conventionally parametrized in terms of a phase angle
$\theta$. According to previous theories, $\theta$ can be expressed as a 3D
Brillouin-zone integral of the Chern-Simons 3-form defined in terms of the
occupied Bloch functions. Such an expression is valid only if a smooth and
periodic gauge has been chosen in the entire Brillouin zone, and even then,
convergence with respect to the $\mathbf{k}$-space mesh density can be
difficult to obtain. In order to solve this problem, we propose to relax the
periodicity condition in one direction (say, the $k_z$ direction) so that a
gauge discontinuity is introduced on a 2D $\mathbf{k}$ plane normal to $k_z$.
The total $\theta$ response then has contributions from both the integral of
the Chern-Simons 3-form over the 3D bulk BZ and the gauge discontinuity
expressed as a 2D integral over the $\mathbf{k}$ plane. Sometimes the boundary
plane may be further divided into subregions by 1D "vortex loops" which make a
third kind of contribution to the total $\theta$, expressed as a combination of
Berry phases around the vortex loops. The total $\theta$ thus consists of three
terms which can be expressed as integrals over 3D, 2D and 1D manifolds. When
time-reversal symmetry is present and the gauge in the bulk BZ is chosen to
respect this symmetry, both the 3D and 2D integrals vanish; the entire
contribution then comes from the vortex-loop integral, which is either 0 or
$\pi$ corresponding to the $\mathbb{Z}_2$ classification of 3D time-reversal
invariant insulators. We demonstrate our method by applying it to the
Fu-Kane-Mele model with an applied staggered Zeeman field.
