diff --git a/CMakeLists.txt b/CMakeLists.txt
index dc9a1231b..6a2e5aee9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -397,6 +397,7 @@ add_library(libopenmc SHARED
   src/geometry_aux.cpp
   src/hdf5_interface.cpp
   src/lattice.cpp
+  src/material.cpp
   src/math_functions.cpp
   src/message_passing.cpp
   src/mgxs.cpp
diff --git a/src/api.F90 b/src/api.F90
index cc9120c7b..6a6c10671 100644
--- a/src/api.F90
+++ b/src/api.F90
@@ -214,7 +214,7 @@ contains
         if (p % material == MATERIAL_VOID) then
           id = 0
         else
-          id = materials(p % material) % id
+          id = materials(p % material) % id()
         end if
       end if
       instance = p % cell_instance - 1
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index 71c458105..5a2dca282 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -16,11 +16,11 @@ module geometry_header
   implicit none
 
   interface
-    function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr)
+    function cell_pointer(cell_ind) bind(C) result(ptr)
       import C_PTR, C_INT32_T
       integer(C_INT32_T), intent(in), value :: cell_ind
       type(C_PTR)                           :: ptr
-    end function cell_pointer_c
+    end function cell_pointer
 
     function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id)
       import C_PTR, C_INT32_T
@@ -110,12 +110,11 @@ module geometry_header
       integer(HID_T), intent(in), value :: group
     end subroutine cell_to_hdf5_c
 
-    function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') &
-         result(ptr)
+    function lattice_pointer(lat_ind) bind(C) result(ptr)
       import C_PTR, C_INT32_T
       integer(C_INT32_T), intent(in), value :: lat_ind
       type(C_PTR)                           :: ptr
-    end function lattice_pointer_c
+    end function lattice_pointer
 
     function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id)
       import C_PTR, C_INT32_T
@@ -571,7 +570,7 @@ contains
     ! Extend the C++ cells array and get pointers to the C++ objects
     call extend_cells_c(n)
     do i = n_cells - n, n_cells
-      cells(i) % ptr = cell_pointer_c(i - 1)
+      cells(i) % ptr = cell_pointer(i - 1)
     end do
 
     err = 0
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index da33a6627..aef9a14c0 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -87,6 +87,11 @@ module input_xml
       type(C_PTR) :: node_ptr
     end subroutine read_settings
 
+    subroutine read_materials(node_ptr) bind(C)
+      import C_PTR
+      type(C_PTR) :: node_ptr
+    end subroutine read_materials
+
     function find_root_universe() bind(C) result(root)
       import C_INT32_T
       integer(C_INT32_T) :: root
@@ -1051,7 +1056,7 @@ contains
     allocate(surfaces(n_surfaces))
 
     do i = 1, n_surfaces
-      surfaces(i) % ptr = surface_pointer_c(i - 1);
+      surfaces(i) % ptr = surface_pointer(i - 1);
 
       if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
 
@@ -1095,7 +1100,7 @@ contains
     do i = 1, n_cells
       c => cells(i)
 
-      c % ptr = cell_pointer_c(i - 1)
+      c % ptr = cell_pointer(i - 1)
 
       ! Initialize distribcell instances and distribcell index
       c % distribcell_index = NONE
@@ -1309,7 +1314,7 @@ contains
     RECT_LATTICES: do i = 1, n_rlats
       allocate(RectLattice::lattices(i) % obj)
       lat => lattices(i) % obj
-      lat % ptr = lattice_pointer_c(i - 1)
+      lat % ptr = lattice_pointer(i - 1)
       select type(lat)
       type is (RectLattice)
 
@@ -1325,7 +1330,7 @@ contains
     HEX_LATTICES: do i = 1, n_hlats
       allocate(HexLattice::lattices(n_rlats + i) % obj)
       lat => lattices(n_rlats + i) % obj
-      lat % ptr = lattice_pointer_c(n_rlats + i - 1)
+      lat % ptr = lattice_pointer(n_rlats + i - 1)
       select type (lat)
       type is (HexLattice)
 
@@ -1540,6 +1545,8 @@ contains
     call doc % load_file(filename)
     root = doc % document_element()
 
+    call read_materials(root % ptr)
+
     ! Get pointer to list of XML <material>
     call get_node_list(root, "material", node_mat_list)
 
@@ -1556,16 +1563,11 @@ contains
     do i = 1, n_materials
       mat => materials(i)
 
+      mat % ptr = material_pointer(i - 1)
+
       ! Get pointer to i-th material node
       node_mat = node_mat_list(i)
 
-      ! Copy material id
-      if (check_for_node(node_mat, "id")) then
-        call get_node_value(node_mat, "id", mat % id)
-      else
-        call fatal_error("Must specify id of material in materials XML file")
-      end if
-
       ! Check if material is depletable
       if (check_for_node(node_mat, "depletable")) then
         call get_node_value(node_mat, "depletable", temp_str)
@@ -1574,9 +1576,9 @@ contains
       end if
 
       ! Check to make sure 'id' hasn't been used
-      if (material_dict % has(mat % id)) then
+      if (material_dict % has(mat % id())) then
         call fatal_error("Two or more materials use the same unique ID: " &
-             // to_str(mat % id))
+             // to_str(mat % id()))
       end if
 
       ! Copy material name
@@ -1596,7 +1598,7 @@ contains
         node_dens = node_mat % child("density")
       else
         call fatal_error("Must specify density element in material " &
-             // trim(to_str(mat % id)))
+             // trim(to_str(mat % id())))
       end if
 
       ! Copy units
@@ -1626,7 +1628,7 @@ contains
         sum_density = .false.
         if (val <= ZERO) then
           call fatal_error("Need to specify a positive density on material " &
-               // trim(to_str(mat % id)) // ".")
+               // trim(to_str(mat % id())) // ".")
         end if
 
         ! Adjust material density based on specified units
@@ -1641,7 +1643,7 @@ contains
           mat % density = 1.0e-24_8 * val
         case default
           call fatal_error("Unkwown units '" // trim(units) &
-               // "' specified on material " // trim(to_str(mat % id)))
+               // "' specified on material " // trim(to_str(mat % id())))
         end select
       end if
 
@@ -1650,7 +1652,7 @@ contains
 
       if (size(node_ele_list) > 0) then
         call fatal_error("Unable to add an element to material " &
-             // trim(to_str(mat % id)) // " since the element option has &
+             // trim(to_str(mat % id())) // " since the element option has &
              &been removed from the xml input. Elements can only be added via &
              &the Python API, which will expand elements into their natural &
              &nuclides.")
@@ -1663,7 +1665,7 @@ contains
       if (.not. check_for_node(node_mat, "nuclide") .and. &
            .not. check_for_node(node_mat, "macroscopic")) then
         call fatal_error("No macroscopic data or nuclides specified on &
-             &material " // trim(to_str(mat % id)))
+             &material " // trim(to_str(mat % id())))
       end if
 
       ! Create list of macroscopic x/s based on those specified, just treat
@@ -1677,7 +1679,7 @@ contains
                          & mode!")
       else if (size(node_macro_list) > 1) then
         call fatal_error("Only one macroscopic object permitted per material, " &
-             // trim(to_str(mat % id)))
+             // trim(to_str(mat % id())))
       else if (size(node_macro_list) == 1) then
 
         node_nuc = node_macro_list(1)
@@ -1685,7 +1687,7 @@ contains
         ! Check for empty name on nuclide
         if (.not. check_for_node(node_nuc, "name")) then
           call fatal_error("No name specified on macroscopic data in material " &
-               // trim(to_str(mat % id)))
+               // trim(to_str(mat % id())))
         end if
 
         ! store nuclide name
@@ -1715,7 +1717,7 @@ contains
           ! Check for empty name on nuclide
           if (.not. check_for_node(node_nuc, "name")) then
             call fatal_error("No name specified on nuclide in material " &
-                 // trim(to_str(mat % id)))
+                 // trim(to_str(mat % id())))
           end if
 
           ! store nuclide name
@@ -1854,7 +1856,7 @@ contains
       if (.not. (all(mat % atom_density >= ZERO) .or. &
            all(mat % atom_density <= ZERO))) then
         call fatal_error("Cannot mix atom and weight percents in material " &
-             // to_str(mat % id))
+             // to_str(mat % id()))
       end if
 
       ! Determine density if it is a sum value
@@ -1935,7 +1937,7 @@ contains
       end if
 
       ! Add material to dictionary
-      call material_dict % set(mat % id, i)
+      call material_dict % set(mat % id(), i)
     end do
 
     ! Set total number of nuclides and S(a,b) tables
diff --git a/src/material.cpp b/src/material.cpp
new file mode 100644
index 000000000..49f0e2415
--- /dev/null
+++ b/src/material.cpp
@@ -0,0 +1,66 @@
+#include "material.h"
+
+#include "error.h"
+#include "xml_interface.h"
+
+//TODO: remove this include
+#include <iostream>
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+std::vector<Material*> global_materials;
+
+//==============================================================================
+// Material implementation
+//==============================================================================
+
+Material::Material(pugi::xml_node material_node)
+{
+  if (check_for_node(material_node, "id")) {
+    id = stoi(get_node_value(material_node, "id"));
+  } else {
+    fatal_error("Must specify id of material in materials XML file.");
+  }
+
+  std::cout << "Reading material with id " << id << std::endl;
+}
+
+//==============================================================================
+// Non-method functions
+//==============================================================================
+
+extern "C" void
+read_materials(pugi::xml_node* node)
+{
+  // Loop over XML material elements and populate the array.
+  for (pugi::xml_node material_node: node->children("material")) {
+    global_materials.push_back(new Material(material_node));
+  }
+}
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
+
+extern "C" {
+  Material* material_pointer(int32_t indx) {return global_materials[indx];}
+
+  int32_t material_id(Material* mat) {return mat->id;}
+
+  void material_set_id(Material* mat, int32_t id) {mat->id = id;}
+
+  void extend_materials_c(int32_t n)
+  {
+    global_materials.reserve(global_materials.size() + n);
+    for (int32_t i = 0; i < n; i++) {
+      global_materials.push_back(new Material());
+    }
+  }
+}
+
+} // namespace openmc
diff --git a/src/material.h b/src/material.h
new file mode 100644
index 000000000..f2e3df2fb
--- /dev/null
+++ b/src/material.h
@@ -0,0 +1,33 @@
+#ifndef OPENMC_CELL_H
+#define OPENMC_CELL_H
+
+#include <vector>
+
+#include "pugixml.hpp"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+class Material;
+extern std::vector<Material*> global_materials;
+
+//==============================================================================
+//! A substance with constituent nuclides and thermal scattering data
+//==============================================================================
+
+class Material
+{
+public:
+  int32_t id; //!< Unique ID
+
+  Material() {};
+
+  explicit Material(pugi::xml_node material_node);
+};
+
+} // namespace openmc
+#endif // OPENMC_CELL_H
diff --git a/src/material_header.F90 b/src/material_header.F90
index 03ecf69d9..5bfbadfcf 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -27,13 +27,39 @@ module material_header
   public :: openmc_material_set_density
   public :: openmc_material_set_densities
   public :: openmc_material_set_id
+  public :: material_pointer
+
+  interface
+    function material_pointer(mat_ind) bind(C) result(ptr)
+      import C_PTR, C_INT32_T
+      integer(C_INT32_T), intent(in), value :: mat_ind
+      type(C_PTR)                           :: ptr
+    end function material_pointer
+
+    function material_id_c(mat_ptr) bind(C, name='material_id') result(id)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: mat_ptr
+      integer(C_INT32_T)             :: id
+    end function material_id_c
+
+    subroutine material_set_id_c(mat_ptr, id) bind(C, name='material_set_id')
+      import C_PTR, C_INT32_T
+      type(C_PTR),        intent(in), value :: mat_ptr
+      integer(C_INT32_T), intent(in), value :: id
+    end subroutine material_set_id_c
+
+    subroutine extend_materials_c(n) bind(C)
+      import C_INT32_t
+      integer(C_INT32_T), intent(in), value :: n
+    end subroutine extend_materials_c
+  end interface
 
 !===============================================================================
 ! MATERIAL describes a material by its constituent nuclides
 !===============================================================================
 
   type, public :: Material
-    integer              :: id              ! unique identifier
+    type(C_PTR) :: ptr
     character(len=104)   :: name = ""       ! User-defined name
     integer              :: n_nuclides = 0  ! number of nuclides
     integer, allocatable :: nuclide(:)      ! index in nuclides array
@@ -67,6 +93,8 @@ module material_header
     logical, allocatable :: p0(:)
 
   contains
+    procedure :: id => material_id
+    procedure :: set_id => material_set_id
     procedure :: set_density => material_set_density
     procedure :: init_nuclide_index => material_init_nuclide_index
     procedure :: assign_sab_tables => material_assign_sab_tables
@@ -88,6 +116,18 @@ contains
 ! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm.
 !===============================================================================
 
+  function material_id(this) result(id)
+    class(Material), intent(in) :: this
+    integer(C_INT32_T)          :: id
+    id = material_id_c(this % ptr)
+  end function material_id
+
+  subroutine material_set_id(this, id)
+    class(Material),    intent(in) :: this
+    integer(C_INT32_T), intent(in) :: id
+    call material_set_id_c(this % ptr, id)
+  end subroutine material_set_id
+
   function material_set_density(this, density) result(err)
     class(Material), intent(inout) :: this
     real(8), intent(in) :: density
@@ -185,7 +225,7 @@ contains
       if (.not. found) then
         call fatal_error("S(a,b) table " // trim(this % &
              sab_names(k)) // " did not match any nuclide on material " &
-             // trim(to_str(this % id)))
+             // trim(to_str(this % id())))
       end if
     end do ASSIGN_SAB
 
@@ -195,7 +235,7 @@ contains
         if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then
           call fatal_error(trim( &
                nuclides(this % nuclide(i_sab_nuclides % data(j))) % name) &
-               // " in material " // trim(to_str(this % id)) // " was found &
+               // " in material " // trim(to_str(this % id())) // " was found &
                &in multiple S(a,b) tables. Each nuclide can only appear in &
                &one S(a,b) table per material.")
         end if
@@ -460,6 +500,7 @@ contains
     integer(C_INT32_T), optional, intent(out) :: index_end
     integer(C_INT) :: err
 
+    integer :: i
     type(Material), allocatable :: temp(:) ! temporary materials array
 
     if (n_materials == 0) then
@@ -481,6 +522,12 @@ contains
     if (present(index_end)) index_end = n_materials + n
     n_materials = n_materials + n
 
+    ! Extend the C++ materials array and get pointers to the C++ objects
+    call extend_materials_c(n)
+    do i = n_materials - n, n_materials
+      materials(i) % ptr = material_pointer(i - 1)
+    end do
+
     err = 0
   end function openmc_extend_materials
 
@@ -606,7 +653,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= size(materials)) then
-      id = materials(index) % id
+      id = materials(index) % id()
       err = 0
     else
       err = E_OUT_OF_BOUNDS
@@ -622,7 +669,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= n_materials) then
-      materials(index) % id = id
+      call materials(index) % set_id(id)
       call material_dict % set(id, index)
       err = 0
     else
diff --git a/src/plot.F90 b/src/plot.F90
index d8db97ec6..2fd898206 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -95,7 +95,7 @@ contains
             id = -1
           else
             rgb = pl % colors(p % material) % rgb
-            id = materials(p % material) % id
+            id = materials(p % material) % id()
           end if
         end associate
       else if (pl % color_by == PLOT_COLOR_CELLS) then
diff --git a/src/summary.F90 b/src/summary.F90
index 062b4efc2..507642f71 100644
--- a/src/summary.F90
+++ b/src/summary.F90
@@ -203,7 +203,7 @@ contains
             call write_dataset(cell_group, "material", MATERIAL_VOID)
           else
             call write_dataset(cell_group, "material", &
-                 materials(c % material(1)) % id)
+                 materials(c % material(1)) % id())
           end if
         else
           allocate(cell_materials(size(c % material)))
@@ -211,7 +211,7 @@ contains
             if (c % material(j) == MATERIAL_VOID) then
               cell_materials(j) = MATERIAL_VOID
             else
-              cell_materials(j) = materials(c % material(j)) % id
+              cell_materials(j) = materials(c % material(j)) % id()
             end if
           end do
           call write_dataset(cell_group, "material", cell_materials)
@@ -333,7 +333,7 @@ contains
     do i = 1, n_materials
       m => materials(i)
       material_group = create_group(materials_group, "material " // &
-           trim(to_str(m%id)))
+           trim(to_str(m%id())))
 
       if (m % depletable) then
         call write_attribute(material_group, "depletable", 1)
diff --git a/src/surface_header.F90 b/src/surface_header.F90
index 916f9c11c..4a5834ec2 100644
--- a/src/surface_header.F90
+++ b/src/surface_header.F90
@@ -8,13 +8,12 @@ module surface_header
   implicit none
 
   interface
-    pure function surface_pointer_c(surf_ind) &
-         bind(C, name='surface_pointer') result(ptr)
+    pure function surface_pointer(surf_ind) bind(C) result(ptr)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), intent(in), value :: surf_ind
       type(C_PTR)                       :: ptr
-    end function surface_pointer_c
+    end function surface_pointer
 
     pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id)
       use ISO_C_BINDING
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 2c660ce87..94abb42c0 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -3043,7 +3043,7 @@ contains
 
           case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
                 SCORE_NU_FISSION)
-            if (materials(p % material) % id == deriv % diff_material) then
+            if (materials(p % material) % id() == deriv % diff_material) then
               score = score * (flux_deriv + ONE &
                    / materials(p % material) % density_gpcc)
             else
@@ -3064,7 +3064,7 @@ contains
 
           case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
                 SCORE_NU_FISSION)
-            if (materials(p % material) % id == deriv % diff_material) then
+            if (materials(p % material) % id() == deriv % diff_material) then
               score = score * (flux_deriv + ONE &
                    / materials(p % material) % density_gpcc)
             else
@@ -3106,7 +3106,7 @@ contains
 
           case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
                 SCORE_NU_FISSION)
-            if (materials(p % material) % id == deriv % diff_material &
+            if (materials(p % material) % id() == deriv % diff_material &
                  .and. p % event_nuclide == deriv % diff_nuclide) then
               associate(mat => materials(p % material))
                 ! Search for the index of the perturbed nuclide.
@@ -3128,7 +3128,7 @@ contains
 
         case (ESTIMATOR_COLLISION)
           scoring_diff_nuclide = &
-               (materials(p % material) % id == deriv % diff_material) &
+               (materials(p % material) % id() == deriv % diff_material) &
                .and. (i_nuclide == deriv % diff_nuclide)
 
           select case (score_bin)
@@ -3138,7 +3138,7 @@ contains
 
           case (SCORE_TOTAL)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % total /= ZERO) then
               score = score * (flux_deriv &
                    + micro_xs(deriv % diff_nuclide) % total &
@@ -3152,7 +3152,7 @@ contains
 
           case (SCORE_SCATTER)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % total - material_xs % absorption /= ZERO) then
               score = score * (flux_deriv &
                    + (micro_xs(deriv % diff_nuclide) % total &
@@ -3168,7 +3168,7 @@ contains
 
           case (SCORE_ABSORPTION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % absorption /= ZERO) then
               score = score * (flux_deriv &
                    + micro_xs(deriv % diff_nuclide) % absorption &
@@ -3182,7 +3182,7 @@ contains
 
           case (SCORE_FISSION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % fission /= ZERO) then
               score = score * (flux_deriv &
                    + micro_xs(deriv % diff_nuclide) % fission &
@@ -3196,7 +3196,7 @@ contains
 
           case (SCORE_NU_FISSION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % nu_fission /= ZERO) then
               score = score * (flux_deriv &
                    + micro_xs(deriv % diff_nuclide) % nu_fission &
@@ -3242,7 +3242,7 @@ contains
             score = score * flux_deriv
 
           case (SCORE_TOTAL)
-            if (materials(p % material) % id == deriv % diff_material .and. &
+            if (materials(p % material) % id() == deriv % diff_material .and. &
                  micro_xs(p % event_nuclide) % total > ZERO) then
               associate(mat => materials(p % material))
                 ! Search for the index of the perturbed nuclide.
@@ -3267,7 +3267,7 @@ contains
             end if
 
           case (SCORE_SCATTER)
-            if (materials(p % material) % id == deriv % diff_material .and. &
+            if (materials(p % material) % id() == deriv % diff_material .and. &
                  (micro_xs(p % event_nuclide) % total &
                  - micro_xs(p % event_nuclide) % absorption) > ZERO) then
               associate(mat => materials(p % material))
@@ -3295,7 +3295,7 @@ contains
             end if
 
           case (SCORE_ABSORPTION)
-            if (materials(p % material) % id == deriv % diff_material .and. &
+            if (materials(p % material) % id() == deriv % diff_material .and. &
                  micro_xs(p % event_nuclide) % absorption > ZERO) then
               associate(mat => materials(p % material))
                 ! Search for the index of the perturbed nuclide.
@@ -3320,7 +3320,7 @@ contains
             end if
 
           case (SCORE_FISSION)
-            if (materials(p % material) % id == deriv % diff_material .and. &
+            if (materials(p % material) % id() == deriv % diff_material .and. &
                  micro_xs(p % event_nuclide) % fission > ZERO) then
               associate(mat => materials(p % material))
                 ! Search for the index of the perturbed nuclide.
@@ -3345,7 +3345,7 @@ contains
             end if
 
           case (SCORE_NU_FISSION)
-            if (materials(p % material) % id == deriv % diff_material .and. &
+            if (materials(p % material) % id() == deriv % diff_material .and. &
                  micro_xs(p % event_nuclide) % nu_fission > ZERO) then
               associate(mat => materials(p % material))
                 ! Search for the index of the perturbed nuclide.
@@ -3385,7 +3385,7 @@ contains
 
           case (SCORE_TOTAL)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % total > ZERO) then
               cum_dsig = ZERO
               associate(mat => materials(p % material))
@@ -3404,7 +3404,7 @@ contains
               end associate
               score = score * (flux_deriv &
                    + cum_dsig / material_xs % total)
-            else if (materials(p % material) % id == deriv % diff_material &
+            else if (materials(p % material) % id() == deriv % diff_material &
                  .and. material_xs % total > ZERO) then
               dsig_t = ZERO
               associate (nuc => nuclides(i_nuclide))
@@ -3423,7 +3423,7 @@ contains
 
           case (SCORE_SCATTER)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  (material_xs % total - material_xs % absorption) > ZERO) then
               cum_dsig = ZERO
               associate(mat => materials(p % material))
@@ -3444,7 +3444,7 @@ contains
               end associate
               score = score * (flux_deriv + cum_dsig &
                    / (material_xs % total - material_xs % absorption))
-            else if ( materials(p % material) % id == deriv % diff_material &
+            else if ( materials(p % material) % id() == deriv % diff_material &
                  .and. (material_xs % total - material_xs % absorption) > ZERO)&
                  then
               dsig_t = ZERO
@@ -3466,7 +3466,7 @@ contains
 
           case (SCORE_ABSORPTION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % absorption > ZERO) then
               cum_dsig = ZERO
               associate(mat => materials(p % material))
@@ -3485,7 +3485,7 @@ contains
               end associate
               score = score * (flux_deriv &
                    + cum_dsig / material_xs % absorption)
-            else if (materials(p % material) % id == deriv % diff_material &
+            else if (materials(p % material) % id() == deriv % diff_material &
                  .and. material_xs % absorption > ZERO) then
               dsig_a = ZERO
               associate (nuc => nuclides(i_nuclide))
@@ -3504,7 +3504,7 @@ contains
 
           case (SCORE_FISSION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % fission > ZERO) then
               cum_dsig = ZERO
               associate(mat => materials(p % material))
@@ -3523,7 +3523,7 @@ contains
               end associate
               score = score * (flux_deriv &
                    + cum_dsig / material_xs % fission)
-            else if (materials(p % material) % id == deriv % diff_material &
+            else if (materials(p % material) % id() == deriv % diff_material &
                  .and. material_xs % fission > ZERO) then
               dsig_f = ZERO
               associate (nuc => nuclides(i_nuclide))
@@ -3542,7 +3542,7 @@ contains
 
           case (SCORE_NU_FISSION)
             if (i_nuclide == -1 .and. &
-                 materials(p % material) % id == deriv % diff_material .and. &
+                 materials(p % material) % id() == deriv % diff_material .and. &
                  material_xs % nu_fission > ZERO) then
               cum_dsig = ZERO
               associate(mat => materials(p % material))
@@ -3563,7 +3563,7 @@ contains
               end associate
               score = score * (flux_deriv &
                    + cum_dsig / material_xs % nu_fission)
-            else if (materials(p % material) % id == deriv % diff_material &
+            else if (materials(p % material) % id() == deriv % diff_material &
                  .and. material_xs % nu_fission > ZERO) then
               dsig_f = ZERO
               associate (nuc => nuclides(i_nuclide))
@@ -3614,7 +3614,7 @@ contains
 
         case (DIFF_DENSITY)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material) then
+            if (mat % id() == deriv % diff_material) then
               ! phi is proportional to e^(-Sigma_tot * dist)
               ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
               ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
@@ -3625,7 +3625,7 @@ contains
 
         case (DIFF_NUCLIDE_DENSITY)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material) then
+            if (mat % id() == deriv % diff_material) then
               ! phi is proportional to e^(-Sigma_tot * dist)
               ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
               ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
@@ -3636,7 +3636,7 @@ contains
 
         case (DIFF_TEMPERATURE)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material) then
+            if (mat % id() == deriv % diff_material) then
               do l=1, mat % n_nuclides
                 associate (nuc => nuclides(mat % nuclide(l)))
                   if (nuc % mp_present .and. &
@@ -3690,7 +3690,7 @@ contains
 
         case (DIFF_DENSITY)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material) then
+            if (mat % id() == deriv % diff_material) then
               ! phi is proportional to Sigma_s
               ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
               ! (1 / phi) * (d_phi / d_rho) = 1 / rho
@@ -3701,7 +3701,7 @@ contains
 
         case (DIFF_NUCLIDE_DENSITY)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material &
+            if (mat % id() == deriv % diff_material &
                  .and. p % event_nuclide == deriv % diff_nuclide) then
               ! Find the index in this material for the diff_nuclide.
               do j = 1, mat % n_nuclides
@@ -3722,7 +3722,7 @@ contains
 
         case (DIFF_TEMPERATURE)
           associate (mat => materials(p % material))
-            if (mat % id == deriv % diff_material) then
+            if (mat % id() == deriv % diff_material) then
               do l=1, mat % n_nuclides
                 associate (nuc => nuclides(mat % nuclide(l)))
                   if (mat % nuclide(l) == p % event_nuclide .and. &
diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90
index 443dc6285..51ab26594 100644
--- a/src/tallies/tally_filter_material.F90
+++ b/src/tallies/tally_filter_material.F90
@@ -76,7 +76,7 @@ contains
 
     allocate(material_ids(size(this % materials)))
     do i = 1, size(this % materials)
-      material_ids(i) = materials(this % materials(i)) % id
+      material_ids(i) = materials(this % materials(i)) % id()
     end do
     call write_dataset(filter_group, "bins", material_ids)
   end subroutine to_statepoint_material
@@ -110,7 +110,7 @@ contains
     integer,               intent(in) :: bin
     character(MAX_LINE_LEN)           :: label
 
-    label = "Material " // to_str(materials(this % materials(bin)) % id)
+    label = "Material " // to_str(materials(this % materials(bin)) % id())
   end function text_label_material
 
 !===============================================================================
diff --git a/src/volume_calc.F90 b/src/volume_calc.F90
index 93087df89..e0c3e8859 100644
--- a/src/volume_calc.F90
+++ b/src/volume_calc.F90
@@ -196,7 +196,7 @@ contains
         i_material = p % material
         if (i_material /= MATERIAL_VOID) then
           do i_domain = 1, size(this % domain_id)
-            if (materials(i_material) % id == this % domain_id(i_domain)) then
+            if (materials(i_material) % id() == this % domain_id(i_domain)) then
               call check_hit(i_domain, i_material, indices, hits, n_mat)
             end if
           end do