diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py
index 894e4c79e..7ad5c75d7 100644
--- a/openmc/capi/cell.py
+++ b/openmc/capi/cell.py
@@ -107,12 +107,9 @@ class Cell(_FortranObjectWithID):
 
         if fill_type.value == 1:
             if n.value > 1:
-                #TODO: off-by-one
-                return [Material(index=i+1 if i >= 0 else i)
-                        for i in indices[:n.value]]
+                return [Material(index=i) for i in indices[:n.value]]
             else:
-                #TODO: off-by-one
-                index = indices[0] + 1 if indices[0] >= 0 else indices[0]
+                index = indices[0]
                 return Material(index=index)
         else:
             raise NotImplementedError
diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py
index 0f5631519..0b85e026e 100644
--- a/openmc/capi/filter.py
+++ b/openmc/capi/filter.py
@@ -233,16 +233,13 @@ class MaterialFilter(Filter):
         materials = POINTER(c_int32)()
         n = c_int32()
         _dll.openmc_material_filter_get_bins(self._index, materials, n)
-        #TODO: fix this off-by-one when materials become 0-indexed
-        return [Material(index=materials[i]+1) for i in range(n.value)]
+        return [Material(index=materials[i]) for i in range(n.value)]
 
     @bins.setter
     def bins(self, materials):
         # Get material indices as int32_t[]
         n = len(materials)
-        #TODO: fix this off-by-one when materials become 0-indexed
-        bins = (c_int32*n)(*(m._index-1 for m in materials))
-
+        bins = (c_int32*n)(*(m._index for m in materials))
         _dll.openmc_material_filter_set_bins(self._index, n, bins)
 
 
diff --git a/openmc/capi/material.py b/openmc/capi/material.py
index 4bf19a1e0..981185338 100644
--- a/openmc/capi/material.py
+++ b/openmc/capi/material.py
@@ -227,7 +227,7 @@ class _MaterialMapping(Mapping):
 
     def __iter__(self):
         for i in range(len(self)):
-            yield Material(index=i + 1).id
+            yield Material(index=i).id
 
     def __len__(self):
         return _dll.n_materials()
diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp
index f0ece1a53..2b44154da 100644
--- a/src/bremsstrahlung.cpp
+++ b/src/bremsstrahlung.cpp
@@ -36,9 +36,9 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
   // Get bremsstrahlung data for this material and particle type
   BremsstrahlungData* mat;
   if (p.type == static_cast<int>(ParticleType::positron)) {
-    mat = &model::materials[p.material -1]->ttb_->positron;
+    mat = &model::materials[p.material]->ttb_->positron;
   } else {
-    mat = &model::materials[p.material -1]->ttb_->electron;
+    mat = &model::materials[p.material]->ttb_->electron;
   }
 
   double e = std::log(p.E);
diff --git a/src/cell.cpp b/src/cell.cpp
index 30ca4b7cf..2608d8f6f 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -727,9 +727,8 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
         int i_mat = indices[i];
         if (i_mat == MATERIAL_VOID) {
           c.material_.push_back(MATERIAL_VOID);
-        } else if (i_mat >= 1 && i_mat <= model::materials.size()) {
-          //TODO: off-by-one
-          c.material_.push_back(i_mat - 1);
+        } else if (i_mat >= 0 && i_mat < model::materials.size()) {
+          c.material_.push_back(i_mat);
         } else {
           set_errmsg("Index in materials array is out of bounds.");
           return OPENMC_E_OUT_OF_BOUNDS;
diff --git a/src/geometry.cpp b/src/geometry.cpp
index 38dba164b..fb081c313 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -147,16 +147,10 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
 
       // Set the material and temperature.
       p->last_material = p->material;
-      int32_t mat;
       if (c.material_.size() > 1) {
-        mat = c.material_[p->cell_instance];
+        p->material = c.material_[p->cell_instance];
       } else {
-        mat = c.material_[0];
-      }
-      if (mat == MATERIAL_VOID) {
-        p->material = MATERIAL_VOID;
-      } else {
-        p->material = mat + 1;
+        p->material = c.material_[0];
       }
       p->last_sqrtkT = p->sqrtkT;
       if (c.sqrtkT_.size() > 1) {
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index d38b5f8f5..c118c04ee 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -85,18 +85,17 @@ adjust_indices()
       }
     } else {
       c->type_ = FILL_MATERIAL;
-      for (auto it = c->material_.begin(); it != c->material_.end(); it++) {
-        int32_t mid = *it;
-        if (mid != MATERIAL_VOID) {
-          auto search = model::material_map.find(mid);
-          if (search != model::material_map.end()) {
-            *it = search->second;
-          } else {
+      for (auto& mat_id : c->material_) {
+        if (mat_id != MATERIAL_VOID) {
+          auto search = model::material_map.find(mat_id);
+          if (search == model::material_map.end()) {
             std::stringstream err_msg;
-            err_msg << "Could not find material " << mid
+            err_msg << "Could not find material " << mat_id
                     << " specified on cell " << c->id_;
             fatal_error(err_msg);
           }
+          // Change from ID to index
+          mat_id = search->second;
         }
       }
     }
diff --git a/src/material.cpp b/src/material.cpp
index 3153d585a..8236a4ef0 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -932,7 +932,7 @@ openmc_get_material_index(int32_t id, int32_t* index)
     set_errmsg("No material exists with ID=" + std::to_string(id) + ".");
     return OPENMC_E_INVALID_ID;
   } else {
-    *index = it->second + 1;
+    *index = it->second;
     return 0;
   }
 }
@@ -941,8 +941,8 @@ extern "C" int
 openmc_material_add_nuclide(int32_t index, const char* name, double density)
 {
   int err = 0;
-  if (index >= 1 && index <= model::materials.size()) {
-    Material* m = model::materials[index - 1];
+  if (index >= 0 && index < model::materials.size()) {
+    Material* m = model::materials[index];
 
     // Check if nuclide is already in material
     for (int i = 0; i < m->nuclide_.size(); ++i) {
@@ -987,8 +987,8 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density)
 extern "C" int
 openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    auto& mat = model::materials[index - 1];
+  if (index >= 0 && index < model::materials.size()) {
+    auto& mat = model::materials[index];
     if (!mat->nuclide_.empty()) {
       *nuclides = mat->nuclide_.data();
       *densities = mat->atom_density_.data();
@@ -1007,8 +1007,8 @@ openmc_material_get_densities(int32_t index, int** nuclides, double** densities,
 extern "C" int
 openmc_material_get_fissionable(int32_t index, bool* fissionable)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    *fissionable = model::materials[index - 1]->fissionable_;
+  if (index >= 0 && index < model::materials.size()) {
+    *fissionable = model::materials[index]->fissionable_;
     return 0;
   } else {
     set_errmsg("Index in materials array is out of bounds.");
@@ -1019,8 +1019,8 @@ openmc_material_get_fissionable(int32_t index, bool* fissionable)
 extern "C" int
 openmc_material_get_id(int32_t index, int32_t* id)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    *id = model::materials[index - 1]->id_;
+  if (index >= 0 && index < model::materials.size()) {
+    *id = model::materials[index]->id_;
     return 0;
   } else {
     set_errmsg("Index in materials array is out of bounds.");
@@ -1031,8 +1031,8 @@ openmc_material_get_id(int32_t index, int32_t* id)
 extern "C" int
 openmc_material_get_volume(int32_t index, double* volume)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    Material* m = model::materials[index - 1];
+  if (index >= 0 && index < model::materials.size()) {
+    Material* m = model::materials[index];
     if (m->volume_ >= 0.0) {
       *volume = m->volume_;
       return 0;
@@ -1051,8 +1051,8 @@ openmc_material_get_volume(int32_t index, double* volume)
 extern "C" int
 openmc_material_set_density(int32_t index, double density, const char* units)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    return model::materials[index - 1]->set_density(density, units);
+  if (index >= 0 && index < model::materials.size()) {
+    return model::materials[index]->set_density(density, units);
   } else {
     set_errmsg("Index in materials array is out of bounds.");
     return OPENMC_E_OUT_OF_BOUNDS;
@@ -1062,9 +1062,8 @@ openmc_material_set_density(int32_t index, double density, const char* units)
 extern "C" int
 openmc_material_set_densities(int32_t index, int n, const char** name, const double* density)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    // TODO: off-by-one
-    auto& mat {model::materials[index - 1]};
+  if (index >= 0 && index < model::materials.size()) {
+    auto& mat {model::materials[index]};
     if (n != mat->nuclide_.size()) {
       mat->nuclide_.resize(n);
       mat->atom_density_ = xt::zeros<double>({n});
@@ -1099,9 +1098,9 @@ openmc_material_set_densities(int32_t index, int n, const char** name, const dou
 extern "C" int
 openmc_material_set_id(int32_t index, int32_t id)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    model::materials[index - 1]->id_ = id;
-    model::material_map[id] = index - 1;
+  if (index >= 0 && index < model::materials.size()) {
+    model::materials[index]->id_ = id;
+    model::material_map[id] = index;
     return 0;
   } else {
     set_errmsg("Index in materials array is out of bounds.");
@@ -1112,8 +1111,8 @@ openmc_material_set_id(int32_t index, int32_t id)
 extern "C" int
 openmc_material_set_volume(int32_t index, double volume)
 {
-  if (index >= 1 && index <= model::materials.size()) {
-    Material* m = model::materials[index - 1];
+  if (index >= 0 && index < model::materials.size()) {
+    Material* m = model::materials[index];
     if (volume >= 0.0) {
       m->volume_ = volume;
       return 0;
@@ -1130,9 +1129,8 @@ openmc_material_set_volume(int32_t index, double volume)
 extern "C" int
 openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end)
 {
-  // TODO: off-by-one
-  if (index_start) *index_start = model::materials.size() + 1;
-  if (index_end) *index_end = model::materials.size() + n;
+  if (index_start) *index_start = model::materials.size();
+  if (index_end) *index_end = model::materials.size() + n - 1;
   for (int32_t i = 0; i < n; i++) {
     model::materials.push_back(new Material());
   }
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 402799c50..581e6b178 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -251,7 +251,7 @@ 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)
 {
-  data::macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
+  data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
        nu_fiss_xs);
 }
 
@@ -302,7 +302,7 @@ get_macro_xs(int index, int xstype, int gin, const int* gout,
   } else {
     dg_c_p = dg;
   }
-  return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
+  return data::macro_xs[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
 }
 
 //==============================================================================
diff --git a/src/particle.cpp b/src/particle.cpp
index bcd981b26..ccce1502f 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -84,8 +84,8 @@ Particle::initialize()
   // clear attributes
   surface           = 0;
   cell_born         = C_NONE;
-  material          = 0;
-  last_material     = 0;
+  material          = C_NONE;
+  last_material     = C_NONE;
   last_sqrtkT       = 0;
   wgt               = 1.0;
   last_wgt          = 1.0;
@@ -201,7 +201,7 @@ Particle::transport()
           // If the material is the same as the last material and the
           // temperature hasn't changed, we don't need to lookup cross
           // sections again.
-          model::materials[material - 1]->calculate_xs(*this);
+          model::materials[material]->calculate_xs(*this);
         }
       } else {
         // Get the MG data
@@ -346,7 +346,7 @@ Particle::transport()
 
       // Set last material to none since cross sections will need to be
       // re-evaluated
-      last_material = F90_NONE;
+      last_material = C_NONE;
 
       // Set all uvws to base level -- right now, after a collision, only the
       // base level uvws are changed
@@ -584,9 +584,7 @@ Particle::cross_surface()
     // set new cell value
     coord[0].cell = i_cell;
     cell_instance = 0;
-    // TODO: off-by-one
-    int mat = model::cells[i_cell]->material_[0];
-    material = (mat == MATERIAL_VOID) ? mat : mat + 1;
+    material = model::cells[i_cell]->material_[0];
     sqrtkT = model::cells[i_cell]->sqrtkT_[0];
     return;
   }
diff --git a/src/physics.cpp b/src/physics.cpp
index 4d0fbf27b..8b65fe852 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -425,8 +425,7 @@ int sample_nuclide(const Particle* p)
   double cutoff = prn() * simulation::material_xs.total;
 
   // Get pointers to nuclide/density arrays
-  // TODO: off-by-one
-  const auto& mat {model::materials[p->material - 1]};
+  const auto& mat {model::materials[p->material]};
   int n = mat->nuclide_.size();
 
   double prob = 0.0;
@@ -451,7 +450,7 @@ int sample_element(Particle* p)
   double cutoff = prn() * simulation::material_xs.total;
 
   // Get pointers to elements, densities
-  const auto& mat {model::materials[p->material - 1]};
+  const auto& mat {model::materials[p->material]};
   int n = mat->nuclide_.size();
 
   int i = 0;
@@ -675,7 +674,7 @@ void scatter(Particle* p, int i_nuclide)
 
   // Sample new outgoing angle for isotropic-in-lab scattering
   // TODO: off-by-one
-  const auto& mat {model::materials[p->material - 1]};
+  const auto& mat {model::materials[p->material]};
   if (!mat->p0_.empty()) {
     int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
     if (mat->p0_[i_nuc_mat]) {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index ab5232c2a..87dd37ef8 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -46,7 +46,7 @@ sample_reaction(Particle* p)
   // change when sampling fission sites. The following block handles all
   // absorption (including fission)
 
-  if (model::materials[p->material - 1]->fissionable_) {
+  if (model::materials[p->material]->fissionable_) {
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
       create_fission_sites(
         p, simulation::fission_bank.data(), &simulation::n_bank,
@@ -84,7 +84,7 @@ scatter(Particle* p)
   // TODO: Remove when no longer needed
   int gin = p->last_g - 1;
   int gout = p->g - 1;
-  int i_mat = p->material - 1;
+  int i_mat = p->material;
   data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
 
   // Adjust return value for fortran indexing
@@ -176,7 +176,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
     // the energy in the fission bank
     int dg;
     int gout;
-    data::macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
+    data::macro_xs[p->material].sample_fission_energy(p->g - 1, dg, gout);
     bank_array[i].E = static_cast<double>(gout + 1);
     bank_array[i].delayed_group = dg + 1;
 
diff --git a/src/plot.cpp b/src/plot.cpp
index ec486ae27..6f8ddd5a6 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -676,8 +676,8 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
         rgb = WHITE;
         id = -1;
       } else {
-        rgb = pl.colors_[p.material - 1];
-        id = model::materials[p.material - 1]->id_;
+        rgb = pl.colors_[p.material];
+        id = model::materials[p.material]->id_;
       }
     } else if (PlotColorBy::cells == pl.color_by_) {
       // Assign color based on cell
diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp
index b5b85fd6f..1944b39dd 100644
--- a/src/tallies/derivative.cpp
+++ b/src/tallies/derivative.cpp
@@ -131,8 +131,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
     score *= flux_deriv;
     return;
   }
-  //TODO: off-by-one
-  const Material& material {*model::materials[p->material-1]};
+  const Material& material {*model::materials[p->material]};
   if (material.id_ != deriv.diff_material) {
     score *= flux_deriv;
     return;
@@ -574,8 +573,7 @@ score_track_derivative(const Particle* p, double distance)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material == MATERIAL_VOID) return;
-  //TODO: off-by-one
-  const Material& material {*model::materials[p->material-1]};
+  const Material& material {*model::materials[p->material]};
 
   for (auto& deriv : model::tally_derivs) {
     if (deriv.diff_material != material.id_) continue;
@@ -622,8 +620,7 @@ void score_collision_derivative(const Particle* p)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material == MATERIAL_VOID) return;
-  //TODO: off-by-one
-  const Material& material {*model::materials[p->material-1]};
+  const Material& material {*model::materials[p->material]};
 
   for (auto& deriv : model::tally_derivs) {
     if (deriv.diff_material != material.id_) continue;
diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp
index 05842b056..0dee964de 100644
--- a/src/tallies/filter_material.cpp
+++ b/src/tallies/filter_material.cpp
@@ -42,7 +42,7 @@ void
 MaterialFilter::get_all_bins(const Particle* p, int estimator,
                              FilterMatch& match) const
 {
-  auto search = map_.find(p->material - 1);
+  auto search = map_.find(p->material);
   if (search != map_.end()) {
     //TODO: off-by-one
     match.bins_.push_back(search->second + 1);
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 9edac5ec2..cdd80adba 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -600,7 +600,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           score = 0.;
           // Add up contributions from each nuclide in the material.
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -723,7 +723,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
               {*dynamic_cast<DelayedGroupFilter*>(
               model::tally_filters[i_dg_filt].get())};
             if (p->material != MATERIAL_VOID) {
-              const Material& material {*model::materials[p->material-1]};
+              const Material& material {*model::materials[p->material]};
               for (auto i = 0; i < material.nuclide_.size(); ++i) {
                 auto j_nuclide = material.nuclide_[i];
                 auto atom_density = material.atom_density_(i);
@@ -743,7 +743,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           } else {
             score = 0.;
             if (p->material != MATERIAL_VOID) {
-              const Material& material {*model::materials[p->material-1]};
+              const Material& material {*model::materials[p->material]};
               for (auto i = 0; i < material.nuclide_.size(); ++i) {
                 auto j_nuclide = material.nuclide_[i];
                 auto atom_density = material.atom_density_(i);
@@ -890,7 +890,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
               {*dynamic_cast<DelayedGroupFilter*>(
               model::tally_filters[i_dg_filt].get())};
             if (p->material != MATERIAL_VOID) {
-              const Material& material {*model::materials[p->material-1]};
+              const Material& material {*model::materials[p->material]};
               for (auto i = 0; i < material.nuclide_.size(); ++i) {
                 auto j_nuclide = material.nuclide_[i];
                 auto atom_density = material.atom_density_(i);
@@ -915,7 +915,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           } else {
             score = 0.;
             if (p->material != MATERIAL_VOID) {
-              const Material& material {*model::materials[p->material-1]};
+              const Material& material {*model::materials[p->material]};
               for (auto i = 0; i < material.nuclide_.size(); ++i) {
                 auto j_nuclide = material.nuclide_[i];
                 auto atom_density = material.atom_density_(i);
@@ -986,7 +986,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           }
         } else {
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -1022,7 +1022,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
         } else {
           score = 0.;
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -1096,7 +1096,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
             * atom_density * flux;
         } else {
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -1144,7 +1144,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
         } else {
           score = 0.;
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -1186,7 +1186,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           }
         } else {
           if (p->material != MATERIAL_VOID) {
-            const Material& material {*model::materials[p->material-1]};
+            const Material& material {*model::materials[p->material]};
             for (auto i = 0; i < material.nuclide_.size(); ++i) {
               auto j_nuclide = material.nuclide_[i];
               auto atom_density = material.atom_density_(i);
@@ -1274,7 +1274,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
   // To significantly reduce de-referencing, point matxs to the macroscopic
   // Mgxs for the material of interest
-  data::macro_xs[p->material - 1].set_angle_index(p_uvw);
+  data::macro_xs[p->material].set_angle_index(p_uvw);
 
   // Do same for nucxs, point it to the microscopic nuclide data of interest
   if (i_nuclide >= 0) {
@@ -1937,7 +1937,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux,
   int filter_index)
 {
   const Tally& tally {*model::tallies[i_tally]};
-  const Material& material {*model::materials[p->material-1]};
+  const Material& material {*model::materials[p->material]};
 
   // Score all individual nuclide reaction rates.
   for (auto i = 0; i < material.nuclide_.size(); ++i) {
@@ -2051,12 +2051,12 @@ void score_analog_tally_mg(const Particle* p)
         double atom_density = 0.;
         if (i_nuclide >= 0) {
           //TODO: off-by-one
-          auto j = model::materials[p->material-1]
+          auto j = model::materials[p->material]
             ->mat_nuclide_index_[i_nuclide];
           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);
+          atom_density = model::materials[p->material]->atom_density_(j);
         }
 
         score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
@@ -2110,12 +2110,12 @@ score_tracklength_tally(const Particle* p, double distance)
           if (i_nuclide >= 0) {
             if (p->material != MATERIAL_VOID) {
               //TODO: off-by-one
-              auto j = model::materials[p->material-1]
+              auto j = model::materials[p->material]
                 ->mat_nuclide_index_[i_nuclide];
               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);
+              atom_density = model::materials[p->material]->atom_density_(j);
             }
           }
 
@@ -2180,12 +2180,12 @@ void score_collision_tally(const Particle* p)
           double atom_density = 0.;
           if (i_nuclide >= 0) {
             //TODO: off-by-one
-            auto j = model::materials[p->material-1]
+            auto j = model::materials[p->material]
               ->mat_nuclide_index_[i_nuclide];
             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);
+            atom_density = model::materials[p->material]->atom_density_(j);
           }
 
           //TODO: consider replacing this "if" with pointers or templates
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index e3a9bb479..0f18e245f 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -117,14 +117,11 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       // If this location is not in the geometry at all, move on to next block
       if (!find_cell(&p, false)) continue;
 
-      // TODO: off-by-one
-      int i_material = p.material == MATERIAL_VOID ? p.material : p.material - 1;
-
       if (domain_type_ == FILTER_MATERIAL) {
-        if (i_material != MATERIAL_VOID) {
+        if (p.material != MATERIAL_VOID) {
           for (int i_domain = 0; i_domain < n; i_domain++) {
-            if (model::materials[i_material]->id_ == domain_ids_[i_domain]) {
-              this->check_hit(i_material, indices[i_domain], hits[i_domain]);
+            if (model::materials[p.material]->id_ == domain_ids_[i_domain]) {
+              this->check_hit(p.material, indices[i_domain], hits[i_domain]);
               break;
             }
           }
@@ -133,7 +130,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         for (int level = 0; level < p.n_coord; ++level) {
           for (int i_domain=0; i_domain < n; i_domain++) {
             if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) {
-              this->check_hit(i_material, indices[i_domain], hits[i_domain]);
+              this->check_hit(p.material, indices[i_domain], hits[i_domain]);
               break;
             }
           }
@@ -142,7 +139,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         for (int level = 0; level < p.n_coord; ++level) {
           for (int i_domain = 0; i_domain < n; ++i_domain) {
             if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) {
-              check_hit(i_material, indices[i_domain], hits[i_domain]);
+              check_hit(p.material, indices[i_domain], hits[i_domain]);
               break;
             }
           }