diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9facdf7ba5..b72bb57eb9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -345,7 +345,6 @@ add_library(libopenmc SHARED
   src/photon_physics.F90
   src/physics_common.F90
   src/physics.F90
-  src/physics_mg.F90
   src/plot.F90
   src/plot_header.F90
   src/progress_header.F90
@@ -423,6 +422,8 @@ add_library(libopenmc SHARED
   src/nuclide.cpp
   src/output.cpp
   src/particle.cpp
+  src/physics_common.cpp
+  src/physics_mg.cpp
   src/plot.cpp
   src/position.cpp
   src/pugixml/pugixml_c.cpp
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index c48b212bd7..57d817d23d 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -149,6 +149,8 @@ extern "C" {
   extern int32_t n_tallies;
   extern int32_t n_universes;
   extern bool openmc_simulation_initialized;
+  extern double global_tally_absorption;
+#pragma omp threadprivate(global_tally_absorption)
 
   // Variables that are shared by necessity (can be removed from public header
   // later)
diff --git a/include/openmc/material.h b/include/openmc/material.h
index 1af6f85976..909215ee37 100644
--- a/include/openmc/material.h
+++ b/include/openmc/material.h
@@ -26,6 +26,7 @@ class Material
 public:
   int32_t id; //!< Unique ID
   double volume_ {-1.0}; //!< Volume in [cm^3]
+  bool fissionable {false}; //!< Does this material contain fissionable nuclides
 
   //! \brief Default temperature for cells containing this material.
   //!
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index a3c9c07cf0..8bfa118995 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -22,7 +22,7 @@ namespace openmc {
 //! @return The requested percentile
 //==============================================================================
 
-extern "C" double normal_percentile_c(double p);
+extern "C" double normal_percentile(double p);
 
 //==============================================================================
 //! Calculate the percentile of the Student's t distribution with a specified
@@ -58,7 +58,7 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
 //!   evaluated at x
 //==============================================================================
 
-extern "C" double evaluate_legendre_c(int n, const double data[], double x);
+extern "C" double evaluate_legendre(int n, const double data[], double x);
 
 //==============================================================================
 //! Calculate the n-th order real spherical harmonics for a given angle (in
@@ -144,7 +144,7 @@ Direction rotate_angle(Direction u, double mu, double* phi);
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double maxwell_spectrum_c(double T);
+extern "C" double maxwell_spectrum(double T);
 
 //==============================================================================
 //! Samples an energy from a Watt energy-dependent fission distribution.
@@ -159,7 +159,7 @@ extern "C" double maxwell_spectrum_c(double T);
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double watt_spectrum_c(double a, double b);
+extern "C" double watt_spectrum(double a, double b);
 
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.
diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h
index 02a43fece4..e0addf5bc0 100644
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@@ -18,6 +18,9 @@ extern std::vector<Mgxs> nuclides_MG;
 extern std::vector<Mgxs> macro_xs;
 extern "C" int num_energy_groups;
 
+//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
+// bring over energy_bin_avg, energy_bins, etc, as vectors)
+
 //==============================================================================
 // Mgxs data loading interface methods
 //==============================================================================
@@ -44,13 +47,6 @@ extern "C" void
 calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
      double& total_xs, double& abs_xs, double& nu_fiss_xs);
 
-extern "C" void
-sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
-     double uvw[3]);
-
-extern "C" void
-sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
-
 extern "C" double
 get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
 
diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h
new file mode 100644
index 0000000000..33fb8ed614
--- /dev/null
+++ b/include/openmc/physics_common.h
@@ -0,0 +1,16 @@
+//! \file physics_common.h
+//! A collection of physics methods common to MG, CE, photon, etc.
+
+#ifndef OPENMC_PHYSICS_COMMON_H
+#define OPENMC_PHYSICS_COMMON_H
+
+#include "openmc/particle.h"
+
+namespace openmc {
+
+//! \brief Performs the russian roulette operation for a particle
+extern "C" void
+russian_roulette(Particle* p);
+
+} // namespace openmc
+#endif // OPENMC_PHYSICS_COMMON_H
diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h
new file mode 100644
index 0000000000..14b7e6cbd5
--- /dev/null
+++ b/include/openmc/physics_mg.h
@@ -0,0 +1,60 @@
+//! \file physics_mg.h
+//! Methods needed to perform the collision physics for multi-group mode
+
+#ifndef OPENMC_PHYSICS_MG_H
+#define OPENMC_PHYSICS_MG_H
+
+#include "openmc/capi.h"
+#include "openmc/particle.h"
+#include "openmc/nuclide.h"
+
+namespace openmc {
+
+//TODO: Remove energy_bin_avg and material_xs parameters when they reside on
+// the C-side this should happen after materials, physics, input, and tallies
+// are brought over
+
+//! \brief samples particle behavior after a collision event.
+//! \param p Particle to operate on
+//! \param energy_bin_avg Average energy within each energy bin
+//! \param material_xs The cross section cache for the current material
+extern "C" void
+collision_mg(Particle* p, const double* energy_bin_avg,
+     const MaterialMacroXS* material_xs);
+
+//! \brief samples a reaction type.
+//!
+//! Note that there is special logic when suvival biasing is turned on since
+//! fission and disappearance are treated implicitly.
+//! \param p Particle to operate on
+//! \param energy_bin_avg Average energy within each energy bin
+//! \param material_xs The cross section cache for the current material
+void
+sample_reaction(Particle* p, const double* energy_bin_avg,
+     const MaterialMacroXS* material_xs);
+
+//! \brief Samples the scattering event
+//! \param p Particle to operate on
+//! \param energy_bin_avg Average energy within each energy bin
+void
+scatter(Particle* p, const double* energy_bin_avg);
+
+//! \brief Determines the average total, prompt and delayed neutrons produced
+//! from fission and creates the appropriate bank sites.
+//! \param p Particle to operate on
+//! \param bank_array The particle bank to populate
+//! \param size_bank Number of particles currently in the bank
+//! \param bank_array_size Allocated size of the bank
+//! \param material_xs The cross section cache for the current material
+void
+create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
+     int64_t bank_array_size, const MaterialMacroXS* material_xs);
+
+//! \brief Handles an absorption event
+//! \param p Particle to operate on
+//! \param material_xs The cross section cache for the current material
+void
+absorption(Particle* p, const MaterialMacroXS* material_xs);
+
+} // namespace openmc
+#endif // OPENMC_PHYSICS_MG_H
diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 390ee466ff..a8ea849b3e 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -12,8 +12,8 @@ _dll.t_percentile_c.argtypes = [c_double, c_int]
 _dll.calc_pn_c.restype = None
 _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
-_dll.evaluate_legendre_c.restype = c_double
-_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
+_dll.evaluate_legendre.restype = c_double
+_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double]
 
 _dll.calc_rn_c.restype = None
 _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
@@ -27,11 +27,11 @@ _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
                                 POINTER(c_double)]
-_dll.maxwell_spectrum_c.restype = c_double
-_dll.maxwell_spectrum_c.argtypes = [c_double]
+_dll.maxwell_spectrum.restype = c_double
+_dll.maxwell_spectrum.argtypes = [c_double]
 
-_dll.watt_spectrum_c.restype = c_double
-_dll.watt_spectrum_c.argtypes = [c_double, c_double]
+_dll.watt_spectrum.restype = c_double
+_dll.watt_spectrum.argtypes = [c_double, c_double]
 
 _dll.broaden_wmp_polynomials_c.restype = None
 _dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
@@ -100,9 +100,8 @@ def evaluate_legendre(data, x):
     """
 
     data_arr = np.array(data, dtype=np.float64)
-    return _dll.evaluate_legendre_c(len(data),
-                                    data_arr.ctypes.data_as(POINTER(c_double)),
-                                    x)
+    return _dll.evaluate_legendre(len(data),
+                                  data_arr.ctypes.data_as(POINTER(c_double)), x)
 
 
 def calc_rn(n, uvw):
@@ -182,7 +181,7 @@ def calc_zn_rad(n, rho):
     zn_rad = np.zeros(num_bins, dtype=np.float64)
     _dll.calc_zn_rad_c(n, rho, zn_rad)
     return zn_rad
-    
+
 
 def rotate_angle(uvw0, mu, phi=None):
     """ Rotates direction cosines through a polar angle whose cosine is
@@ -231,7 +230,7 @@ def maxwell_spectrum(T):
 
     """
 
-    return _dll.maxwell_spectrum_c(T)
+    return _dll.maxwell_spectrum(T)
 
 
 def watt_spectrum(a, b):
@@ -251,7 +250,7 @@ def watt_spectrum(a, b):
 
     """
 
-    return _dll.watt_spectrum_c(a, b)
+    return _dll.watt_spectrum(a, b)
 
 
 def broaden_wmp_polynomials(E, dopp, n):
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 01b09e65db..e546b11ddd 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -90,7 +90,7 @@ Maxwell::Maxwell(pugi::xml_node node)
 
 double Maxwell::sample() const
 {
-  return maxwell_spectrum_c(theta_);
+  return maxwell_spectrum(theta_);
 }
 
 //==============================================================================
@@ -110,7 +110,7 @@ Watt::Watt(pugi::xml_node node)
 
 double Watt::sample() const
 {
-  return watt_spectrum_c(a_, b_);
+  return watt_spectrum(a_, b_);
 }
 
 //==============================================================================
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 863ffd8b08..e311e322d2 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -279,7 +279,7 @@ double MaxwellEnergy::sample(double E) const
 
   while (true) {
     // Sample maxwell fission spectrum
-    double E_out = maxwell_spectrum_c(theta);
+    double E_out = maxwell_spectrum(theta);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
@@ -343,7 +343,7 @@ double WattEnergy::sample(double E) const
 
   while (true) {
     // Sample energy-dependent Watt fission spectrum
-    double E_out = watt_spectrum_c(a, b);
+    double E_out = watt_spectrum(a, b);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 420ae5bbb8..989f7b357f 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2986,7 +2986,7 @@ contains
 
         ! Check if material is fissionable
         if (nuclides(materials(i) % nuclide(j)) % fissionable) then
-          materials(i) % fissionable = .true.
+          call materials(i) % set_fissionable(.true.)
         end if
       end do
 
diff --git a/src/material.cpp b/src/material.cpp
index 1a0c2d2078..ecfc7f55e2 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -122,6 +122,13 @@ extern "C" {
     material_map[id] = index - 1;
   }
 
+  bool material_fissionable(Material* mat) {return mat->fissionable;}
+
+  void material_set_fissionable(Material* mat, bool fissionable)
+  {
+    mat->fissionable = fissionable;
+  }
+
   void extend_materials_c(int32_t n)
   {
     materials.reserve(materials.size() + n);
diff --git a/src/material_header.F90 b/src/material_header.F90
index 71e6ec5f1a..b4d0617255 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -51,6 +51,20 @@ module material_header
       integer(C_INT32_T), intent(in), value :: index
     end subroutine material_set_id_c
 
+    function material_fissionable_c(mat_ptr) &
+         bind(C, name='material_fissionable') result(fissionable)
+      import C_PTR, C_BOOL
+      type(C_PTR), intent(in), value :: mat_ptr
+      logical(C_BOOL)                :: fissionable
+    end function material_fissionable_c
+
+    subroutine material_set_fissionable_c(mat_ptr, fissionable) &
+         bind(C, name='material_set_fissionable')
+      import C_PTR, C_BOOL
+      type(C_PTR),        intent(in), value :: mat_ptr
+      logical(C_BOOL),    intent(in), value :: fissionable
+    end subroutine material_set_fissionable_c
+
     subroutine extend_materials_c(n) bind(C)
       import C_INT32_T
       integer(C_INT32_T), intent(in), value :: n
@@ -95,7 +109,6 @@ module material_header
     character(20), allocatable :: sab_names(:) ! name of S(a,b) table
 
     ! Does this material contain fissionable nuclides? Is it depletable?
-    logical :: fissionable = .false.
     logical :: depletable = .false.
 
     ! enforce isotropic scattering in lab for specific nuclides
@@ -105,6 +118,8 @@ module material_header
   contains
     procedure :: id => material_id
     procedure :: set_id => material_set_id
+    procedure :: fissionable => material_fissionable
+    procedure :: set_fissionable => material_set_fissionable
     procedure :: set_density => material_set_density
     procedure :: init_nuclide_index => material_init_nuclide_index
     procedure :: assign_sab_tables => material_assign_sab_tables
@@ -139,6 +154,18 @@ contains
     call material_set_id_c(this % ptr, id, index)
   end subroutine material_set_id
 
+  function material_fissionable(this) result(fissionable)
+    class(Material), intent(in) :: this
+    logical(C_BOOL)             :: fissionable
+    fissionable = material_fissionable_c(this % ptr)
+  end function material_fissionable
+
+  subroutine material_set_fissionable(this, fissionable)
+    class(Material),intent(in) :: this
+    logical,        intent(in) :: fissionable
+    call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL))
+  end subroutine material_set_fissionable
+
   function material_set_density(this, density) result(err)
     class(Material), intent(inout) :: this
     real(8), intent(in) :: density
@@ -685,7 +712,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= size(materials)) then
-      fissionable = materials(index) % fissionable
+      fissionable = materials(index) % fissionable()
       err = 0
     else
       err = E_OUT_OF_BOUNDS
diff --git a/src/math.F90 b/src/math.F90
index 9964922c17..143e21e030 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -12,10 +12,7 @@ module math
   public :: calc_rn
   public :: calc_zn
   public :: calc_zn_rad
-  public :: evaluate_legendre
   public :: rotate_angle
-  public :: maxwell_spectrum
-  public :: watt_spectrum
   public :: faddeeva
   public :: w_derivative
   public :: broaden_wmp_polynomials
@@ -42,16 +39,6 @@ module math
       real(C_DOUBLE), intent(out) :: pnx(n + 1)
     end subroutine calc_pn
 
-    pure function evaluate_legendre_c_intfc(n, data, x) &
-         bind(C, name='evaluate_legendre_c') result(val)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: n
-      real(C_DOUBLE), intent(in) :: data(n)
-      real(C_DOUBLE), value, intent(in) :: x
-      real(C_DOUBLE) :: val
-    end function evaluate_legendre_c_intfc
-
     pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
       use ISO_C_BINDING
       implicit none
@@ -85,23 +72,6 @@ module math
       real(C_DOUBLE), optional, intent(in) :: phi
     end subroutine rotate_angle_c_intfc
 
-    function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') &
-         result(E_out)
-      use ISO_C_BINDING
-      implicit none
-      real(C_DOUBLE), value, intent(in) :: T
-      real(C_DOUBLE) :: E_out
-    end function maxwell_spectrum
-
-    function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') &
-         result(E_out)
-      use ISO_C_BINDING
-      implicit none
-      real(C_DOUBLE), value, intent(in) :: a
-      real(C_DOUBLE), value, intent(in) :: b
-      real(C_DOUBLE) :: E_out
-    end function watt_spectrum
-
     subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
            bind(C, name='broaden_wmp_polynomials_c')
       use ISO_C_BINDING
@@ -157,20 +127,6 @@ module math
 
 contains
 
-!===============================================================================
-! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
-! and the value of x
-!===============================================================================
-
-  pure function evaluate_legendre(data, x) result(val) bind(C)
-    real(C_DOUBLE), intent(in) :: data(:)
-    real(C_DOUBLE), intent(in) :: x
-    real(C_DOUBLE)             :: val
-
-    val = evaluate_legendre_c_intfc(size(data) - 1, data, x)
-
-  end function evaluate_legendre
-
 !===============================================================================
 ! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
 ! mu and through an azimuthal angle sampled uniformly. Note that this is done
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 5c1a2f4dce..8d30ea84d7 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -6,7 +6,7 @@ namespace openmc {
 // Mathematical methods
 //==============================================================================
 
-double normal_percentile_c(double p) {
+double normal_percentile(double p) {
   constexpr double p_low = 0.02425;
   constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
                            -2.759285104469687e2, 1.383577518672690e2,
@@ -79,7 +79,7 @@ double t_percentile_c(double p, int df){
     // 16 (4), pp. 1123-1132 (1987).
     double n = df;
     double k = 1. / (n - 2.);
-    double z = normal_percentile_c(p);
+    double z = normal_percentile(p);
     double z2 = z * z;
     t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 +
          75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) *
@@ -103,7 +103,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
 }
 
 
-double evaluate_legendre_c(int n, const double data[], double x) {
+double evaluate_legendre(int n, const double data[], double x) {
   double pnx[n + 1];
   double val = 0.0;
   calc_pn_c(n, x, pnx);
@@ -658,7 +658,7 @@ Direction rotate_angle(Direction u, double mu, double* phi)
 }
 
 
-double maxwell_spectrum_c(double T) {
+double maxwell_spectrum(double T) {
   // Set the random numbers
   double r1 = prn();
   double r2 = prn();
@@ -674,8 +674,8 @@ double maxwell_spectrum_c(double T) {
 }
 
 
-double watt_spectrum_c(double a, double b) {
-  double w = maxwell_spectrum_c(a);
+double watt_spectrum(double a, double b) {
+  double w = maxwell_spectrum(a);
   double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);
 
   return E_out;
diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90
index 03a0e9402f..90178a015d 100644
--- a/src/mgxs_data.F90
+++ b/src/mgxs_data.F90
@@ -90,7 +90,8 @@ contains
         end if
       end do NUCLIDE_LOOP
 
-      mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide)
+      call mat % set_fissionable( &
+           logical(query_fissionable_c(mat % n_nuclides, mat % nuclide)))
 
     end do MATERIAL_LOOP
 
diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90
index 9dccc92301..77f629b80e 100644
--- a/src/mgxs_interface.F90
+++ b/src/mgxs_interface.F90
@@ -62,26 +62,6 @@ module mgxs_interface
       real(C_DOUBLE),     intent(inout) :: nu_fiss_xs
     end subroutine calculate_xs_c
 
-    subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: i_mat
-      integer(C_INT), value, intent(in) :: gin
-      integer(C_INT),     intent(inout) :: gout
-      real(C_DOUBLE),     intent(inout) :: mu
-      real(C_DOUBLE),     intent(inout) :: wgt
-      real(C_DOUBLE),     intent(inout) :: uvw(1:3)
-    end subroutine sample_scatter_c
-
-    subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: i_mat
-      integer(C_INT), value, intent(in) :: gin
-      integer(C_INT),     intent(inout) :: dg
-      integer(C_INT),     intent(inout) :: gout
-    end subroutine sample_fission_energy_c
-
     subroutine get_name_c(index, name_len, name) bind(C)
       use ISO_C_BINDING
       implicit none
@@ -156,7 +136,7 @@ module mgxs_interface
   real(8), allocatable :: energy_bins(:)
 
   ! Midpoint of the energy group structure
-  real(8), allocatable :: energy_bin_avg(:)
+  real(C_DOUBLE), allocatable :: energy_bin_avg(:)
 
   ! Energy group structure with increasing energy
   real(C_DOUBLE), allocatable, target :: rev_energy_bins(:)
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 0bb464da07..c592d28a28 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -97,36 +97,6 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
 
 //==============================================================================
 
-void
-sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
-                 double uvw[3])
-{
-  int gout_c = gout - 1;
-  macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
-
-  // adjust return value for fortran indexing
-  gout = gout_c + 1;
-
-  // Rotate the angle
-  rotate_angle_c(uvw, mu, nullptr);
-}
-
-//==============================================================================
-
-void
-sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
-{
-  int dg_c = 0;
-  int gout_c = 0;
-  macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
-
-  // adjust return values for fortran indexing
-  dg = dg_c + 1;
-  gout = gout_c + 1;
-}
-
-//==============================================================================
-
 double
 get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
 {
diff --git a/src/physics_common.F90 b/src/physics_common.F90
index 63bf62ce62..0d5b9fb324 100644
--- a/src/physics_common.F90
+++ b/src/physics_common.F90
@@ -1,33 +1,20 @@
 module physics_common
 
-  use constants
+  use, intrinsic :: ISO_C_BINDING
+
   use particle_header,        only: Particle
-  use random_lcg,             only: prn
-  use settings,               only: weight_cutoff, weight_survive
 
   implicit none
 
-contains
-
 !===============================================================================
-! RUSSIAN_ROULETTE
+! RUSSIAN_ROULETTE FROM C
 !===============================================================================
 
-  subroutine russian_roulette(p)
-
-    type(Particle), intent(inout) :: p
-
-    if (p % wgt < weight_cutoff) then
-      if (prn() < p % wgt / weight_survive) then
-        p % wgt = weight_survive
-        p % last_wgt = p % wgt
-      else
-        p % wgt = ZERO
-        p % last_wgt = ZERO
-        p % alive = .false.
-      end if
-    end if
-
-  end subroutine russian_roulette
+  interface
+    subroutine russian_roulette(p) bind(C)
+      import Particle
+      type(Particle), intent(inout) :: p
+    end subroutine russian_roulette
+  end interface
 
 end module physics_common
diff --git a/src/physics_common.cpp b/src/physics_common.cpp
new file mode 100644
index 0000000000..02c77ae7d9
--- /dev/null
+++ b/src/physics_common.cpp
@@ -0,0 +1,26 @@
+#include "openmc/physics_common.h"
+
+#include "openmc/settings.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+//==============================================================================
+// RUSSIAN_ROULETTE
+//==============================================================================
+
+void russian_roulette(Particle* p)
+{
+  if (p->wgt < settings::weight_cutoff) {
+    if (prn() < p->wgt / settings::weight_survive) {
+      p->wgt = settings::weight_survive;
+      p->last_wgt = p->wgt;
+    } else {
+      p->wgt = 0.;
+      p->last_wgt = 0.;
+      p->alive = false;
+    }
+  }
+}
+
+} //namespace openmc
diff --git a/src/physics_mg.F90 b/src/physics_mg.F90
deleted file mode 100644
index 1d60b371ee..0000000000
--- a/src/physics_mg.F90
+++ /dev/null
@@ -1,275 +0,0 @@
-module physics_mg
-  ! This module contains the multi-group specific physics routines so as to not
-  ! hinder performance of the CE versions with multiple if-thens.
-
-  use bank_header
-  use constants
-  use error,                  only: fatal_error, warning, write_message
-  use material_header,        only: Material, materials
-  use math,                   only: rotate_angle
-  use mgxs_interface
-  use message_passing
-  use nuclide_header,         only: material_xs
-  use particle_header
-  use physics_common
-  use random_lcg,             only: prn
-  use settings
-  use simulation_header
-  use string,                 only: to_str
-  use tally_header
-
-  implicit none
-
-contains
-
-!===============================================================================
-! COLLISION_MG samples a nuclide and reaction and then calls the appropriate
-! routine for that reaction
-!===============================================================================
-
-  subroutine collision_mg(p)
-
-    type(Particle), intent(inout) :: p
-
-    ! Add to collision counter for particle
-    p % n_collision = p % n_collision + 1
-
-    ! Sample nuclide/reaction for the material the particle is in
-    call sample_reaction(p)
-
-    ! Display information about collision
-    if (verbosity >= 10 .or. trace) then
-      call write_message("    " // "Energy Group = " // trim(to_str(p % g)))
-    end if
-
-  end subroutine collision_mg
-
-!===============================================================================
-! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
-! each nuclide within a material and then samples a reaction for that nuclide
-! and calls the appropriate routine to process the physics. Note that there is
-! special logic when suvival biasing is turned on since fission and
-! disappearance are treated implicitly.
-!===============================================================================
-
-  subroutine sample_reaction(p)
-
-    type(Particle), intent(inout) :: p
-
-    type(Material), pointer :: mat
-
-    mat => materials(p % material)
-
-    ! Create fission bank sites. Note that while a fission reaction is sampled,
-    ! it never actually "happens", i.e. the weight of the particle does not
-    ! change when sampling fission sites. The following block handles all
-    ! absorption (including fission)
-
-    if (mat % fissionable) then
-      if (run_mode == MODE_EIGENVALUE) then
-        call create_fission_sites(p, fission_bank, n_bank)
-      elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
-        call create_fission_sites(p, p % secondary_bank, p % n_secondary)
-      end if
-    end if
-
-    ! If survival biasing is being used, the following subroutine adjusts the
-    ! weight of the particle. Otherwise, it checks to see if absorption occurs
-
-    if (material_xs % absorption > ZERO) then
-      call absorption(p)
-    else
-      p % absorb_wgt = ZERO
-    end if
-    if (.not. p % alive) return
-
-    ! Sample a scattering reaction and determine the secondary energy of the
-    ! exiting neutron
-    call scatter(p)
-
-    ! Play russian roulette if survival biasing is turned on
-    if (survival_biasing) then
-      call russian_roulette(p)
-      if (.not. p % alive) return
-    end if
-
-  end subroutine sample_reaction
-
-!===============================================================================
-! ABSORPTION
-!===============================================================================
-
-  subroutine absorption(p)
-
-    type(Particle), intent(inout) :: p
-
-    if (survival_biasing) then
-      ! Determine weight absorbed in survival biasing
-      p % absorb_wgt = (p % wgt * &
-                        material_xs % absorption / material_xs % total)
-
-      ! Adjust weight of particle by probability of absorption
-      p % wgt = p % wgt - p % absorb_wgt
-      p % last_wgt = p % wgt
-
-      ! Score implicit absorption estimate of keff
-!$omp atomic
-      global_tallies(RESULT_VALUE, K_ABSORPTION) = &
-           global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * &
-           material_xs % nu_fission / material_xs % absorption
-    else
-      ! See if disappearance reaction happens
-      if (material_xs % absorption > prn() * material_xs % total) then
-        ! Score absorption estimate of keff
-!$omp atomic
-        global_tallies(RESULT_VALUE, K_ABSORPTION) = &
-             global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * &
-             material_xs % nu_fission / material_xs % absorption
-
-        p % alive = .false.
-        p % event = EVENT_ABSORB
-      end if
-    end if
-
-  end subroutine absorption
-
-!===============================================================================
-! SCATTER
-!===============================================================================
-
-  subroutine scatter(p)
-
-    type(Particle), intent(inout)  :: p
-
-    call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
-         p % wgt, p % coord(1) % uvw)
-
-    ! Update energy value for downstream compatability (in tallying)
-    p % E = energy_bin_avg(p % g)
-
-    ! Set event component
-    p % event = EVENT_SCATTER
-
-  end subroutine scatter
-
-!===============================================================================
-! CREATE_FISSION_SITES determines the average total, prompt, and delayed
-! neutrons produced from fission and creates appropriate bank sites.
-!===============================================================================
-
-  subroutine create_fission_sites(p, bank_array, size_bank)
-    type(Particle), intent(inout) :: p
-    type(Bank),     intent(inout) :: bank_array(:)
-    integer(8),     intent(inout) :: size_bank
-
-    integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
-    integer :: i                        ! loop index
-    integer :: dg                       ! delayed group
-    integer :: gout                     ! group out
-    integer :: nu                       ! actual number of neutrons produced
-    real(8) :: nu_t                     ! total nu
-    real(8) :: mu                       ! fission neutron angular cosine
-    real(8) :: phi                      ! fission neutron azimuthal angle
-    real(8) :: weight                   ! weight adjustment for ufs method
-
-    interface
-      function ufs_get_weight(p) result(weight) bind(C)
-        import Particle, C_DOUBLE
-        type(Particle), intent(in) :: p
-        real(C_DOUBLE) :: WEIGHT
-      end function
-    end interface
-
-    ! TODO: Heat generation from fission
-
-    ! If uniform fission source weighting is turned on, we increase of decrease
-    ! the expected number of fission sites produced
-
-    if (ufs) then
-      weight = ufs_get_weight(p)
-    else
-      weight = ONE
-    end if
-
-    ! Determine expected number of neutrons produced
-    nu_t = p % wgt / keff * weight * &
-         material_xs % nu_fission / material_xs % total
-
-    ! Sample number of neutrons produced
-    if (prn() > nu_t - int(nu_t)) then
-      nu = int(nu_t)
-    else
-      nu = int(nu_t) + 1
-    end if
-
-    ! Check for bank size getting hit. For fixed source calculations, this is a
-    ! fatal error. For eigenvalue calculations, it just means that k-effective
-    ! was too high for a single batch.
-    if (size_bank + nu > size(bank_array)) then
-      if (run_mode == MODE_FIXEDSOURCE) then
-        call fatal_error("Secondary particle bank size limit reached. If you &
-             &are running a subcritical multiplication problem, k-effective &
-             &may be too close to one.")
-      else
-        if (master) call warning("Maximum number of sites in fission bank &
-             &reached. This can result in irreproducible results using different &
-             &numbers of processes/threads.")
-      end if
-    end if
-
-    ! Bank source neutrons
-    if (nu == 0 .or. size_bank == size(bank_array)) return
-
-    ! Initialize counter of delayed neutrons encountered for each delayed group
-    ! to zero.
-    nu_d(:) = 0
-
-    p % fission = .true. ! Fission neutrons will be banked
-    do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
-      ! Bank source neutrons by copying particle data
-      bank_array(i) % xyz = p % coord(1) % xyz
-
-      ! Set particle as neutron
-      bank_array(i) % particle = NEUTRON
-
-      ! Set weight of fission bank site
-      bank_array(i) % wgt = ONE/weight
-
-      ! Sample cosine of angle -- fission neutrons are treated as being emitted
-      ! isotropically.
-      mu = TWO * prn() - ONE
-
-      ! Sample azimuthal angle uniformly in [0,2*pi)
-      phi = TWO * PI * prn()
-      bank_array(i) % uvw(1) = mu
-      bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
-      bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
-
-      ! Sample secondary energy distribution for fission reaction and set energy
-      ! in fission bank
-      call sample_fission_energy_c(p % material, p % g, dg, gout)
-
-      bank_array(i) % E = real(gout, 8)
-      bank_array(i) % delayed_group = dg
-
-      ! Set delayed group on particle too
-      p % delayed_group = dg
-
-      ! Increment the number of neutrons born delayed
-      if (p % delayed_group > 0) then
-        nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1
-      end if
-
-    end do
-
-    ! increment number of bank sites
-    size_bank = min(size_bank + nu, int(size(bank_array),8))
-
-    ! Store total weight banked for analog fission tallies
-    p % n_bank   = nu
-    p % wgt_bank = nu/weight
-    p % n_delayed_bank(:) = nu_d(:)
-
-  end subroutine create_fission_sites
-
-end module physics_mg
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
new file mode 100644
index 0000000000..38eb9a5b95
--- /dev/null
+++ b/src/physics_mg.cpp
@@ -0,0 +1,233 @@
+#include "openmc/physics_mg.h"
+
+#include <stdexcept>
+#include <sstream>
+
+#include "xtensor/xarray.hpp"
+
+#include "openmc/constants.h"
+#include "openmc/eigenvalue.h"
+#include "openmc/error.h"
+#include "openmc/material.h"
+#include "openmc/math_functions.h"
+#include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
+#include "openmc/physics_common.h"
+#include "openmc/random_lcg.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+
+namespace openmc {
+
+void
+collision_mg(Particle* p, const double* energy_bin_avg,
+             const MaterialMacroXS* material_xs)
+{
+  // Add to the collision counter for the particle
+  p->n_collision++;
+
+  // Sample the reaction type
+  sample_reaction(p, energy_bin_avg, material_xs);
+
+  // Display information about collision
+  if ((settings::verbosity >= 10) || (openmc_trace)) {
+    std::stringstream msg;
+    msg << "    Energy Group = " << p->g;
+    write_message(msg, 1);
+  }
+}
+
+void
+sample_reaction(Particle* p, const double* energy_bin_avg,
+                const MaterialMacroXS* material_xs)
+{
+  // Create fission bank sites. Note that while a fission reaction is sampled,
+  // it never actually "happens", i.e. the weight of the particle does not
+  // change when sampling fission sites. The following block handles all
+  // absorption (including fission)
+
+  if (materials[p->material - 1]->fissionable) {
+    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+      Bank* result_bank;
+      int64_t result_bank_size;
+      // Get pointer to fission bank from Fortran side
+      openmc_fission_bank(&result_bank, &result_bank_size);
+      create_fission_sites(p, result_bank, &n_bank, result_bank_size,
+                           material_xs);
+    } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
+               (settings::create_fission_neutrons)) {
+      create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
+                           MAX_SECONDARY, material_xs);
+    }
+  }
+
+  // If survival biasing is being used, the following subroutine adjusts the
+  // weight of the particle. Otherwise, it checks to see if absorption occurs.
+  if (material_xs->absorption > 0.) {
+    absorption(p, material_xs);
+  } else {
+    p->absorb_wgt = 0.;
+  }
+  if (!p->alive) return;
+
+  // Sample a scattering event to determine the energy of the exiting neutron
+  scatter(p, energy_bin_avg);
+
+  // Play Russian roulette if survival biasing is turned on
+  if (settings::survival_biasing) {
+    russian_roulette(p);
+    if (!p->alive) return;
+  }
+}
+
+void
+scatter(Particle* p, const double* energy_bin_avg)
+{
+  // Adjust indices for Fortran to C++ indexing
+  // TODO: Remove when no longer needed
+  int gin = p->last_g - 1;
+  int gout = p->g - 1;
+  int i_mat = p->material - 1;
+  macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
+
+  // Adjust return value for fortran indexing
+  // TODO: Remove when no longer needed
+  p->g = gout + 1;
+
+  // Rotate the angle
+  rotate_angle_c(p->coord[0].uvw, p->mu, nullptr);
+
+  // Update energy value for downstream compatability (in tallying)
+  p->E = energy_bin_avg[gout];
+
+  // Set event component
+  p->event = EVENT_SCATTER;
+}
+
+void
+create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
+     int64_t bank_array_size, const MaterialMacroXS* material_xs)
+{
+  // TODO: Heat generation from fission
+
+  // If uniform fission source weighting is turned on, we increase or decrease
+  // the expected number of fission sites produced
+  double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
+
+  // Determine the expected number of neutrons produced
+  double nu_t = p->wgt / openmc_keff * weight * material_xs->nu_fission /
+       material_xs->total;
+
+  // Sample the number of neutrons produced
+  int nu = static_cast<int>(nu_t);
+  if (prn() <= (nu_t - int(nu_t))) {
+    nu++;
+  }
+
+  // Check for the bank size getting hit. For fixed source calculations, this
+  // is a fatal error; for eigenvalue calculations, it just means that k-eff
+  // was too high for a single batch.
+  if (*size_bank + nu > bank_array_size) {
+    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+      throw std::runtime_error{"Secondary particle bank size limit reached."
+           " If you are running a subcritical multiplication problem,"
+           " k-effective may be too close to one."};
+    } else {
+      if (mpi::master) {
+        std::stringstream msg;
+        msg << "Maximum number of sites in fission bank reached. This can"
+             " result in irreproducible results using different numbers of"
+             " processes/threads.";
+        warning(msg);
+      }
+    }
+  }
+
+  // Begin banking the source neutrons
+  // First, if our bank is full then don't continue
+  if ((nu == 0) || (*size_bank == bank_array_size)) return;
+
+  // Initialize the counter of delayed neutrons encountered for each delayed
+  // group.
+  double nu_d[MAX_DELAYED_GROUPS] = {0.};
+
+  p->fission = true;
+  for (size_t i = static_cast<size_t>(*size_bank);
+       i < static_cast<size_t>(std::min(*size_bank + nu, bank_array_size)); i++) {
+    // Bank source neutrons by copying the particle data
+    bank_array[i].xyz[0] = p->coord[0].xyz[0];
+    bank_array[i].xyz[1] = p->coord[0].xyz[1];
+    bank_array[i].xyz[2] = p->coord[0].xyz[2];
+
+    // Set that the bank particle is a neutron
+    bank_array[i].particle = static_cast<int>(ParticleType::neutron);
+
+    // Set the weight of the fission bank site
+    bank_array[i].wgt = 1. / weight;
+
+    // Sample the cosine of the angle, assuming fission neutrons are emitted
+    // isotropically
+    double mu = 2. * prn() - 1.;
+
+    // Sample the azimuthal angle uniformly in [0, 2.pi)
+    double phi = 2. * PI * prn();
+    bank_array[i].uvw[0] = mu;
+    bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi);
+    bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi);
+
+    // Sample secondary energy distribution for the fission reaction and set
+    // the energy in the fission bank
+    int dg;
+    int gout;
+    macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
+    bank_array[i].E = static_cast<double>(gout + 1);
+    bank_array[i].delayed_group = dg + 1;
+
+    // Set the delayed group on the particle as well
+    p->delayed_group = dg + 1;
+
+    // Increment the number of neutrons born delayed
+    if (p->delayed_group > 0) {
+      nu_d[dg]++;
+    }
+  }
+
+  // Increment number of bank sites
+  *size_bank = std::min(*size_bank + nu, bank_array_size);
+
+  // Store the total weight banked for analog fission tallies
+  p->n_bank = nu;
+  p->wgt_bank = nu / weight;
+  for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
+    p->n_delayed_bank[d] = nu_d[d];
+  }
+}
+
+void
+absorption(Particle* p, const MaterialMacroXS* material_xs)
+{
+  if (settings::survival_biasing) {
+    // Determine weight absorbed in survival biasing
+    p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total;
+
+    // Adjust weight of particle by the probability of absorption
+    p->wgt -= p->absorb_wgt;
+    p->last_wgt = p->wgt;
+
+    // Score implicit absorpion estimate of keff
+#pragma omp atomic
+    global_tally_absorption += p->absorb_wgt * material_xs->nu_fission /
+         material_xs->absorption;
+  } else {
+    if (material_xs->absorption > prn() * material_xs->total) {
+#pragma omp atomic
+      global_tally_absorption += p->wgt * material_xs->nu_fission /
+           material_xs->absorption;
+      p->alive = false;
+      p->event = EVENT_ABSORB;
+    }
+
+  }
+}
+
+} //namespace openmc
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index 8d71386243..21e7f9df7a 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -316,8 +316,8 @@ ScattDataLegendre::update_max_val()
         }
 
         // Calculate probability
-        double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
-                                       dist[gin][i_gout].data(), mu);
+        double f = evaluate_legendre(dist[gin][i_gout].size() - 1,
+                                     dist[gin][i_gout].data(), mu);
 
         // if this is a new maximum, store it
         if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f;
@@ -339,8 +339,8 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
     f = 0.;
   } else {
     int i_gout = gout - gmin[gin];
-    f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
-                            dist[gin][i_gout].data(), mu);
+    f = evaluate_legendre(dist[gin][i_gout].size() - 1,
+                          dist[gin][i_gout].data(), mu);
   }
   return f;
 }
@@ -881,8 +881,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
       tab.fmu[gin][i_gout].resize(n_mu);
       for (int imu = 0; imu < n_mu; imu++) {
         tab.fmu[gin][i_gout][imu] =
-             evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1,
-                                 leg.dist[gin][i_gout].data(), tab.mu[imu]);
+             evaluate_legendre(leg.dist[gin][i_gout].size() - 1,
+                               leg.dist[gin][i_gout].data(), tab.mu[imu]);
       }
 
       // Ensure positivity
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 45f6d3cb86..dabd9de3a8 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -32,13 +32,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
   double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
 
   // x is essentially a Maxwellian distribution
-  double x = maxwell_spectrum_c(1.0);
+  double x = maxwell_spectrum(1.0);
 
   double y;
   double r1, r2, r3, r4, r5, r6;
   switch (n_bodies_) {
   case 3:
-    y = maxwell_spectrum_c(1.0);
+    y = maxwell_spectrum(1.0);
     break;
   case 4:
     r1 = prn();
diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90
index db84056812..c6238a88d9 100644
--- a/src/tallies/tally_header.F90
+++ b/src/tallies/tally_header.F90
@@ -137,7 +137,7 @@ module tally_header
   ! invalidation. Thus, we use threadprivate variables to accumulate global
   ! tallies and then reduce at the end of a generation.
   real(C_DOUBLE), public :: global_tally_collision   = ZERO
-  real(C_DOUBLE), public :: global_tally_absorption  = ZERO
+  real(C_DOUBLE), public, bind(C) :: global_tally_absorption  = ZERO
   real(C_DOUBLE), public :: global_tally_tracklength = ZERO
   real(C_DOUBLE), public :: global_tally_leakage     = ZERO
 !$omp threadprivate(global_tally_collision, global_tally_absorption, &
diff --git a/src/tracking.F90 b/src/tracking.F90
index 8cf2fc05c7..94c009be00 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -17,7 +17,6 @@ module tracking
   use nuclide_header
   use particle_header
   use physics,            only: collision
-  use physics_mg,         only: collision_mg
   use random_lcg,         only: prn, prn_set_stream
   use settings
   use simulation_header
@@ -33,6 +32,16 @@ module tracking
 
   implicit none
 
+  interface
+    subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C)
+      import Particle, C_DOUBLE, MaterialMacroXS
+      type(Particle),            intent(inout) :: p
+      real(C_DOUBLE),            intent(in)    :: energy_bin_avg(*)
+      type(MaterialMacroXS),     intent(in)    :: material_xs
+    end subroutine collision_mg
+
+  end interface
+
 contains
 
 !===============================================================================
@@ -226,7 +235,7 @@ contains
         if (run_CE) then
           call collision(p)
         else
-          call collision_mg(p)
+          call collision_mg(p, energy_bin_avg, material_xs)
         end if
 
         ! Score collision estimator tallies -- this is done after a collision