diff --git a/src/material.cpp b/src/material.cpp
index 51b273ca9..510cea966 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -636,7 +636,7 @@ void Material::init_nuclide_index()
   int n = settings::run_CE ?
     data::nuclides.size() : data::nuclides_MG.size();
   mat_nuclide_index_.resize(n);
-  std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), -1);
+  std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
   for (int i = 0; i < nuclide_.size(); ++i) {
     mat_nuclide_index_[nuclide_[i]] = i;
   }
@@ -1135,11 +1135,6 @@ extern "C" {
     return model::materials[i_mat - 1]->density_gpcc_;
   }
 
-  int material_nuclide_index(int32_t i_mat, int i_nuc)
-  {
-    return model::materials[i_mat - 1]->mat_nuclide_index_[i_nuc - 1] + 1;
-  }
-
   void material_calculate_xs(const Particle* p)
   {
     model::materials[p->material - 1]->calculate_xs(*p);
diff --git a/src/material_header.F90 b/src/material_header.F90
index 81baeb861..cda79ec60 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -11,7 +11,6 @@ module material_header
   public :: material_id
   public :: material_nuclide
   public :: material_nuclide_size
-  public :: material_nuclide_index
   public :: material_atom_density
   public :: material_density_gpcc
 
@@ -40,13 +39,6 @@ module material_header
       integer(C_INT) :: n
     end function
 
-    function material_nuclide_index(i_mat, i_nuc) bind(C) result(idx)
-      import C_INT32_T, C_INT
-      integer(C_INT32_T), value :: i_mat
-      integer(C_INT), value :: i_nuc
-      integer(C_INT) :: idx
-    end function
-
     function material_atom_density(i_mat, idx) bind(C) result(density)
       import C_INT32_T, C_INT, C_DOUBLE
       integer(C_INT32_T), value :: i_mat
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 1199f6ba6..febabb0e3 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -3,6 +3,7 @@
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/error.h"
+#include "openmc/material.h"
 #include "openmc/message_passing.h"
 #include "openmc/nuclide.h"
 #include "openmc/settings.h"
@@ -44,10 +45,6 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index,
 extern "C" void
 score_all_nuclides(Particle* p, int i_tally, double flux, int filter_index);
 
-extern "C" int material_nuclide_index(int i_material, int i_nuclide);
-
-extern "C" double material_atom_density(int i_material, int i);
-
 //==============================================================================
 // Global variable definitions
 //==============================================================================
@@ -729,9 +726,13 @@ score_analog_tally_mg(Particle* p)
 
         double atom_density = 0.;
         if (i_nuclide > 0) {
-          auto j = material_nuclide_index(p->material, i_nuclide);
-          if (j == 0) continue;
-          atom_density = material_atom_density(p->material, j);
+          //TODO: off-by-one
+          auto j = model::materials[p->material-1]
+            ->mat_nuclide_index_[i_nuclide-1];
+          if (j == C_NONE) continue;
+          //atom_density = material_atom_density(p->material, j);
+          //TODO: off-by-one
+          atom_density = model::materials[p->material-1]->atom_density_(j);
         }
 
         score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
@@ -791,9 +792,13 @@ score_tracklength_tally(Particle* p, double distance)
           double atom_density = 0.;
           if (i_nuclide > 0) {
             if (p->material != MATERIAL_VOID) {
-              auto j = material_nuclide_index(p->material, i_nuclide);
-              if (j == 0) continue;
-              atom_density = material_atom_density(p->material, j);
+              //TODO: off-by-one
+              auto j = model::materials[p->material-1]
+                ->mat_nuclide_index_[i_nuclide-1];
+              if (j == C_NONE) continue;
+              //atom_density = material_atom_density(p->material, j);
+              //TODO: off-by-one
+              atom_density = model::materials[p->material-1]->atom_density_(j);
             }
           }
 
@@ -867,9 +872,13 @@ score_collision_tally(Particle* p)
 
           double atom_density = 0.;
           if (i_nuclide > 0) {
-            auto j = material_nuclide_index(p->material, i_nuclide);
-            if (j == 0) continue;
-            atom_density = material_atom_density(p->material, j);
+            //TODO: off-by-one
+            auto j = model::materials[p->material-1]
+              ->mat_nuclide_index_[i_nuclide-1];
+            if (j == C_NONE) continue;
+            //atom_density = material_atom_density(p->material, j);
+            //TODO: off-by-one
+            atom_density = model::materials[p->material-1]->atom_density_(j);
           }
 
           //TODO: consider replacing this "if" with pointers or templates