diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index 1ee9ceba0..6deb6c01d 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -98,10 +98,7 @@ public:
   Nuclide(hid_t group, const double* temperature, int n, int i_nuclide);
 
   //! Initialize logarithmic grid for energy searches
-  //! \param E_min Minimum energy in [eV]
-  //! \param E_max Maximum energy in [eV]
-  //! \param M Number of equally log-spaced bins
-  void init_grid(double E_min, double E_max, int M);
+  void init_grid();
 
   void calculate_xs(int i_sab, double E, int i_log_union,
     double sqrtkT, double sab_frac);
@@ -160,11 +157,11 @@ public:
 private:
   void create_derived();
 
-  static constexpr int XS_TOTAL {0};
-  static constexpr int XS_ABSORPTION {1};
-  static constexpr int XS_FISSION {2};
-  static constexpr int XS_NU_FISSION {3};
-  static constexpr int XS_PHOTON_PROD {4};
+  static int XS_TOTAL;
+  static int XS_ABSORPTION;
+  static int XS_FISSION;
+  static int XS_NU_FISSION;
+  static int XS_PHOTON_PROD;
 };
 
 //==============================================================================
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index 66596b76d..6e9a947d0 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -137,18 +137,6 @@ public:
   std::vector<ThermalData> data_;
 };
 
-//==============================================================================
-// Fortran compatibility functions
-//==============================================================================
-
-extern "C" {
-  ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature,
-    int n, int method, double tolerance, const double* minmax);
-  void sab_free(ThermalScattering* data);
-  bool sab_has_nuclide(ThermalScattering* data, const char* name);
-  double sab_threshold(ThermalScattering* data);
-}
-
 } // namespace openmc
 
 #endif // OPENMC_THERMAL_H
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index f0796c3f6..8e90ae03e 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2157,6 +2157,9 @@ contains
     end if
 
     ! Set up logarithmic grid for nuclides
+    do i = 1, size(nuclides)
+      call nuclides(i) % init_grid()
+    end do
     log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins
 
     do i = 1, size(materials)
diff --git a/src/material_header.F90 b/src/material_header.F90
index 13914cc59..5daf9c42c 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -443,7 +443,7 @@ contains
       ! Calculate microscopic cross section for this nuclide
       if (p % E /= micro_xs(i_nuclide) % last_E &
            .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
-           .or. i_sab /= micro_xs(i_nuclide) % index_sab &
+           .or. i_sab /= micro_xs(i_nuclide) % index_sab + 1 &
            .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
         call nuclides(i_nuclide) % calculate_xs(i_sab, p % E, i_grid, &
              p % sqrtkT, sab_frac)
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 011364588..753c0af67 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -39,7 +39,14 @@ MaterialMacroXS material_xs;
 // Nuclide implementation
 //==============================================================================
 
+int Nuclide::XS_TOTAL {0};
+int Nuclide::XS_ABSORPTION {1};
+int Nuclide::XS_FISSION {2};
+int Nuclide::XS_NU_FISSION {3};
+int Nuclide::XS_PHOTON_PROD {4};
+
 Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
+  : i_nuclide_{i_nuclide}
 {
   // Get name of nuclide from group, removing leading '/'
   name_ = object_name(group).substr(1);
@@ -48,7 +55,6 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
   read_attribute(group, "A", A_);
   read_attribute(group, "metastable", metastable_);
   read_attribute(group, "atomic_weight_ratio", awr_);
-  i_nuclide_ = i_nuclide;
 
   // Determine temperatures available
   hid_t kT_group = open_group(group, "kTs");
@@ -382,8 +388,13 @@ void Nuclide::create_derived()
   }
 }
 
-void Nuclide::init_grid(double E_min, double E_max, int M)
+void Nuclide::init_grid()
 {
+  int neutron = static_cast<int>(ParticleType::neutron) - 1;
+  double E_min = data::energy_min[neutron];
+  double E_max = data::energy_max[neutron];
+  int M = settings::n_log_bins;
+
   // Determine equal-logarithmic energy spacing
   double spacing = std::log(E_max/E_min)/M;
 
@@ -398,10 +409,10 @@ void Nuclide::init_grid(double E_min, double E_max, int M)
     // equal-logarithmic grid
     int j = 0;
     for (int k = 0; k <= M; ++k) {
-      while (std::log(grid.energy[j]/E_min) <= umesh(k)) {
+      while (std::log(grid.energy[j + 1]/E_min) <= umesh(k)) {
         // Ensure that for isotopes where maxval(grid.energy) << E_max that
         // there are no out-of-bounds issues.
-        if (j == grid.energy.size()) break;
+        if (j + 1 == grid.energy.size()) break;
         ++j;
       }
       grid.grid_index[k] = j;
@@ -453,7 +464,7 @@ void Nuclide::calculate_elastic_xs() const
 {
   // Get temperature index, grid index, and interpolation factor
   auto& micro = simulation::micro_xs[i_nuclide_];
-  int i_temp = micro.index_temp - 1;
+  int i_temp = micro.index_temp;
   int i_grid = micro.index_grid - 1;
   double f = micro.interp_factor;
 
@@ -589,7 +600,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
       int i_high = grid.grid_index[i_log_union + 1] + 1;
 
       // Perform binary search over reduced range
-      i_grid = lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E) - 1;
+      i_grid = i_low + lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E);
     }
 
     // check for rare case where two energy points are the same
@@ -600,7 +611,8 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
       (grid.energy[i_grid + 1]- grid.energy[i_grid]);
 
     micro_xs.index_temp = i_temp;
-    micro_xs.index_grid = i_grid;
+    // TODO: off-by-one
+    micro_xs.index_grid = i_grid + 1;
     micro_xs.interp_factor = f;
 
     // Calculate microscopic nuclide total cross section
@@ -871,6 +883,8 @@ extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature,
   return data::nuclides.back().get();
 }
 
+extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); }
+
 extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx)
 {
   return nuc->reactions_[i_rx-1].get();
@@ -882,10 +896,15 @@ extern "C" void nuclide_calculate_xs_c(Nuclide* nuc, int i_sab, double E,
   nuc->calculate_xs(i_sab, E, i_log_union, sqrtkT, sab_frac);
 }
 
+extern "C" void nuclide_calculate_elastic_xs_c(Nuclide* nuc)
+{
+  nuc->calculate_elastic_xs();
+}
+
 extern "C" void nuclides_clear() { data::nuclides.clear(); }
 
-
 extern "C" NuclideMicroXS* micro_xs_ptr();
+
 void set_micro_xs()
 {
 #pragma omp parallel
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 04c543522..d6b564c82 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -94,6 +94,7 @@ module nuclide_header
   contains
     procedure :: clear => nuclide_clear
     procedure :: from_hdf5 => nuclide_from_hdf5
+    procedure :: init_grid => nuclide_init_grid
     procedure :: nu    => nuclide_nu
     procedure, private :: create_derived => nuclide_create_derived
     procedure :: calculate_xs => nuclide_calculate_xs
@@ -673,6 +674,17 @@ contains
 
   end function nuclide_nu
 
+  subroutine nuclide_init_grid(this)
+    class(nuclide), intent(inout) :: this
+    interface
+      subroutine nuclide_init_grid_c(ptr) bind(C)
+        import C_PTR
+        type(C_PTR), value :: ptr
+      end subroutine
+    end interface
+    call nuclide_init_grid_c(this % ptr)
+  end subroutine
+
 !===============================================================================
 ! NUCLIDE_CALCULATE_XS determines microscopic cross sections for the nuclide
 ! at the energy of the given particle
@@ -689,8 +701,9 @@ contains
     real(8), intent(in) :: sab_frac    ! fraction of atoms affected by S(a,b)
 
     interface
-      subroutine nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C)
-        import C_INT, C_DOUBLE
+      subroutine nuclide_calculate_xs_c(ptr, i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C)
+        import C_PTR, C_INT, C_DOUBLE
+        type(C_PTR), value :: ptr
         integer(C_INT), value :: i_sab
         real(C_DOUBLE), value :: E
         integer(C_INT), value :: i_log_union
@@ -699,7 +712,7 @@ contains
       end subroutine
     end interface
 
-    call nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac)
+    call nuclide_calculate_xs_c(this % ptr, i_sab - 1, E, i_log_union, sqrtkT, sab_frac)
   end subroutine nuclide_calculate_xs
 
 !===============================================================================
@@ -709,25 +722,15 @@ contains
 ! to calculate it early to adjust the total cross section correctly.
 !===============================================================================
 
-  subroutine nuclide_calculate_elastic_xs(this, micro_xs)
+  subroutine nuclide_calculate_elastic_xs(this)
     class(Nuclide),       intent(in)    :: this
-    type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
-
-    integer :: i_temp
-    integer :: i_grid
-    real(8) :: f
-
-    ! Get temperature index, grid index, and interpolation factor
-    i_temp =  micro_xs % index_temp
-    i_grid =  micro_xs % index_grid
-    f      =  micro_xs % interp_factor
-
-    if (i_temp > 0) then
-      associate (rx => this % reactions(1))
-        micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + &
-             f * rx % xs(i_temp, i_grid + 1)
-      end associate
-    end if
+    interface
+      subroutine nuclide_calculate_elastic_xs_c(ptr) bind(C)
+        import C_PTR
+        type(C_PTR), value :: ptr
+      end subroutine
+    end interface
+    call nuclide_calculate_elastic_xs_c(this % ptr)
   end subroutine nuclide_calculate_elastic_xs
 
 !===============================================================================
diff --git a/src/physics.cpp b/src/physics.cpp
index e3af8bdef..5c8d2ccf3 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -518,7 +518,7 @@ Reaction* sample_fission(int i_nuclide, double E)
   }
 
   // Get grid index and interpolatoin factor and sample fission cdf
-  int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1;
+  int i_temp = simulation::micro_xs[i_nuclide].index_temp;
   int i_grid = simulation::micro_xs[i_nuclide].index_grid;
   double f = simulation::micro_xs[i_nuclide].interp_factor;
   double cutoff = prn() * simulation::micro_xs[i_nuclide].fission;
@@ -542,8 +542,7 @@ Reaction* sample_fission(int i_nuclide, double E)
 void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product)
 {
   // Get grid index and interpolation factor and sample photon production cdf
-  // TODO: off-by-one
-  int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1;
+  int i_temp = simulation::micro_xs[i_nuclide].index_temp;
   int i_grid = simulation::micro_xs[i_nuclide].index_grid;
   double f = simulation::micro_xs[i_nuclide].interp_factor;
   double cutoff = prn() * simulation::micro_xs[i_nuclide].photon_prod;
@@ -616,7 +615,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
   // Get pointer to nuclide and grid index/interpolation factor
   const auto& nuc {data::nuclides[i_nuclide]};
   const auto& micro {simulation::micro_xs[i_nuclide]};
-  int i_temp =  micro.index_temp - 1;
+  int i_temp =  micro.index_temp;
   int i_grid =  micro.index_grid - 1;
   double f = micro.interp_factor;
 
diff --git a/src/sab_header.F90 b/src/sab_header.F90
index 061da5d1c..b1e49d293 100644
--- a/src/sab_header.F90
+++ b/src/sab_header.F90
@@ -21,7 +21,6 @@ module sab_header
     type(C_PTR) :: ptr
   contains
     procedure :: from_hdf5
-    procedure :: free
     procedure :: has_nuclide
     procedure :: threshold
   end type SAlphaBeta
@@ -44,9 +43,7 @@ module sab_header
       type(C_PTR) :: ptr
     end function
 
-    subroutine sab_free(ptr) bind(C)
-      import C_PTR
-      type(C_PTR), value :: ptr
+    subroutine sab_clear() bind(C)
     end subroutine
 
     function sab_has_nuclide(ptr, name) result(val) bind(C)
@@ -87,11 +84,6 @@ contains
     end if
   end subroutine from_hdf5
 
-  subroutine free(this)
-    class(SAlphaBeta), intent(inout) :: this
-    call sab_free(this % ptr)
-  end subroutine
-
   function has_nuclide(this, name) result(val)
     class(SAlphaBeta), intent(in) :: this
     character(len=*), intent(in) :: name
@@ -111,14 +103,9 @@ contains
 !===============================================================================
 
   subroutine free_memory_sab()
-    integer :: i
     n_sab_tables = 0
-    if (allocated(sab_tables)) then
-      do i = 1, size(sab_tables)
-        call sab_tables(i) % free()
-      end do
-      deallocate(sab_tables)
-    end if
+    call sab_clear()
+    if (allocated(sab_tables)) deallocate(sab_tables)
     call sab_dict % clear()
   end subroutine free_memory_sab
 
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 1889d0dcb..e014e7402 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -958,7 +958,7 @@ contains
         else
           if (i_nuclide > 0) then
             if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
-              call nuclides(i_nuclide) % calculate_elastic_xs(micro_xs(i_nuclide))
+              call nuclides(i_nuclide) % calculate_elastic_xs()
             end if
             score = micro_xs(i_nuclide) % elastic * atom_density * flux
           else
@@ -971,7 +971,7 @@ contains
                 ! Get index in nuclides array
                 i_nuc = materials(p % material) % nuclide(l)
                 if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then
-                  call nuclides(i_nuc) % calculate_elastic_xs(micro_xs(i_nuc))
+                  call nuclides(i_nuc) % calculate_elastic_xs()
                 end if
 
                 score = score + micro_xs(i_nuc) % elastic * atom_density_ * flux
@@ -1126,7 +1126,7 @@ contains
               if (m /= 0) then
                 ! Retrieve temperature and energy grid index and interpolation
                 ! factor
-                i_temp = micro_xs(i_nuclide) % index_temp
+                i_temp = micro_xs(i_nuclide) % index_temp + 1
                 if (i_temp > 0) then
                   i_energy = micro_xs(i_nuclide) % index_grid
                   f = micro_xs(i_nuclide) % interp_factor
@@ -1159,7 +1159,7 @@ contains
                   if (m /= 0) then
                     ! Retrieve temperature and energy grid index and
                     ! interpolation factor
-                    i_temp = micro_xs(i_nuc) % index_temp
+                    i_temp = micro_xs(i_nuc) % index_temp + 1
                     if (i_temp > 0) then
                       i_energy = micro_xs(i_nuc) % index_grid
                       f = micro_xs(i_nuc) % interp_factor
diff --git a/src/thermal.cpp b/src/thermal.cpp
index 48780e4eb..7fed2c0f1 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -590,7 +590,7 @@ ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
 // Fortran compatibility functions
 //==============================================================================
 
-ThermalScattering*
+extern "C" ThermalScattering*
 sab_from_hdf5(hid_t group, const double* temperature, int n,
     int method, double tolerance, const double* minmax)
 {
@@ -604,13 +604,16 @@ sab_from_hdf5(hid_t group, const double* temperature, int n,
   return data::thermal_scatt.back().get();
 }
 
-void sab_free(ThermalScattering* data) { delete data; }
+extern "C" void sab_clear() { data::thermal_scatt.clear(); }
 
-bool sab_has_nuclide(ThermalScattering* data, const char* name)
+extern "C" bool sab_has_nuclide(ThermalScattering* data, const char* name)
 {
   return data->has_nuclide(name);
 }
 
-double sab_threshold(ThermalScattering* data) { return data->threshold(); }
+extern "C" double sab_threshold(ThermalScattering* data)
+{
+  return data->threshold();
+}
 
 } // namespace openmc