diff --git a/BEAVRS/beavrs/builder.py b/BEAVRS/beavrs/builder.py
index be23b1a..d5386b3 100644
--- a/BEAVRS/beavrs/builder.py
+++ b/BEAVRS/beavrs/builder.py
@@ -61,6 +61,7 @@ class BEAVRS(object):
         self.dd_mesh_lower_left = None
         self.dd_mesh_upper_right = None
         self.dd_nodemap = None
+        self.settings_volumes = None
         self.dd_truncate = False
         self.dd_interactions = False
 
@@ -112,20 +113,30 @@ class BEAVRS(object):
     def write_openmc_geometry(self):
         self.openmc_geometry.export_to_xml()
 
+    @property
+    def materials(self):
+        return openmc.Materials(self.mats.values())
+
     def write_openmc_materials(self):
-        materials_file = openmc.Materials(self.mats.values())
-        materials_file.export_to_xml()
+        self.materials.export_to_xml()
+
+    @property
+    def plots(self):
+        plots = Plots(self.mats)
+        plot_file = openmc.Plots(plots.plots)
+        return plot_file
 
     def write_openmc_plots(self):
-        self.plots = Plots(self.mats)
-        plot_file = openmc.Plots(self.plots.plots)
-        plot_file.export_to_xml()
+        self.plots.export_to_xml()
 
-    def write_openmc_settings(self):
+    @property
+    def settings(self):
         settings_file = openmc.Settings()
         settings_file.batches = self.settings_batches
         settings_file.inactive = self.settings_inactive
         settings_file.particles = self.settings_particles
+        if self.settings_volumes:
+            settings_file.volume_calculations = self.settings_volumes
         if self.dd_mesh_dimension:
             settings_file.dd_mesh_dimension = self.dd_mesh_dimension
             settings_file.dd_mesh_lower_left = self.dd_mesh_lower_left
@@ -139,7 +150,10 @@ class BEAVRS(object):
         output = {'tallies': self.settings_output_tallies,
                   'summary': self.settings_summary}
         settings_file.output = output
-        settings_file.export_to_xml()
+        return settings_file
+
+    def write_openmc_settings(self):
+        self.settings.export_to_xml()
 
     def write_openmc_tallies(self):
         if len(self.tallies) == 0: return
@@ -150,6 +164,9 @@ class BEAVRS(object):
             tallies_file.append(tally)
         tallies_file.export_to_xml()
 
+    def set_volumes(self, samples=100000):
+        self.settings_volumes = openmc.VolumeCalculation(domains=[self.main_universe], samples=samples, lower_left=(-170, -170, 0), upper_right=(170, 170, 420))
+
     def set_params(self, particles=None, batches=None, inactive=None):
         if particles: self.settings_particles = particles
         if batches: self.settings_batches = batches
diff --git a/BEAVRS/beavrs/materials.py b/BEAVRS/beavrs/materials.py
index 445df48..73c40a1 100644
--- a/BEAVRS/beavrs/materials.py
+++ b/BEAVRS/beavrs/materials.py
@@ -145,6 +145,7 @@ def openmc_materials(ppm):
         mat_name = 'Fuel {0:1.1f}%'.format(name)
         mats[mat_name] = openmc.Material(name=mat_name)
         mats[mat_name].temperature = 300
+        mats[mat_name].volume = (np.pi*c.pelletOR**2)*(c.fuel_ActiveFuel_top-c.fuel_ActiveFuel_bot)
         mats[mat_name].set_density('g/cc', den)
         mats[mat_name].add_element('O', a_O, 'ao')
         mats[mat_name].add_element('U', a_U, 'ao', enrichment=enr*100)
diff --git a/BEAVRS/make_beavrs.py b/BEAVRS/make_beavrs.py
index 2fa5234..210c2b8 100755
--- a/BEAVRS/make_beavrs.py
+++ b/BEAVRS/make_beavrs.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 
-import os
+import openmc, os
+import openmc.deplete
 from beavrs.builder import BEAVRS
 from optparse import OptionParser
 
@@ -17,14 +18,40 @@ p.add_option('-d', '--2d', action='store_true', dest='is_2d',
 p.add_option('-s', '--symmetric', action='store_true', dest='is_symmetric',
              default=False, help='Create octant-symmetric input files,' \
              + ' not symmetric by default')
+p.add_option('-e', '--deplete', action='store_true', dest='run_deplete', default=False,
+             help='Run the depletion code for this setup')
 (options, args) = p.parse_args()
 
 if not len(args) == 0:
     p.print_help()
 
+chain_path = "/opt/xdata/endfb-vii.1-hdf5/chain_endfb71_pwr.xml"
+
 b = BEAVRS(is_symmetric=options.is_symmetric, is_2d=options.is_2d)
 b.write_openmc_geometry()
 b.write_openmc_materials()
 b.write_openmc_plots()
 b.write_openmc_settings()
 
+# openmc.calculate_volumes()
+
+###############################################################################
+#                      Depletion settings
+###############################################################################
+
+if options.run_deplete:
+    model = openmc.model.Model(geometry=b.openmc_geometry, materials=b.materials, settings=b.settings, tallies=b.tallies, plots=plots)
+
+    ## setting the transport operator
+    operator = openmc.deplete.CoupledOperator(model,chain_path,diff_burnable_mats=False,normalization_mode='fission-q',fission_q={"U235": 202.27e6})
+
+    ## setting the system linear power [W]
+    power = [45.50,91.01,136.51,182.02,227.52,273.03,324.45,324.45,273.03,273.03,227.52,182.02,136.51,91.01,45.50]
+    time_steps = [8640.00,8640.00,8640.00,8640.00,8640.00,8640.00,4268160.00,45792000.00,1728000.00,46915200.00,8640.00,8640.00,8640.00,8640.00,8640.00]
+
+    print(time_steps)
+
+    ## depleting usin a first-order predictor algorithm
+    integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power, timestep_units =  's')
+
+    integrator.integrate()