diff --git a/Depletion/pin_dep.py b/Depletion/pin_dep.py
new file mode 100644
index 0000000..950556b
--- /dev/null
+++ b/Depletion/pin_dep.py
@@ -0,0 +1,274 @@
+import openmc
+import openmc.deplete
+import numpy as np
+import matplotlib.pyplot as plt
+
+###############################################################################
+#                      Simulation Input File Parameters
+###############################################################################
+
+## OpenMC simulation parameters
+batches = 300
+inactive = 50
+particles = 100000
+
+## Depletion simulation parameters
+ts = 40*24*60*60 # s
+td = 10*24*60*60 # s
+
+time_steps = np.array([ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td,ts,ts,ts,ts,ts,ts,ts,ts,ts,ts,td])
+
+chain_file = '../chain_endfb71.xml'
+power = 17373 # W/cm, for 2D simulations only (use W for 3D)
+
+###############################################################################
+#                              Define materials
+###############################################################################
+
+## Instantiate some Materials and register the appropriate Nuclides
+zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
+zircaloy.set_density('g/cm3', 6.49012)
+zircaloy.temperature = 600.0
+zircaloy.add_nuclide('Sn112', 1.37041E-04, 'wo')
+zircaloy.add_nuclide('Sn114', 9.34709E-05, 'wo')
+zircaloy.add_nuclide('Sn115', 4.93221E-05, 'wo')
+zircaloy.add_nuclide('Sn116', 2.12757E-03, 'wo')
+zircaloy.add_nuclide('Sn117', 1.13348E-03, 'wo')
+zircaloy.add_nuclide('Sn118', 3.60520E-03, 'wo')
+zircaloy.add_nuclide('Sn119', 1.28950E-03, 'wo')
+zircaloy.add_nuclide('Sn120', 4.93338E-03, 'wo')
+zircaloy.add_nuclide('Sn122', 7.12568E-04, 'wo')
+zircaloy.add_nuclide('Sn124', 9.05713E-04, 'wo')
+zircaloy.add_nuclide('Fe54', 7.34249E-05, 'wo')
+zircaloy.add_nuclide('Fe56', 1.19419E-03, 'wo')
+zircaloy.add_nuclide('Fe57', 2.80871E-05, 'wo')
+zircaloy.add_nuclide('Fe58', 3.77460E-06, 'wo')
+zircaloy.add_nuclide('Cr50', 4.17870E-05, 'wo')
+zircaloy.add_nuclide('Cr52', 8.37035E-04, 'wo')
+zircaloy.add_nuclide('Cr53', 9.67293E-05, 'wo')
+zircaloy.add_nuclide('Cr54', 2.44824E-05, 'wo')
+zircaloy.add_nuclide('Ni58', 4.03133E-04, 'wo')
+zircaloy.add_nuclide('Ni60', 1.60602E-04, 'wo')
+zircaloy.add_nuclide('Ni61', 7.09947E-06, 'wo')
+zircaloy.add_nuclide('Ni62', 2.29762E-05, 'wo')
+zircaloy.add_nuclide('Ni64', 6.07653E-06, 'wo')
+zircaloy.add_nuclide('Zr90', 4.97990E-01, 'wo')
+zircaloy.add_nuclide('Zr91', 1.09810E-01, 'wo')
+zircaloy.add_nuclide('Zr92', 1.69691E-01, 'wo')
+zircaloy.add_nuclide('Zr94', 1.75712E-01, 'wo')
+zircaloy.add_nuclide('Zr96', 2.89111E-02, 'wo')
+
+
+uo2 = openmc.Material(name='Fuel Batch 1')
+uo2.set_density('g/cc' ,10.07)
+uo2.temperature = 900.0
+uo2.add_nuclide('O16'  ,1.18530E-01,'wo')
+uo2.add_nuclide('U234' ,3.92254E-04,'wo')
+uo2.add_nuclide('U235' ,4.36326E-02,'wo')
+uo2.add_nuclide('U238' ,8.37445E-01,'wo')
+
+uo2.depletable = True
+
+
+borated_water = openmc.Material(material_id=4, name='Borated water')
+borated_water.set_density('g/cm3', 0.7245)
+borated_water.temperature = 600.0
+borated_water.add_nuclide('B10',8.40178E-05,'wo')
+borated_water.add_nuclide('B11',3.71839E-04,'wo')
+borated_water.add_nuclide('H1', 1.11832E-01,'wo')
+borated_water.add_nuclide('H2', 3.35290E-05,'wo')
+borated_water.add_nuclide('O16',8.85414E-01,'wo' )
+borated_water.add_nuclide('O17',2.26386E-03,'wo' )
+borated_water.add_s_alpha_beta('c_H_in_H2O')
+
+
+
+###############################################################################
+#                             Create geometry
+###############################################################################
+
+## Instantiate ZCylinder surfaces
+fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.41266, name='Fuel OR')
+clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.474364, name='Clad IR')
+
+
+left = openmc.XPlane(surface_id=4, x0=-0.632458, name='left')
+right = openmc.XPlane(surface_id=5, x0=0.632458, name='right')
+bottom = openmc.YPlane(surface_id=6, y0=-0.632458, name='bottom')
+top = openmc.YPlane(surface_id=7, y0=0.632458, name='top')
+up = openmc.ZPlane(surface_id=8, z0=50.0, name='up')
+down = openmc.ZPlane(surface_id=9, z0=-50.0, name='down')
+
+left.boundary_type   = 'reflective'
+right.boundary_type  = 'reflective'
+top.boundary_type    = 'reflective'
+bottom.boundary_type = 'reflective'
+up.boundary_type     = 'reflective'
+down.boundary_type   = 'reflective'
+
+
+
+## Instantiate Cells
+fuel = openmc.Cell(cell_id=1, name='cell 1')
+clad = openmc.Cell(cell_id=3, name='cell 3')
+water = openmc.Cell(cell_id=4, name='cell 4')
+
+## Use surface half-spaces to define regions
+fuel.region = -fuel_or & -up & +down
+clad.region = +fuel_or & -clad_ir  & -up & +down
+water.region = +clad_ir & +left & -right & +bottom & -top & -up & +down
+
+## Register Materials with Cells
+fuel.fill = uo2
+clad.fill = zircaloy
+water.fill = borated_water
+
+## Instantiate Universe
+root = openmc.Universe(universe_id=0, name='root universe')
+
+## Register Cells with Universe
+root.add_cells([fuel, clad, water])
+
+## Instantiate a Geometry, register the root Universe
+geometry = openmc.Geometry(root)
+
+###############################################################################
+#                     Set volumes of depletable materials
+###############################################################################
+
+## Compute cell areas
+## area = {}
+## area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2
+
+## Set materials volume for depletion. Set to an area for 2D simulations
+uo2.volume = 53.49763956
+
+###############################################################################
+#                     Transport calculation settings
+###############################################################################
+
+## Instantiate a Settings object, set all runtime parameters, and export to XML
+settings_file = openmc.Settings()
+settings_file.batches = batches
+settings_file.inactive = inactive
+settings_file.particles = particles
+settings_file.photon_transport = True
+## Create an initial uniform spatial source distribution over fissionable zones
+bounds = [-0.632458, -0.632458, -50, 0.632458, 0.632458, 50]
+# uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
+uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+settings_file.source = openmc.source.IndependentSource(space=uniform_dist)
+## settings_file.electron_treatment = 'ttb'
+## settings_file.cutoff = {'energy_photon' : 1000.0}
+
+#entropy_mesh = openmc.RegularMesh()
+#entropy_mesh.lower_left = [-0.39218, -0.39218, -50]
+#entropy_mesh.upper_right = [0.39218, 0.39218, 50]
+#entropy_mesh.dimension = [10, 10, 1]
+#settings_file.entropy_mesh = entropy_mesh
+
+###############################################################################
+#                   Initialize and run depletion calculation
+###############################################################################
+tallies_file = openmc.Tallies()
+t1 = openmc.Tally(name='edep')
+t1.filters = [openmc.MaterialFilter([uo2,zircaloy,borated_water])]
+#t1.filters.append(openmc.ParticleFilter(['neutron']))
+#,'photon','electron','positron']))
+t1.scores = ['flux','fission', 'nu-fission','heating-local','heating']
+tallies_file.append(t1)
+tallies_file.export_to_xml()
+
+
+fission_q={"U232":     193.0442774*1E6,
+           "U233":     199.7961836*1E6,
+           "U234":     200.632851*1E6,
+           "U235":     202.27*1E6,
+           "U236":     203.4043119*1E6,
+           "U237":     196.429643*1E6,
+           "U238":     206.8513817*1E6,
+           "Np237":    205.3704803*1E6,
+           "Np238":    208.6993709*1E6,
+           "Pu238":    209.5400125*1E6,
+           "Pu239":    208.0185328*1E6,
+           "Pu240":    208.6125667*1E6,
+           "Pu241":    211.2377153*1E6,
+           "Pu242":    212.0721865*1E6,
+           "Am241":    211.2167986*1E6,
+           "Am242_m1": 215.1453706*1E6,
+           "Am243":    212.952779*1E6,
+           "Cm242":    212.7864913*1E6,
+           "Cm243":    213.375296*1E6,
+           "Cm244":    217.9267669*1E6,
+           "Cm245":    214.6240222*1E6,
+           "Cm246":    220.179494*1E6,
+           "Th232":    197.1083894*1E6,
+           "Pa231":    194.0999425*1E6}
+
+## op = openmc.deplete.Operator(geometry, settings_file, chain_file,fission_q=fission_q)
+op = openmc.deplete.Operator(geometry, settings_file, chain_file, energy_mode="energy-deposition")
+
+## Perform simulation using the predictor algorithm
+integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
+integrator.integrate()
+
+## openmc.run(mpi_args=['mpiexec', '-n', '32', '--bind-to', 'core'])
+
+###############################################################################
+#                    Read depletion calculation results
+###############################################################################
+
+## Open results file
+#results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
+
+## Obtain K_eff as a function of time
+#time, keff = results.get_eigenvalue()
+
+## Obtain U235 concentration as a function of time
+#time, n_U235 = results.get_atoms('1', 'U235')
+
+## Obtain Xe135 absorption as a function of time
+#time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
+# t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atom/cm3", time_units="d")
+
+###############################################################################
+#                            Generate plots
+###############################################################################
+
+#plt.figure()
+#plt.plot(time/(24*60*60), keff, label="K-effective")
+#plt.xlabel("Time (days)")
+#plt.ylabel("Keff")
+#plt.show()
+
+#plt.figure()
+#plt.plot(time/(24*60*60), n_U235, label="U 235")
+#plt.xlabel("Time (days)")
+#plt.ylabel("n U5 (-)")
+#plt.show()
+
+#plt.figure()
+#plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
+#plt.xlabel("Time (days)")
+#plt.ylabel("RR (-)")
+#plt.show()
+#plt.close('all')
+#tallies_file = openmc.Tallies()
+#energy_filter = openmc.EnergyFilter([0., 20.0e6])
+
+## Instantiate flux Tally in moderator and fuel
+#tally = openmc.Tally(name='flux')
+#tally.filters = [openmc.CellFilter(fuel)]
+#tally.filters.append(energy_filter)
+#tally.scores = ['flux']
+#tallies_file.append(tally)
+
+## Instantiate reaction rate Tally in moderator
+#tally = openmc.Tally(name='fuel rxn rates')
+#tally.filters = [openmc.CellFilter(fuel)]
+#tally.filters.append(energy_filter)
+#tally.scores = ['fission']
+#tally.nuclides = ['O16', 'H1']
+#tallies_file.append(tally)
+
+#tallies_file.export_to_xml()
diff --git a/MSRE/build-it.py b/MSRE/build-it.py
new file mode 100644
index 0000000..7a70941
--- /dev/null
+++ b/MSRE/build-it.py
@@ -0,0 +1,119 @@
+import openmc
+
+## Materials
+
+## fuel (partially enriched uranium):
+
+## 1200 F pg. 17 MSRE Design and Operations, part III, Nuclear Analysis
+temperature = 1200
+fuel = openmc.Material(name="MSRE Fuel")
+fuel.set_density("g/cm3", 2.146)
+fuel.temperature = (1200-32)*5/9
+fuel.add_nuclide('Li6', 0.0005, 'wo')
+fuel.add_nuclide('Li7', 10.9, 'wo')
+fuel.add_nuclide('F19', 66.8, 'wo')
+fuel.add_nuclide('Be9', 6.27, 'wo')
+fuel.add_element('Zr', 10.92, 'wo')
+fuel.add_nuclide('U235', 1.67, 'wo')
+fuel.add_nuclide('U238', 3.44, 'wo')
+
+## Moderator graphite
+## p. 87 msr operations (robertson) part i
+## 1200 F, pg. 17 MSRE Design and Operations, part iii, nuclear analysis
+moder = openmc.Material(name="Graphite Moderator")
+moder.set_density("g/cc", 1.86)
+moder.temperature = (1200-32)*5/9
+moder.add_nuclide('B10', 1.592e-5, 'wo')
+moder.add_nuclide('B11', 6.408e-5, 'wo')
+moder.add_element('V', 0.0009, 'wo')
+moder.add_element('S', 0.0005, 'wo')
+moder.add_element('C', 99.9852, 'wo')
+moder.add_s_alpha_beta(name='C_graphite')
+## ignoring the oxygen b/c low content and XS
+
+## Hastelloy Tank
+tank = openmc.Material(name="MSRE Tank")
+tank.set_density("g/cc", 8.86)
+tank.temperature = (1200-32)*5/9
+tank.add_nuclide('Si30', 0.00030872, 'wo')
+tank.add_nuclide('W186', 0.0014215, 'wo')
+tank.add_nuclide('Mn55', 0.008, 'wo')
+tank.add_nuclide('W184', 0.001532, 'wo')
+tank.add_nuclide('W183', 0.0007155, 'wo')
+tank.add_nuclide('W182', 0.001325, 'wo')
+tank.add_nuclide('Ni58', 0.4721201092, 'wo')
+tank.add_nuclide('Mo92', 0.023744, 'wo')
+tank.add_nuclide('Fe58', 0.000141, 'wo')
+tank.add_nuclide('Mo94', 0.0148, 'wo')
+tank.add_nuclide('Mo95', 0.025472, 'wo')
+tank.add_nuclide('Mo96', 0.026688, 'wo')
+tank.add_nuclide('Mo97', 0.01528, 'wo')
+tank.add_nuclide('Mo98', 0.038608, 'wo')
+tank.add_nuclide('Fe56', 0.045877, 'wo')
+tank.add_nuclide('Fe57', 0.0010595, 'wo')
+tank.add_nuclide('Fe54', 0.0029225, 'wo')
+tank.add_nuclide('Mo100', 0.015408, 'wo')
+tank.add_nuclide('Si28', 0.00922297, 'wo')
+tank.add_nuclide('Si29', 0.00046832, 'wo')
+tank.add_nuclide('Cr50', 0.0030415, 'wo')
+tank.add_nuclide('Cr52', 0.0586523, 'wo')
+tank.add_nuclide('Cr53', 0.0066507, 'wo')
+tank.add_nuclide('Cr54', 0.0016555, 'wo')
+tank.add_nuclide('Cu65', 0.00107905, 'wo')
+tank.add_nuclide('Cu63', 0.00242095, 'wo')
+tank.add_nuclide('Ni64', 0.0064191286, 'wo')
+tank.add_nuclide('Ni61', 0.0079053205, 'wo')
+tank.add_nuclide('Ni60', 0.1818598208, 'wo')
+tank.add_nuclide('Ni62', 0.025205621, 'wo')
+
+## Control rod absorber material
+## operations report i; p. 102
+## 70wtpct gad III oxide, 30wtpct al III oxide
+ctrlPois = openmc.Material(name='Control Rod Absorbers')
+ctrlPois.set_density("g/cc", 5.873)
+ctrlPois.temperature = (1200-32)*5/9
+ctrlPois.add_element('Gd', 0.003862068965517241, 'ao')
+ctrlPois.add_nuclide('Al27', 0.005884660651235779, 'ao')
+ctrlPois.add_nuclide('O16', 0.014620094425129529, 'ao')
+
+air = openmc.Material(name='Air')
+air.set_density("g/cc", 0.001225)
+air.temperature = (1200-32)*5/9
+air.add_nuclide('O16', 0.21, 'ao')
+air.add_nuclide('N14', 0.79, 'ao')
+
+surf11 = openmc.Surface(id='11',bound=inf)
+surf12 = openmc.Surface(id='12',bound=inf)
+surf13 = openmc.Surface(id='13',bound=inf)
+surf14 = openmc.Surface(id='14',bound=inf)
+
+cell1 = openmc.Cell(id='111', material=fuel, bound=-11)
+cell1 = openmc.Cell(id='112', material=moder, bound=-12)
+cell1 = openmc.Cell(id='113', material=tank, bound=-13)
+cell1 = openmc.Cell(id='114', material=ctrlPois, bound=-14)
+cell1 = openmc.Cell(id='115', material=air, bound=-15)
+
+##### Single Unit Cell
+
+## 4 Bounding planes of stringer
+surf1 = openmc.XPlane(x0=-2.54)
+surf2 = openmc.XPlane(x0=2.54)
+surf3 = openmc.YPlane(y0=-2.54)
+surf4 = openmc.YPlane(y0=2.54)
+
+## 4 square cylinders form slot edges
+surf5 = openmc.ZCylinder()
+
+
+
+
+
+
+
+
+## MSR Cluster
+msr = openmc.Universe(name='MSR2G-partially-enriched-U-full-core')
+
+
+
+
diff --git a/Rep_Na_3_55.py b/Rep_Na_3_55.py
new file mode 100644
index 0000000..13174a1
--- /dev/null
+++ b/Rep_Na_3_55.py
@@ -0,0 +1,306 @@
+## BISON depletion capability comparison
+## Rep_Na_3 assessment case, fuel pellet mesh with 55 radial divisions.
+## BISON solves the depletion for regions with the same radius but different volumes.
+## The number of radial divisions in the mesh controls the number of the depletion regions
+
+try:
+    from mpi4py import MPI
+    rank = MPI.COMM_WORLD.Get_rank()
+except ModuleNotFoundError:
+    rank = 0
+
+import math, sys
+import openmc
+import openmc.deplete
+import numpy as np
+from matplotlib import pyplot
+
+######################################################################
+#                      Operational Parameters
+######################################################################
+
+fuel_temp = 978.31
+gap_temp = 663.69
+water_den = 0.7051
+ppm_boron = 500
+water_temp = 585.25
+
+######################################################################
+#                      Geometric Parameters
+######################################################################
+
+Rfo = 0.40956		## pellet outer radius
+Rci = 0.41776		## clad inner radius
+Rco = 0.47736		## clad outer radius
+hpitch = 0.63		## half-pitch between rods
+fheight = 1.3589	## fuel height
+nr = 55				## number of radial divisions
+
+######################################################################
+#                      Chain Parameters
+######################################################################
+
+#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_simple.xml"
+#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_casl_pwr.xml"
+#chain_path = "/storage/home/lca5209/OpenMC/chains/chain_endfb71_pwr.xml"
+chain_path = "/opt/xdata/endfb-vii.1-hdf5/chain_endfb71_pwr.xml"
+
+######################################################################
+#                      Settings
+######################################################################
+
+batches = 100
+inactive = 20
+particles = 50000
+seed = np.random.randint(2147483647)*2+1
+
+######################################################################
+#                      Setting the materials
+######################################################################
+## fuel = 3.0 wt% enriched uranium dioxide
+
+fuel = openmc.Material(material_id=1,name="UO2")
+
+fuel.add_nuclide("U238",0.84223593,'wo')
+fuel.add_nuclide("U235",0.03918626,'wo')
+fuel.add_nuclide("O16",0.11857781,'wo')
+fuel.set_density("g/cc",10.394075)
+fuel.temperature = fuel_temp
+fuel.volume = (math.pi*Rfo**2)*fheight
+
+mat_list = []
+
+for i in range(nr):
+    mat = fuel.clone()
+    mat.volume = (math.pi*((((i+1)/55)*Rfo)**2-(((i)/55)*Rfo)**2))*fheight
+    mat_list.append(mat)
+
+## gap = helium
+
+gap = openmc.Material(name='He')
+gap.set_density('g/cc',0.001598)
+gap.add_element('He', 1.0,'wo')
+gap.temperature = gap_temp
+mat_list.append(gap)
+
+## clad = zircaloy
+
+clad = openmc.Material(name='Zircaloy_4')
+clad.set_density('g/cc', 6.56)
+clad.add_element('Sn', 0.014  , 'wo')
+clad.add_element('Fe', 0.00165, 'wo')
+clad.add_element('Cr', 0.001  , 'wo')
+clad.add_element('Zr', 0.98335, 'wo')
+mat_list.append(clad)
+
+## coolant = borated water
+
+coolant = openmc.Material(name='Borated_water')
+coolant.set_density('g/cc', water_den)
+coolant.add_nuclide('B10',ppm_boron*0.199E-06,'wo')
+coolant.add_nuclide('B11',ppm_boron*0.801E-06,'wo')
+coolant.add_element('H',0.1108612,'wo')
+coolant.add_element('O',0.8886388,'wo')
+coolant.remove_nuclide('O18')
+coolant.temperature = water_temp
+coolant.add_s_alpha_beta('c_H_in_H2O')
+mat_list.append(coolant)
+
+sys.stdout.flush()
+try:
+    MPI.COMM_WORLD.Barrier()
+except NameError:
+    pass
+if rank == 0:
+    print(mat_list)
+sys.stdout.flush()
+
+materials = openmc.Materials(mat_list)
+materials.cross_sections="/opt/xdata/endfb-vii.1-hdf5/cross_sections.xml"
+materials.export_to_xml()
+
+######################################################################
+#                      Building the geometry
+######################################################################
+
+## Fuel rings, gap, and cladding
+
+cylinder = []
+
+counter_1 = 0
+for i in range(nr):
+    c = openmc.ZCylinder(surface_id=(i+1), x0=0, y0=0, r=((i+1)/55)*Rfo, name='Fuel '+ str(i+1) + ' OR')
+    cylinder.append(c)
+    counter_1 +=1
+cylinder_g = openmc.ZCylinder(surface_id=counter_1+1, x0=0, y0=0, r=Rci, name='Gap OR')
+cylinder_c = openmc.ZCylinder(surface_id=counter_1+2, x0=0, y0=0, r=Rco, name='Clad OR')
+cylinder.append(cylinder_g)
+cylinder.append(cylinder_c)
+
+## Box surrounding fuel rod
+left = openmc.XPlane(surface_id=counter_1+3, x0=-hpitch, name='left')
+right = openmc.XPlane(surface_id=counter_1+4, x0=hpitch, name='right')
+bottom = openmc.YPlane(surface_id=counter_1+5, y0=-hpitch, name='bottom')
+top = openmc.YPlane(surface_id=counter_1+6, y0=hpitch, name='top')
+
+## Surfaces above/below fuel
+fsouth = openmc.ZPlane(surface_id=counter_1+7, z0=-fheight/2.0, name='fsouth')
+fnorth = openmc.ZPlane(surface_id=counter_1+8, z0=fheight/2.0, name='fnorth')
+
+## Boundary Conditions
+left.boundary_type = 'reflective'
+right.boundary_type = 'reflective'
+top.boundary_type = 'reflective'
+bottom.boundary_type = 'reflective'
+fsouth.boundary_type = 'reflective'
+fnorth.boundary_type = 'reflective'
+
+cell = []
+
+counter_2 = 0
+for i in range(nr):
+    f = openmc.Cell(cell_id=(i+1), name='fuel '+str(i+1))
+    cell.append(f)
+    counter_2 += 1
+gap_cell = openmc.Cell(cell_id=counter_2+1, name='Gap cell')
+clad_cell = openmc.Cell(cell_id=counter_2+2, name='Clad cell')
+water_cell = openmc.Cell(cell_id=counter_2+3, name='Water cell')
+cell.append(gap_cell)
+cell.append(clad_cell)
+cell.append(water_cell)
+
+for i in range(len(cell)-1):
+    if i==0:
+        cell[i].region = -cylinder[i] & -fnorth & +fsouth
+    elif i==(len(cell)-2):
+        cell[i+1].region = +cylinder[i] & +left & -right & +bottom & -top & -fnorth & +fsouth
+        cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
+    else:
+        cell[i].region = +cylinder[(i-1)] & -cylinder[i] & -fnorth & +fsouth
+
+count_3 = 0
+for i in range(nr):
+    cell[i].fill = mat_list[i]
+    count_3 +=1
+cell[count_3].fill = gap
+cell[count_3+1].fill = clad
+cell[count_3+2].fill = coolant
+
+## Instantiate Universe
+root = openmc.Universe(universe_id=0, name='root universe')
+
+## Register Cells with Universe
+root.add_cells(cell)
+
+## Instantiate a Geometry, register the root Universe, and export to XML
+geometry = openmc.Geometry(root)
+geometry.export_to_xml()
+
+plot = openmc.Plot.from_geometry(geometry)
+plot.pixels = (250, 250)
+plot.filename = 'pinplot_55_regions'
+plot.width = (hpitch*2, hpitch*2)
+plot.pixels = (500, 500)
+plot.color_by = 'cell'
+plot.colors = {cell[-3]: 'purple',cell[-2]: 'grey',
+               cell[-1]: 'blue'}
+plots = openmc.Plots([plot])
+plots.export_to_xml()
+openmc.plot_geometry()
+
+###############################################################################
+#                   Exporting to OpenMC settings.xml file
+###############################################################################
+
+## Instantiate a Settings object, set all runtime parameters, and export to XML
+settings = openmc.Settings()
+settings.batches = batches
+settings.inactive = inactive
+settings.particles = particles
+settings.seed = seed
+settings.verbosity = 7
+settings.photon_transport = False
+
+## Cross-section temperature information
+settings.temperature = {
+    'method': 'interpolation',
+    'tolerance': 100 ,
+    'multipole': True
+}
+
+## Create an initial uniform spatial source distribution over fissionable zones
+bounds = [-hpitch, -hpitch, -fheight/2, hpitch, hpitch, fheight/2]
+uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
+settings.source = openmc.source.IndependentSource(space=uniform_dist)
+
+entropy_mesh = openmc.RegularMesh()
+entropy_mesh.lower_left = [-hpitch, -hpitch, -1.e50]
+entropy_mesh.upper_right = [hpitch, hpitch, 1.e50]
+entropy_mesh.dimension = [10, 10, 1]
+
+settings.entropy_mesh = entropy_mesh
+settings.export_to_xml()
+
+###############################################################################
+#                   Exporting to OpenMC tallies.xml file
+###############################################################################
+
+## Instantiate a tally mesh
+mesh = openmc.RegularMesh()
+mesh.type = 'regular'
+mesh.dimension = [1, 1, 100]
+mesh.lower_left = [-hpitch, -hpitch, -fheight/2]
+mesh.upper_right = [hpitch, hpitch, fheight/2]
+
+## Instantiate some tally Filters
+energy_filter = openmc.EnergyFilter([0., 0.625, 20.e6])
+energy_filter2 = openmc.EnergyFilter(np.logspace(-3, 7, num=100))
+mesh_filter = openmc.MeshFilter(mesh)
+
+## Instantiate the Tally
+tally = openmc.Tally(tally_id=1, name='Spatial Tally 1')
+tally.filters = [energy_filter, mesh_filter]
+tally.scores = ['flux']#, 'fission', 'nu-fission']
+
+## Instantiate the Tally
+tallyspec = openmc.Tally(tally_id=2, name='Energy Spectrum Tally 1')
+tallyspec.filters = [energy_filter2]
+tallyspec.scores = ['flux']#, 'fission', 'nu-fission']
+
+particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
+mat_filter = openmc.MaterialFilter(materials)
+
+tallydep = openmc.Tally(tally_id=3, name='Spatial Tally 2')
+tallydep.filters = [mat_filter, particle_filter]
+tallydep.scores = ['fission', 'heating', 'heating-local']
+
+## Instantiate a Tallies collection and export to XML
+tallies = openmc.Tallies([tally, tallyspec,tallydep])
+tallies.export_to_xml()
+
+###############################################################################
+#                      Depletion settings
+###############################################################################
+
+model = openmc.model.Model(geometry=geometry, materials=materials, settings=settings, tallies=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]
+
+sys.stdout.flush()
+try:
+    MPI.COMM_WORLD.Barrier()
+except NameError:
+    pass
+if rank == 0:
+    print(time_steps)
+sys.stdout.flush()
+
+## depleting usin a first-order predictor algorithm
+integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power, timestep_units =  's')
+
+integrator.integrate()
diff --git a/mc-performance/my_attempt.py b/mc-performance/my_attempt.py
new file mode 100644
index 0000000..793a99f
--- /dev/null
+++ b/mc-performance/my_attempt.py
@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+
+import openmc
+
+## Monte Carlo Performance Benchmark, version 1.2, July 2011
+
+fuel = openmc.Material(name="Fuel", material_id=1)
+fuel.set_density('g/cm3', 10.062)
+fuel.add_nuclide("U234", 4.9476e-6)
+fuel.add_nuclide("U235", 4.8218e-4)
+fuel.add_nuclide("U236", 9.0402e-5)
+fuel.add_nuclide("U238", 2.1504e-2)
+fuel.add_nuclide("Np237", 7.3733e-6)
+fuel.add_nuclide("Pu238", 1.5148e-6)
+fuel.add_nuclide("Pu239", 1.3955e-4)
+fuel.add_nuclide("Pu240", 3.4405e-5)
+fuel.add_nuclide("Pu241", 2.1439e-5)
+fuel.add_nuclide("Pu242", 3.7422e-6)
+fuel.add_nuclide("Am241", 4.5041e-7)
+fuel.add_nuclide("Am242_m1", 9.2301e-9)
+fuel.add_nuclide("Am243", 4.7878e-7)
+fuel.add_nuclide("Cm242", 1.0485e-7)
+fuel.add_nuclide("Cm243", 1.4268e-9)
+fuel.add_nuclide("Cm244", 8.8756e-8)
+fuel.add_nuclide("Cm245", 3.5285e-9)
+fuel.add_nuclide("Mo95", 2.6497e-5)
+fuel.add_nuclide("Tc99", 3.2772e-5)
+fuel.add_nuclide("Ru101", 3.0742e-5)
+fuel.add_nuclide("Ru103", 2.3505e-6)
+fuel.add_nuclide("Ag109", 2.0009e-6)
+fuel.add_nuclide("Xe135", 1.0801e-8)
+fuel.add_nuclide("Cs133", 3.4612e-5)
+fuel.add_nuclide("Nd143", 2.6078e-5)
+fuel.add_nuclide("Nd145", 1.9898e-5)
+fuel.add_nuclide("Sm147", 1.6128e-6)
+fuel.add_nuclide("Sm149", 1.1627e-7)
+fuel.add_nuclide("Sm150", 7.1727e-6)
+fuel.add_nuclide("Sm151", 5.4947e-7)
+fuel.add_nuclide("Sm152", 3.0221e-6)
+fuel.add_nuclide("Eu153", 2.6209e-6)
+fuel.add_nuclide("Gd155", 1.5369e-9)
+fuel.add_nuclide("O16", 4.5737e-2)
+
+clad = openmc.Material(name='Cladding', material_id=2)
+clad.set_density('g/cm3', 5.77)
+clad.add_element("Zr", 1.0)
+
+cold_water = openmc.Material(name='Cold borated water', material_id=3)
+cold_water.set_density("atom/b-cm", 0.07416)
+cold_water.add_nuclide("H1", 2.0)
+cold_water.add_nuclide("O16", 1.0)
+cold_water.add_nuclide("B10", 6.490e-4)
+cold_water.add_nuclide("B11", 2.689e-3)
+cold_water.add_s_alpha_beta('c_H_in_H2O')
+
+hot_water = openmc.Material(name='Hot borated water', material_id=4)
+hot_water.set_density("atom/b-cm", 0.06614)
+hot_water.add_nuclide("H1", 2.0)
+hot_water.add_nuclide("O16", 1.0)
+hot_water.add_nuclide("B10", 6.490e-4)
+hot_water.add_nuclide("B11", 2.689e-3)
+hot_water.add_s_alpha_beta('c_H_in_H2O')
+
+rpv_steel = openmc.Material(name='Reactor pressure vessel steel',
+                            material_id=5)
+rpv_steel.set_density("g/cm3", 7.9)
+
+
+
+
+
+
+lower_rad_ref = openmc.Material(name='Lower radial reflector',
+                                material_id=6)
+lower_rad_ref.set_density("g/cm3", 4.32)
+
+
+
+
+
+
+upper_rad_ref = openmc.Material(name='Upper radial reflector /'
+                                'Top plate region', material_id=7)
+upper_rad_ref.set_density("g/cm3", 4.28)
+
+
+
+
+bot_plate = openmc.Material(name='Bottom plate region', material_id=8)
+bot_plate.set_density("g/cm3", 7.184)
+
+
+
+
+
+
+bot_nozzle = openmc.Material(name='Bottom nozzle region', material_id=9)
+bot_nozzle.set_density("g/cm3", 2.53)
+
+
+
+top_nozzle = openmc.Material(name='Top nozzle region', material_id=10)
+top_nozzle.set_density("g/cm3", 1.746)
+
+
+
+
+
+
+top_fa = openmc.Material(name='Top of fuel assemblies', material_id=11)
+top_fa.set_density("g/cm3", 3.044)
+
+
+
+
+
+bot_fa = openmc.Material(name='Bottom of fuel assemblies', material_id=12)
+bot_fa.set_density("g/cm3", 1.762)
+
+
+
+
+# Define the materials file.
+materials = openmc.Materials()
+materials.default_xs = '71c'
+materials += [fuel, , , , , , , , , , , , bot_fa]
+materials.export_to_xml()
+
+
+# Define surfaces.
\ No newline at end of file
diff --git a/own_make/build_AP1000.py b/own_make/build_AP1000.py
new file mode 100644
index 0000000..4e2f0fb
--- /dev/null
+++ b/own_make/build_AP1000.py
@@ -0,0 +1,401 @@
+import openmc
+from openmc.data import atomic_mass, atomic_weight
+
+MB10 = atomic_mass('B10')
+MB11 = atomic_mass('B11')
+MU234 = atomic_mass('U234')
+MU235 = atomic_mass('U235')
+MU238 = atomic_mass('U238')
+
+materials = openmc.Materials()
+air = openmc.Material(name='Air')
+air.temperature = 300
+air.set_density('g/cc', 0.000616)
+air.add_element('O', 0.2095, 'ao')
+air.add_element('N', 0.7809, 'ao')
+air.add_element('Ar', 0.00933, 'ao')
+air.add_element('C', 0.00027, 'ao')
+
+## Fuels
+'''
+fuel_16 = openmc.Material(name='UO2 (1.6%)')
+fuel_16.set_density('g/cm3', 10.31341)
+fuel_16.temperature = 300
+fuel_16.add_nuclide('U234', 3.0131e-06)
+fuel_16.add_nuclide('U235', 3.7503e-04)
+fuel_16.add_nuclide('U238', 2.2625e-02)
+fuel_16.add_nuclide('O16', 4.5897e-02)
+fuel_16.add_nuclide('O17', 1.7436e-05)
+fuel_16.add_nuclide('O18', 9.2032e-05)
+
+fuel_24 = openmc.Material(name='UO2 (2.4%)')
+fuel_24.set_density('g/cm3', 10.29748)
+fuel_24.temperature = 300
+fuel_24.add_nuclide('U234', 4.4842e-06)
+fuel_24.add_nuclide('U235', 5.5814e-04)
+fuel_24.add_nuclide('U238', 2.2407e-02)
+fuel_24.add_nuclide('O16', 4.5830e-02)
+fuel_24.add_nuclide('O17', 1.7411e-05)
+fuel_24.add_nuclide('O18', 9.1898e-05)
+
+fuel_31 = openmc.Material(name='UO2 (3.1%)')
+fuel_31.set_density('g/cm3', 10.30166)
+fuel_31.temperature = 300
+fuel_31.add_nuclide('U234', 5.7987e-06)
+fuel_31.add_nuclide('U235', 7.2175e-04)
+fuel_31.add_nuclide('U238', 2.2253e-02)
+fuel_31.add_nuclide('O16', 4.5853e-02)
+fuel_31.add_nuclide('O17', 1.7420e-05)
+fuel_31.add_nuclide('O18', 9.1942e-05)
+
+fuel_32 = openmc.Material(name='UO2 (3.2%)')
+fuel_32.set_density('g/cm3', 10.34115)
+fuel_32.temperature = 300
+fuel_32.add_nuclide('U234', 5.9959e-06)
+fuel_32.add_nuclide('U235', 7.4630e-04)
+fuel_32.add_nuclide('U238', 2.2317e-02)
+fuel_32.add_nuclide('O16', 4.6029e-02)
+fuel_32.add_nuclide('O17', 1.7487e-05)
+fuel_32.add_nuclide('O18', 9.2296e-05)
+
+fuel_34 = openmc.Material(name='UO2 (3.4%)')
+fuel_34.set_density('g/cm3', 10.35917)
+fuel_34.temperature = 300
+fuel_34.add_nuclide('U234', 6.4018e-06)
+fuel_34.add_nuclide('U235', 7.9681e-04)
+fuel_34.add_nuclide('U238', 2.2307e-02)
+fuel_34.add_nuclide('O16', 4.6110e-02)
+fuel_34.add_nuclide('O17', 1.7517e-05)
+fuel_34.add_nuclide('O18', 9.2459e-05)
+'''
+
+## Specify enrichments to be calculated
+enrs = [0.0161006, 0.0239993, 0.0310221, 0.0319547, 0.0340585]
+dens = [10.31341, 10.29748, 10.30166, 10.34115, 10.35917]
+
+# Loop around enrichments
+for enr, den in zip(enrs, dens):
+
+    # Calculate molar mass of Uranium
+    enr_25 = enr
+    enr_24 = 0.008*enr_25
+    enr_28 = 1.0 - (enr_24 + enr_25)
+    MU = 1.0/(enr_24/MU234 + enr_25/MU235 + enr_28/MU238)
+
+    # Determine molar mass of UO2
+    MUO2 = MU + 2.0*atomic_weight('O')
+
+    # Compute weight percent of U in UO2
+    wUpUO2 = MU/MUO2
+
+    # Calculate Uranium atom fraction
+    a_U = wUpUO2*MUO2/MU
+
+    # Calculate Oxygen atom fraction
+    a_O = (1.0 - wUpUO2)*MUO2/atomic_weight('O')
+
+    # Create material
+    name = enr*100
+    mat_name = 'Fuel {0:1.1f}%'.format(name)
+    mat = openmc.Material(name=mat_name)
+    mat.temperature = 300
+    mat.set_density('g/cc', den)
+    mat.add_element('O', a_O, 'ao')
+    mat.add_element('U', a_U, 'ao', enrichment=enr*100)
+    materials += mat
+
+
+borosilicate_glass = openmc.Material(name='Borosilicate Glass')
+borosilicate_glass.set_density('g/cm3', 2.26)
+borosilicate_glass.temperature = 300
+borosilicate_glass.add_nuclide('Al27', 1.7352e-03)
+borosilicate_glass.add_nuclide('B10', 9.6506e-04)
+borosilicate_glass.add_nuclide('B11', 3.9189e-03)
+borosilicate_glass.add_nuclide('O16', 4.6514e-02)
+borosilicate_glass.add_nuclide('O17', 1.7671e-05)
+borosilicate_glass.add_nuclide('O18', 9.3268e-05)
+borosilicate_glass.add_nuclide('Si28', 1.6926e-02)
+borosilicate_glass.add_nuclide('Si29', 8.5944e-04)
+borosilicate_glass.add_nuclide('Si30', 5.6654e-04)
+
+ag_in_cd_rods = openmc.Material(name='Ag-In-Cd Control Rod')
+ag_in_cd_rods.set_density('g/cm3', 10.16)
+ag_in_cd_rods.temperature = 300
+ag_in_cd_rods.add_nuclide('Ag107', 2.3523e-02)
+ag_in_cd_rods.add_nuclide('Ag109', 2.1854e-02)
+ag_in_cd_rods.add_nuclide('Cd106', 3.3882e-05)
+ag_in_cd_rods.add_nuclide('Cd108', 2.4166e-05)
+ag_in_cd_rods.add_nuclide('Cd110', 3.3936e-04)
+ag_in_cd_rods.add_nuclide('Cd111', 3.4821e-04)
+ag_in_cd_rods.add_nuclide('Cd112', 6.5611e-04)
+ag_in_cd_rods.add_nuclide('Cd113', 3.3275e-04)
+ag_in_cd_rods.add_nuclide('Cd114', 7.8252e-04)
+ag_in_cd_rods.add_nuclide('Cd116', 2.0443e-04)
+ag_in_cd_rods.add_nuclide('In113', 3.4219e-04)
+ag_in_cd_rods.add_nuclide('In115', 7.6511e-03)
+
+b4c = openmc.Material(name='B4C Control Rod')
+b4c.set_density('g/cm3', 1.76)
+b4c.temperature = 300
+b4c.add_nuclide('B10', 1.5206e-02)
+b4c.add_nuclide('B11', 6.1514e-02)
+b4c.add_nuclide('C12', 1.8972e-02)
+b4c.add_nuclide('C13', 2.1252e-04)
+
+helium = openmc.Material(name='Helium')
+helium.set_density('g/cm3', 0.0015981)
+helium.temperature = 300
+helium.add_nuclide('He3', 4.8089e-10)
+helium.add_nuclide('He4', 2.4044e-04)
+
+inconel_718 = openmc.Material(name='Inconel 718')
+inconel_718.set_density('g/cm3', 8.2)
+inconel_718.temperature = 300
+inconel_718.add_nuclide('Cr50', 7.8239e-04)
+inconel_718.add_nuclide('Cr52', 1.5088e-02)
+inconel_718.add_nuclide('Cr53', 1.7108e-03)
+inconel_718.add_nuclide('Cr54', 4.2586e-04)
+inconel_718.add_nuclide('Fe54', 1.4797e-03)
+inconel_718.add_nuclide('Fe56', 2.3229e-02)
+inconel_718.add_nuclide('Fe57', 5.3645e-04)
+inconel_718.add_nuclide('Fe58', 7.1392e-05)
+inconel_718.add_nuclide('Mn55', 7.8201e-04)
+inconel_718.add_nuclide('Ni58', 2.9320e-02)
+inconel_718.add_nuclide('Ni60', 1.1294e-02)
+inconel_718.add_nuclide('Ni61', 4.9094e-04)
+inconel_718.add_nuclide('Ni62', 1.5653e-03)
+inconel_718.add_nuclide('Ni64', 3.9864e-04)
+inconel_718.add_nuclide('Si28', 5.6757e-04)
+inconel_718.add_nuclide('Si29', 2.8820e-05)
+inconel_718.add_nuclide('Si30', 1.8998e-05)
+
+ss304 = openmc.Material(name='Stainless Steel 304')
+ss304.set_density('g/cm3', 8.03)
+ss304.temperature = 300
+ss304.add_nuclide('Cr50', 7.6778e-04)
+ss304.add_nuclide('Cr52', 1.4806e-02)
+ss304.add_nuclide('Cr53', 1.6789e-03)
+ss304.add_nuclide('Cr54', 4.1791e-04)
+ss304.add_nuclide('Fe54', 3.4620e-03)
+ss304.add_nuclide('Fe56', 5.4345e-02)
+ss304.add_nuclide('Fe57', 1.2551e-03)
+ss304.add_nuclide('Fe58', 1.6703e-04)
+ss304.add_nuclide('Mn55', 1.7604e-03)
+ss304.add_nuclide('Ni58', 5.6089e-03)
+ss304.add_nuclide('Ni60', 2.1605e-03)
+ss304.add_nuclide('Ni61', 9.3917e-05)
+ss304.add_nuclide('Ni62', 2.9945e-04)
+ss304.add_nuclide('Ni64', 7.6261e-05)
+ss304.add_nuclide('Si28', 9.5281e-04)
+ss304.add_nuclide('Si29', 4.8381e-05)
+ss304.add_nuclide('Si30', 3.1893e-05)
+
+zirc_4 = openmc.Material()
+zirc_4.set_density('g/cm3', 6.55)
+zirc_4.temperature = 300
+zirc_4.add_nuclide('Cr50', 3.2962e-06)
+zirc_4.add_nuclide('Cr52', 6.3564e-05)
+zirc_4.add_nuclide('Cr53', 7.2076e-06)
+zirc_4.add_nuclide('Cr54', 1.7941e-06)
+zirc_4.add_nuclide('Fe54', 8.6698e-06)
+zirc_4.add_nuclide('Fe56', 1.3610e-04)
+zirc_4.add_nuclide('Fe57', 3.1431e-06)
+zirc_4.add_nuclide('Fe58', 4.1829e-07)
+zirc_4.add_nuclide('O16', 3.0744e-04)
+zirc_4.add_nuclide('O17', 1.1680e-07)
+zirc_4.add_nuclide('O18', 6.1648e-07)
+zirc_4.add_nuclide('Sn112', 4.6735e-06)
+zirc_4.add_nuclide('Sn114', 3.1799e-06)
+zirc_4.add_nuclide('Sn115', 1.6381e-06)
+zirc_4.add_nuclide('Sn116', 7.0055e-05)
+zirc_4.add_nuclide('Sn117', 3.7003e-05)
+zirc_4.add_nuclide('Sn118', 1.1669e-04)
+zirc_4.add_nuclide('Sn119', 4.1387e-05)
+zirc_4.add_nuclide('Sn120', 1.5697e-04)
+zirc_4.add_nuclide('Sn122', 2.2308e-05)
+zirc_4.add_nuclide('Sn124', 2.7897e-05)
+zirc_4.add_nuclide('Zr90', 2.1828e-02)
+zirc_4.add_nuclide('Zr91', 4.7601e-03)
+zirc_4.add_nuclide('Zr92', 7.2759e-03)
+zirc_4.add_nuclide('Zr94', 7.3734e-03)
+zirc_4.add_nuclide('Zr96', 1.1879e-03)
+
+borated_water = openmc.Material()
+borated_water.set_density('g/cm3', 0.740582068)
+borated_water.temperature = 300
+borated_water.add_nuclide('B10', 7.9714e-06)
+borated_water.add_nuclide('B11', 3.2247e-05)
+borated_water.add_nuclide('H1', 4.9456e-02)
+borated_water.add_nuclide('H2', 7.7035e-06)
+borated_water.add_nuclide('O16', 2.4673e-02)
+borated_water.add_nuclide('O17', 9.3734e-06)
+borated_water.add_nuclide('O18', 4.9474e-05)
+
+plate_borated_water = openmc.Material()
+plate_borated_water.set_density('g/cm3', 0.981002532)
+plate_borated_water.add_nuclide('B10', 1.0559e-05)
+plate_borated_water.add_nuclide('B11', 4.2716e-05)
+plate_borated_water.add_nuclide('H1', 6.5512e-02)
+plate_borated_water.add_nuclide('H2', 1.0204e-05)
+plate_borated_water.add_nuclide('O16', 3.2683e-02)
+plate_borated_water.add_nuclide('O17', 1.2416e-05)
+plate_borated_water.add_nuclide('O18', 6.5535e-05)
+
+plate_stainless = openmc.Material()
+plate_stainless.set_density('g/cm3', 3.68384807)
+plate_stainless.temperature = 300
+plate_stainless.add_nuclide('Cr50', 3.5223e-04)
+plate_stainless.add_nuclide('Cr52', 6.7924e-03)
+plate_stainless.add_nuclide('Cr53', 7.7020e-04)
+plate_stainless.add_nuclide('Cr54', 1.9172e-04)
+plate_stainless.add_nuclide('Fe54', 1.5882e-03)
+plate_stainless.add_nuclide('Fe56', 2.4931e-02)
+plate_stainless.add_nuclide('Fe57', 5.7578e-04)
+plate_stainless.add_nuclide('Fe58', 7.6625e-05)
+plate_stainless.add_nuclide('Mn55', 8.0762e-04)
+plate_stainless.add_nuclide('Ni58', 2.5731e-03)
+plate_stainless.add_nuclide('Ni60', 9.9117e-04)
+plate_stainless.add_nuclide('Ni61', 4.3085e-05)
+plate_stainless.add_nuclide('Ni62', 1.3738e-04)
+plate_stainless.add_nuclide('Ni64', 3.4985e-05)
+plate_stainless.add_nuclide('Si28', 4.3711e-04)
+plate_stainless.add_nuclide('Si29', 2.2195e-05)
+plate_stainless.add_nuclide('Si30', 1.4631e-05)
+
+carbon_steel = openmc.Material()
+carbon_steel.set_density('g/cm3', 7.8)
+carbon_steel.temperature = 300
+carbon_steel.add_nuclide('Al27', 4.3523e-05)
+carbon_steel.add_nuclide('B10', 2.5833e-06)
+carbon_steel.add_nuclide('B11', 1.0450e-05)
+carbon_steel.add_nuclide('C12', 1.0442e-03)
+carbon_steel.add_nuclide('C13', 1.1697e-05)
+carbon_steel.add_nuclide('Ca40', 1.7043e-05)
+carbon_steel.add_nuclide('Ca42', 1.1375e-07)
+carbon_steel.add_nuclide('Ca43', 2.3734e-08)
+carbon_steel.add_nuclide('Ca44', 3.6673e-07)
+carbon_steel.add_nuclide('Ca46', 7.0322e-10)
+carbon_steel.add_nuclide('Ca48', 3.2875e-08)
+carbon_steel.add_nuclide('Cr50', 1.3738e-05)
+carbon_steel.add_nuclide('Cr52', 2.6493e-04)
+carbon_steel.add_nuclide('Cr53', 3.0041e-05)
+carbon_steel.add_nuclide('Cr54', 7.4778e-06)
+carbon_steel.add_nuclide('Cu63', 1.0223e-04)
+carbon_steel.add_nuclide('Cu65', 4.5608e-05)
+carbon_steel.add_nuclide('Fe54', 4.7437e-03)
+carbon_steel.add_nuclide('Fe56', 7.4465e-02)
+carbon_steel.add_nuclide('Fe57', 1.7197e-03)
+carbon_steel.add_nuclide('Fe58', 2.2886e-04)
+carbon_steel.add_nuclide('Mn55', 6.4126e-04)
+carbon_steel.add_nuclide('Mo100', 2.9814e-05)
+carbon_steel.add_nuclide('Mo92', 4.4822e-05)
+carbon_steel.add_nuclide('Mo94', 2.8110e-05)
+carbon_steel.add_nuclide('Mo95', 4.8567e-05)
+carbon_steel.add_nuclide('Mo96', 5.1015e-05)
+carbon_steel.add_nuclide('Mo97', 2.9319e-05)
+carbon_steel.add_nuclide('Mo98', 7.4327e-05)
+carbon_steel.add_nuclide('Nb93', 5.0559e-06)
+carbon_steel.add_nuclide('Ni58', 4.0862e-04)
+carbon_steel.add_nuclide('Ni60', 1.5740e-04)
+carbon_steel.add_nuclide('Ni61', 6.8420e-06)
+carbon_steel.add_nuclide('Ni62', 2.1815e-05)
+carbon_steel.add_nuclide('Ni64', 5.5557e-06)
+carbon_steel.add_nuclide('P31', 3.7913e-05)
+carbon_steel.add_nuclide('S32', 3.4808e-05)
+carbon_steel.add_nuclide('S33', 2.7420e-07)
+carbon_steel.add_nuclide('S34', 1.5368e-06)
+carbon_steel.add_nuclide('S36', 5.3398e-09)
+carbon_steel.add_nuclide('Si28', 6.1702e-04)
+carbon_steel.add_nuclide('Si29', 3.1330e-05)
+carbon_steel.add_nuclide('Si30', 2.0653e-05)
+carbon_steel.add_nuclide('Ti46', 1.2144e-06)
+carbon_steel.add_nuclide('Ti47', 1.0952e-06)
+carbon_steel.add_nuclide('Ti48', 1.0851e-05)
+carbon_steel.add_nuclide('Ti49', 7.9634e-07)
+carbon_steel.add_nuclide('Ti50', 7.6249e-07)
+carbon_steel.add_nuclide('V50', 1.1526e-07)
+carbon_steel.add_nuclide('V51', 4.5989e-05)
+
+"""
+materials += [fuel_16, fuel_24, fuel_31, fuel_32, fuel_34]
+"""
+materials += [air, borosilicate_glass, ag_in_cd_rods, b4c, helium,
+              inconel_718, ss304, zirc_4, borated_water, plate_borated_water,
+              plate_stainless, carbon_steel]
+materials.cross_sections="/opt/xdata/cross_sections.xml"
+materials.export_to_xml()
+
+## Fuel Pin
+fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
+clad_ir = openmc.ZCylinder(x0=0, y0=0, r=0.40005, name='Clad IR')
+clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
+
+## Upper Fuel Pin Plenum
+spring = openmc.ZCylinder(x0=0, y0=0, r=0.06459, name='Spring')
+clad_ir = openmc.ZCylinder(x0=0, y0=0, r=0.40005, name='Clad IR')
+clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
+
+## Empty Guide Tube above Dashpot
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_or = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide OR')
+
+## Empty Guide Tube at Dashpot
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.50419, name='Guide IR')
+guide_or = openmc.ZCylinder(x0=0, y0=0, r=0.54610, name='Guide OR')
+
+## Instrument Tube Pin Geometry
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
+## Bare Instrument Thimble Pin
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+
+## BP above Dashpoint
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.21400, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.23051, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.24130, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.42672, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
+## BP Plennum
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.21400, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.23051, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.43688, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.50419, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.54610, name='Guide IR')
+
+## Control Rod Pin Upper Geometry
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.37338, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
+## Control Rod Pin Lower Geometry
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38227, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
+## Control Rod Pin Spacer
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.37845, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
+## Control Rod Plenum Geometry
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.06459, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.38608, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.48387, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.56134, name='Guide IR')
+guide_ir = openmc.ZCylinder(x0=0, y0=0, r=0.60198, name='Guide IR')
+
diff --git a/own_make/pin_cell.py b/own_make/pin_cell.py
new file mode 100644
index 0000000..2a3adb6
--- /dev/null
+++ b/own_make/pin_cell.py
@@ -0,0 +1,653 @@
+from numbers import Integral
+
+import numpy as np
+
+import openmc
+import openmc.model
+
+
+def pwr_pin_cell():
+    """Create a PWR pin-cell model.
+
+    This model is a single fuel pin with 2.4 w/o enriched UO2 corresponding to a
+    beginning-of-cycle condition and borated water. The specifications are from
+    the `BEAVRS <http://crpg.mit.edu/research/beavrs>`_ benchmark. Note that the
+    number of particles/batches is initially set very low for testing purposes.
+
+    Returns
+    -------
+    model : openmc.model.Model
+        A PWR pin-cell model
+
+    """
+    model = openmc.model.Model()
+
+    # Define materials.
+    fuel = openmc.Material(name='UO2 (2.4%)')
+    fuel.set_density('g/cm3', 10.29769)
+    fuel.add_nuclide('U234', 4.4843e-6)
+    fuel.add_nuclide('U235', 5.5815e-4)
+    fuel.add_nuclide('U238', 2.2408e-2)
+    fuel.add_nuclide('O16', 4.5829e-2)
+
+    clad = openmc.Material(name='Zircaloy')
+    clad.set_density('g/cm3', 6.55)
+    clad.add_nuclide('Zr90', 2.1827e-2)
+    clad.add_nuclide('Zr91', 4.7600e-3)
+    clad.add_nuclide('Zr92', 7.2758e-3)
+    clad.add_nuclide('Zr94', 7.3734e-3)
+    clad.add_nuclide('Zr96', 1.1879e-3)
+
+    hot_water = openmc.Material(name='Hot borated water')
+    hot_water.set_density('g/cm3', 0.740582)
+    hot_water.add_nuclide('H1', 4.9457e-2)
+    hot_water.add_nuclide('O16', 2.4672e-2)
+    hot_water.add_nuclide('B10', 8.0042e-6)
+    hot_water.add_nuclide('B11', 3.2218e-5)
+    hot_water.add_s_alpha_beta('c_H_in_H2O')
+
+    # Define the materials file.
+    model.materials = (fuel, clad, hot_water)
+
+    # Instantiate ZCylinder surfaces
+    pitch = 1.26
+    fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
+    clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
+    left = openmc.XPlane(x0=-pitch/2, name='left', boundary_type='reflective')
+    right = openmc.XPlane(x0=pitch/2, name='right', boundary_type='reflective')
+    bottom = openmc.YPlane(y0=-pitch/2, name='bottom',
+                           boundary_type='reflective')
+    top = openmc.YPlane(y0=pitch/2, name='top', boundary_type='reflective')
+
+    # Instantiate Cells
+    fuel_pin = openmc.Cell(name='Fuel', fill=fuel)
+    cladding = openmc.Cell(name='Cladding', fill=clad)
+    water = openmc.Cell(name='Water', fill=hot_water)
+
+    # Use surface half-spaces to define regions
+    fuel_pin.region = -fuel_or
+    cladding.region = +fuel_or & -clad_or
+    water.region = +clad_or & +left & -right & +bottom & -top
+
+    # Create root universe
+    model.geometry.root_universe = openmc.Universe(0, name='root universe')
+    model.geometry.root_universe.add_cells([fuel_pin, cladding, water])
+
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    model.settings.particles = 100
+    model.settings.source = openmc.Source(space=openmc.stats.Box(
+        [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
+
+    plot = openmc.Plot.from_geometry(model.geometry)
+    plot.pixels = (300, 300)
+    plot.color_by = 'material'
+    model.plots.append(plot)
+
+    return model
+
+
+
+def pwr_core():
+    """Create a PWR full-core model.
+
+    This model is the OECD/NEA Monte Carlo Performance benchmark which is a
+    grossly simplified pressurized water reactor (PWR) with 241 fuel
+    assemblies. Note that the number of particles/batches is initially set very
+    low for testing purposes.
+
+    Returns
+    -------
+    model : openmc.model.Model
+        Full-core PWR model
+
+    """
+    model = openmc.model.Model()
+
+    # Define materials.
+    fuel = openmc.Material(1, name='UOX fuel')
+    fuel.set_density('g/cm3', 10.062)
+    fuel.add_nuclide('U234', 4.9476e-6)
+    fuel.add_nuclide('U235', 4.8218e-4)
+    fuel.add_nuclide('U238', 2.1504e-2)
+    fuel.add_nuclide('Xe135', 1.0801e-8)
+    fuel.add_nuclide('O16', 4.5737e-2)
+
+    clad = openmc.Material(2, name='Zircaloy')
+    clad.set_density('g/cm3', 5.77)
+    clad.add_nuclide('Zr90', 0.5145)
+    clad.add_nuclide('Zr91', 0.1122)
+    clad.add_nuclide('Zr92', 0.1715)
+    clad.add_nuclide('Zr94', 0.1738)
+    clad.add_nuclide('Zr96', 0.0280)
+
+    cold_water = openmc.Material(3, name='Cold borated water')
+    cold_water.set_density('atom/b-cm', 0.07416)
+    cold_water.add_nuclide('H1', 2.0)
+    cold_water.add_nuclide('O16', 1.0)
+    cold_water.add_nuclide('B10', 6.490e-4)
+    cold_water.add_nuclide('B11', 2.689e-3)
+    cold_water.add_s_alpha_beta('c_H_in_H2O')
+
+    hot_water = openmc.Material(4, name='Hot borated water')
+    hot_water.set_density('atom/b-cm', 0.06614)
+    hot_water.add_nuclide('H1', 2.0)
+    hot_water.add_nuclide('O16', 1.0)
+    hot_water.add_nuclide('B10', 6.490e-4)
+    hot_water.add_nuclide('B11', 2.689e-3)
+    hot_water.add_s_alpha_beta('c_H_in_H2O')
+
+    rpv_steel = openmc.Material(5, name='Reactor pressure vessel steel')
+    rpv_steel.set_density('g/cm3', 7.9)
+    rpv_steel.add_nuclide('Fe54', 0.05437098, 'wo')
+    rpv_steel.add_nuclide('Fe56', 0.88500663, 'wo')
+    rpv_steel.add_nuclide('Fe57', 0.0208008, 'wo')
+    rpv_steel.add_nuclide('Fe58', 0.00282159, 'wo')
+    rpv_steel.add_nuclide('Ni58', 0.0067198, 'wo')
+    rpv_steel.add_nuclide('Ni60', 0.0026776, 'wo')
+    rpv_steel.add_nuclide('Mn55', 0.01, 'wo')
+    rpv_steel.add_nuclide('Cr52', 0.002092475, 'wo')
+    rpv_steel.add_nuclide('C0', 0.0025, 'wo')
+    rpv_steel.add_nuclide('Cu63', 0.0013696, 'wo')
+
+    lower_rad_ref = openmc.Material(6, name='Lower radial reflector')
+    lower_rad_ref.set_density('g/cm3', 4.32)
+    lower_rad_ref.add_nuclide('H1', 0.0095661, 'wo')
+    lower_rad_ref.add_nuclide('O16', 0.0759107, 'wo')
+    lower_rad_ref.add_nuclide('B10', 3.08409e-5, 'wo')
+    lower_rad_ref.add_nuclide('B11', 1.40499e-4, 'wo')
+    lower_rad_ref.add_nuclide('Fe54', 0.035620772088, 'wo')
+    lower_rad_ref.add_nuclide('Fe56', 0.579805982228, 'wo')
+    lower_rad_ref.add_nuclide('Fe57', 0.01362750048, 'wo')
+    lower_rad_ref.add_nuclide('Fe58', 0.001848545204, 'wo')
+    lower_rad_ref.add_nuclide('Ni58', 0.055298376566, 'wo')
+    lower_rad_ref.add_nuclide('Mn55', 0.0182870, 'wo')
+    lower_rad_ref.add_nuclide('Cr52', 0.145407678031, 'wo')
+    lower_rad_ref.add_s_alpha_beta('c_H_in_H2O')
+
+    upper_rad_ref = openmc.Material(7, name='Upper radial reflector / Top plate region')
+    upper_rad_ref.set_density('g/cm3', 4.28)
+    upper_rad_ref.add_nuclide('H1', 0.0086117, 'wo')
+    upper_rad_ref.add_nuclide('O16', 0.0683369, 'wo')
+    upper_rad_ref.add_nuclide('B10', 2.77638e-5, 'wo')
+    upper_rad_ref.add_nuclide('B11', 1.26481e-4, 'wo')
+    upper_rad_ref.add_nuclide('Fe54', 0.035953677186, 'wo')
+    upper_rad_ref.add_nuclide('Fe56', 0.585224740891, 'wo')
+    upper_rad_ref.add_nuclide('Fe57', 0.01375486056, 'wo')
+    upper_rad_ref.add_nuclide('Fe58', 0.001865821363, 'wo')
+    upper_rad_ref.add_nuclide('Ni58', 0.055815129186, 'wo')
+    upper_rad_ref.add_nuclide('Mn55', 0.0184579, 'wo')
+    upper_rad_ref.add_nuclide('Cr52', 0.146766614995, 'wo')
+    upper_rad_ref.add_s_alpha_beta('c_H_in_H2O')
+
+    bot_plate = openmc.Material(8, name='Bottom plate region')
+    bot_plate.set_density('g/cm3', 7.184)
+    bot_plate.add_nuclide('H1', 0.0011505, 'wo')
+    bot_plate.add_nuclide('O16', 0.0091296, 'wo')
+    bot_plate.add_nuclide('B10', 3.70915e-6, 'wo')
+    bot_plate.add_nuclide('B11', 1.68974e-5, 'wo')
+    bot_plate.add_nuclide('Fe54', 0.03855611055, 'wo')
+    bot_plate.add_nuclide('Fe56', 0.627585036425, 'wo')
+    bot_plate.add_nuclide('Fe57', 0.014750478, 'wo')
+    bot_plate.add_nuclide('Fe58', 0.002000875025, 'wo')
+    bot_plate.add_nuclide('Ni58', 0.059855207342, 'wo')
+    bot_plate.add_nuclide('Mn55', 0.0197940, 'wo')
+    bot_plate.add_nuclide('Cr52', 0.157390026871, 'wo')
+    bot_plate.add_s_alpha_beta('c_H_in_H2O')
+
+    bot_nozzle = openmc.Material(9, name='Bottom nozzle region')
+    bot_nozzle.set_density('g/cm3', 2.53)
+    bot_nozzle.add_nuclide('H1', 0.0245014, 'wo')
+    bot_nozzle.add_nuclide('O16', 0.1944274, 'wo')
+    bot_nozzle.add_nuclide('B10', 7.89917e-5, 'wo')
+    bot_nozzle.add_nuclide('B11', 3.59854e-4, 'wo')
+    bot_nozzle.add_nuclide('Fe54', 0.030411411144, 'wo')
+    bot_nozzle.add_nuclide('Fe56', 0.495012237964, 'wo')
+    bot_nozzle.add_nuclide('Fe57', 0.01163454624, 'wo')
+    bot_nozzle.add_nuclide('Fe58', 0.001578204652, 'wo')
+    bot_nozzle.add_nuclide('Ni58', 0.047211231662, 'wo')
+    bot_nozzle.add_nuclide('Mn55', 0.0156126, 'wo')
+    bot_nozzle.add_nuclide('Cr52', 0.124142524198, 'wo')
+    bot_nozzle.add_s_alpha_beta('c_H_in_H2O')
+
+    top_nozzle = openmc.Material(10, name='Top nozzle region')
+    top_nozzle.set_density('g/cm3', 1.746)
+    top_nozzle.add_nuclide('H1', 0.0358870, 'wo')
+    top_nozzle.add_nuclide('O16', 0.2847761, 'wo')
+    top_nozzle.add_nuclide('B10', 1.15699e-4, 'wo')
+    top_nozzle.add_nuclide('B11', 5.27075e-4, 'wo')
+    top_nozzle.add_nuclide('Fe54', 0.02644016154, 'wo')
+    top_nozzle.add_nuclide('Fe56', 0.43037146399, 'wo')
+    top_nozzle.add_nuclide('Fe57', 0.0101152584, 'wo')
+    top_nozzle.add_nuclide('Fe58', 0.00137211607, 'wo')
+    top_nozzle.add_nuclide('Ni58', 0.04104621835, 'wo')
+    top_nozzle.add_nuclide('Mn55', 0.0135739, 'wo')
+    top_nozzle.add_nuclide('Cr52', 0.107931450781, 'wo')
+    top_nozzle.add_s_alpha_beta('c_H_in_H2O')
+
+    top_fa = openmc.Material(11, name='Top of fuel assemblies')
+    top_fa.set_density('g/cm3', 3.044)
+    top_fa.add_nuclide('H1', 0.0162913, 'wo')
+    top_fa.add_nuclide('O16', 0.1292776, 'wo')
+    top_fa.add_nuclide('B10', 5.25228e-5, 'wo')
+    top_fa.add_nuclide('B11', 2.39272e-4, 'wo')
+    top_fa.add_nuclide('Zr90', 0.43313403903, 'wo')
+    top_fa.add_nuclide('Zr91', 0.09549277374, 'wo')
+    top_fa.add_nuclide('Zr92', 0.14759527104, 'wo')
+    top_fa.add_nuclide('Zr94', 0.15280552077, 'wo')
+    top_fa.add_nuclide('Zr96', 0.02511169542, 'wo')
+    top_fa.add_s_alpha_beta('c_H_in_H2O')
+
+    bot_fa = openmc.Material(12, name='Bottom of fuel assemblies')
+    bot_fa.set_density('g/cm3', 1.762)
+    bot_fa.add_nuclide('H1', 0.0292856, 'wo')
+    bot_fa.add_nuclide('O16', 0.2323919, 'wo')
+    bot_fa.add_nuclide('B10', 9.44159e-5, 'wo')
+    bot_fa.add_nuclide('B11', 4.30120e-4, 'wo')
+    bot_fa.add_nuclide('Zr90', 0.3741373658, 'wo')
+    bot_fa.add_nuclide('Zr91', 0.0824858164, 'wo')
+    bot_fa.add_nuclide('Zr92', 0.1274914944, 'wo')
+    bot_fa.add_nuclide('Zr94', 0.1319920622, 'wo')
+    bot_fa.add_nuclide('Zr96', 0.0216912612, 'wo')
+    bot_fa.add_s_alpha_beta('c_H_in_H2O')
+
+    # Define the materials file.
+    model.materials = (fuel, clad, cold_water, hot_water, rpv_steel,
+                       lower_rad_ref, upper_rad_ref, bot_plate,
+                       bot_nozzle, top_nozzle, top_fa, bot_fa)
+
+    # Define surfaces.
+    s1 = openmc.ZCylinder(r=0.41, surface_id=1)
+    s2 = openmc.ZCylinder(r=0.475, surface_id=2)
+    s3 = openmc.ZCylinder(r=0.56, surface_id=3)
+    s4 = openmc.ZCylinder(r=0.62, surface_id=4)
+    s5 = openmc.ZCylinder(r=187.6, surface_id=5)
+    s6 = openmc.ZCylinder(r=209.0, surface_id=6)
+    s7 = openmc.ZCylinder(r=229.0, surface_id=7)
+    s8 = openmc.ZCylinder(r=249.0, surface_id=8, boundary_type='vacuum')
+
+    s31 = openmc.ZPlane(z0=-229.0, surface_id=31, boundary_type='vacuum')
+    s32 = openmc.ZPlane(z0=-199.0, surface_id=32)
+    s33 = openmc.ZPlane(z0=-193.0, surface_id=33)
+    s34 = openmc.ZPlane(z0=-183.0, surface_id=34)
+    s35 = openmc.ZPlane(z0=0.0, surface_id=35)
+    s36 = openmc.ZPlane(z0=183.0, surface_id=36)
+    s37 = openmc.ZPlane(z0=203.0, surface_id=37)
+    s38 = openmc.ZPlane(z0=215.0, surface_id=38)
+    s39 = openmc.ZPlane(z0=223.0, surface_id=39, boundary_type='vacuum')
+
+    # Define pin cells.
+    fuel_cold = openmc.Universe(name='Fuel pin, cladding, cold water',
+                                universe_id=1)
+    c21 = openmc.Cell(cell_id=21, fill=fuel, region=-s1)
+    c22 = openmc.Cell(cell_id=22, fill=clad, region=+s1 & -s2)
+    c23 = openmc.Cell(cell_id=23, fill=cold_water, region=+s2)
+    fuel_cold.add_cells((c21, c22, c23))
+
+    tube_cold = openmc.Universe(name='Instrumentation guide tube, '
+                                'cold water', universe_id=2)
+    c24 = openmc.Cell(cell_id=24, fill=cold_water, region=-s3)
+    c25 = openmc.Cell(cell_id=25, fill=clad, region=+s3 & -s4)
+    c26 = openmc.Cell(cell_id=26, fill=cold_water, region=+s4)
+    tube_cold.add_cells((c24, c25, c26))
+
+    fuel_hot = openmc.Universe(name='Fuel pin, cladding, hot water',
+                               universe_id=3)
+    c27 = openmc.Cell(cell_id=27, fill=fuel, region=-s1)
+    c28 = openmc.Cell(cell_id=28, fill=clad, region=+s1 & -s2)
+    c29 = openmc.Cell(cell_id=29, fill=hot_water, region=+s2)
+    fuel_hot.add_cells((c27, c28, c29))
+
+    tube_hot = openmc.Universe(name='Instrumentation guide tube, hot water',
+                               universe_id=4)
+    c30 = openmc.Cell(cell_id=30, fill=hot_water, region=-s3)
+    c31 = openmc.Cell(cell_id=31, fill=clad, region=+s3 & -s4)
+    c32 = openmc.Cell(cell_id=32, fill=hot_water, region=+s4)
+    tube_hot.add_cells((c30, c31, c32))
+
+    # Set positions occupied by guide tubes
+    tube_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8, 11, 14,
+                       2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
+    tube_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8, 8,
+                       11, 11, 11, 11, 11, 13, 13, 14, 14, 14])
+
+    # Define fuel lattices.
+    l100 = openmc.RectLattice(name='Fuel assembly (lower half)', lattice_id=100)
+    l100.lower_left = (-10.71, -10.71)
+    l100.pitch = (1.26, 1.26)
+    l100.universes = np.tile(fuel_cold, (17, 17))
+    l100.universes[tube_x, tube_y] = tube_cold
+
+    l101 = openmc.RectLattice(name='Fuel assembly (upper half)', lattice_id=101)
+    l101.lower_left = (-10.71, -10.71)
+    l101.pitch = (1.26, 1.26)
+    l101.universes = np.tile(fuel_hot, (17, 17))
+    l101.universes[tube_x, tube_y] = tube_hot
+
+    # Define assemblies.
+    fa_cw = openmc.Universe(name='Water assembly (cold)', universe_id=5)
+    c50 = openmc.Cell(cell_id=50, fill=cold_water, region=+s34 & -s35)
+    fa_cw.add_cell(c50)
+
+    fa_hw = openmc.Universe(name='Water assembly (hot)', universe_id=7)
+    c70 = openmc.Cell(cell_id=70, fill=hot_water, region=+s35 & -s36)
+    fa_hw.add_cell(c70)
+
+    fa_cold = openmc.Universe(name='Fuel assembly (cold)', universe_id=6)
+    c60 = openmc.Cell(cell_id=60, fill=l100, region=+s34 & -s35)
+    fa_cold.add_cell(c60)
+
+    fa_hot = openmc.Universe(name='Fuel assembly (hot)', universe_id=8)
+    c80 = openmc.Cell(cell_id=80, fill=l101, region=+s35 & -s36)
+    fa_hot.add_cell(c80)
+
+    # Define core lattices
+    l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200)
+    l200.lower_left = (-224.91, -224.91)
+    l200.pitch = (21.42, 21.42)
+    l200.universes = [
+        [fa_cw]*21,
+        [fa_cw]*21,
+        [fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
+        [fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
+        [fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
+        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
+        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
+        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
+        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
+        [fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
+        [fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
+        [fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
+        [fa_cw]*21,
+        [fa_cw]*21]
+
+    l201 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=201)
+    l201.lower_left = (-224.91, -224.91)
+    l201.pitch = (21.42, 21.42)
+    l201.universes = [
+        [fa_hw]*21,
+        [fa_hw]*21,
+        [fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
+        [fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
+        [fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
+        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
+        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
+        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
+        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
+        [fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
+        [fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
+        [fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
+        [fa_hw]*21,
+        [fa_hw]*21]
+
+    # Define root universe.
+    root = openmc.Universe(universe_id=0, name='root universe')
+    c1 = openmc.Cell(cell_id=1, fill=l200, region=-s6 & +s34 & -s35)
+    c2 = openmc.Cell(cell_id=2, fill=l201, region=-s6 & +s35 & -s36)
+    c3 = openmc.Cell(cell_id=3, fill=bot_plate, region=-s7 & +s31 & -s32)
+    c4 = openmc.Cell(cell_id=4, fill=bot_nozzle, region=-s5 & +s32 & -s33)
+    c5 = openmc.Cell(cell_id=5, fill=bot_fa, region=-s5 & +s33 & -s34)
+    c6 = openmc.Cell(cell_id=6, fill=top_fa, region=-s5 & +s36 & -s37)
+    c7 = openmc.Cell(cell_id=7, fill=top_nozzle, region=-s5 & +s37 & -s38)
+    c8 = openmc.Cell(cell_id=8, fill=upper_rad_ref, region=-s7 & +s38 & -s39)
+    c9 = openmc.Cell(cell_id=9, fill=bot_nozzle, region=+s6 & -s7 & +s32 & -s38)
+    c10 = openmc.Cell(cell_id=10, fill=rpv_steel, region=+s7 & -s8 & +s31 & -s39)
+    c11 = openmc.Cell(cell_id=11, fill=lower_rad_ref, region=+s5 & -s6 & +s32 & -s34)
+    c12 = openmc.Cell(cell_id=12, fill=upper_rad_ref, region=+s5 & -s6 & +s36 & -s38)
+    root.add_cells((c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12))
+
+    # Assign root universe to geometry
+    model.geometry.root_universe = root
+
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    model.settings.particles = 100
+    model.settings.source = openmc.Source(space=openmc.stats.Box(
+        [-160, -160, -183], [160, 160, 183]))
+
+    plot = openmc.Plot()
+    # plot.basis = 'yz'
+    plot.type = 'voxel'
+    plot.origin = (0, 0, 0)
+    plot.width = (500, 500, 500)
+    plot.pixels = (250, 250, 250)
+    # plot.color_by = 'material'
+    model.plots.append(plot)
+
+    return model
+
+
+
+def pwr_assembly():
+    """Create a PWR assembly model.
+
+    This model is a reflected 17x17 fuel assembly from the the `BEAVRS
+    <http://crpg.mit.edu/research/beavrs>`_ benchmark. The fuel is 2.4 w/o
+    enriched UO2 corresponding to a beginning-of-cycle condition. Note that the
+    number of particles/batches is initially set very low for testing purposes.
+
+    Returns
+    -------
+    model : openmc.model.Model
+        A PWR assembly model
+
+    """
+
+    model = openmc.model.Model()
+
+    # Define materials.
+    fuel = openmc.Material(name='Fuel')
+    fuel.set_density('g/cm3', 10.29769)
+    fuel.add_nuclide('U234', 4.4843e-6)
+    fuel.add_nuclide('U235', 5.5815e-4)
+    fuel.add_nuclide('U238', 2.2408e-2)
+    fuel.add_nuclide('O16', 4.5829e-2)
+
+    clad = openmc.Material(name='Cladding')
+    clad.set_density('g/cm3', 6.55)
+    clad.add_nuclide('Zr90', 2.1827e-2)
+    clad.add_nuclide('Zr91', 4.7600e-3)
+    clad.add_nuclide('Zr92', 7.2758e-3)
+    clad.add_nuclide('Zr94', 7.3734e-3)
+    clad.add_nuclide('Zr96', 1.1879e-3)
+
+    hot_water = openmc.Material(name='Hot borated water')
+    hot_water.set_density('g/cm3', 0.740582)
+    hot_water.add_nuclide('H1', 4.9457e-2)
+    hot_water.add_nuclide('O16', 2.4672e-2)
+    hot_water.add_nuclide('B10', 8.0042e-6)
+    hot_water.add_nuclide('B11', 3.2218e-5)
+    hot_water.add_s_alpha_beta('c_H_in_H2O')
+
+    # Define the materials file.
+    model.materials = (fuel, clad, hot_water)
+
+    # Instantiate ZCylinder surfaces
+    fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR')
+    clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR')
+
+    # Create boundary planes to surround the geometry
+    pitch = 21.42
+    min_x = openmc.XPlane(x0=-pitch/2, boundary_type='reflective')
+    max_x = openmc.XPlane(x0=+pitch/2, boundary_type='reflective')
+    min_y = openmc.YPlane(y0=-pitch/2, boundary_type='reflective')
+    max_y = openmc.YPlane(y0=+pitch/2, boundary_type='reflective')
+
+    # Create a fuel pin universe
+    fuel_pin_universe = openmc.Universe(name='Fuel Pin')
+    fuel_cell = openmc.Cell(name='fuel', fill=fuel, region=-fuel_or)
+    clad_cell = openmc.Cell(name='clad', fill=clad, region=+fuel_or & -clad_or)
+    hot_water_cell = openmc.Cell(name='hot water', fill=hot_water, region=+clad_or)
+    fuel_pin_universe.add_cells([fuel_cell, clad_cell, hot_water_cell])
+
+
+    # Create a control rod guide tube universe
+    guide_tube_universe = openmc.Universe(name='Guide Tube')
+    gt_inner_cell = openmc.Cell(name='guide tube inner water', fill=hot_water,
+                                region=-fuel_or)
+    gt_clad_cell = openmc.Cell(name='guide tube clad', fill=clad,
+                               region=+fuel_or & -clad_or)
+    gt_outer_cell = openmc.Cell(name='guide tube outer water', fill=hot_water,
+                                region=+clad_or)
+    guide_tube_universe.add_cells([gt_inner_cell, gt_clad_cell, gt_outer_cell])
+
+    # Create fuel assembly Lattice
+    assembly = openmc.RectLattice(name='Fuel Assembly')
+    assembly.pitch = (pitch/17, pitch/17)
+    assembly.lower_left = (-pitch/2, -pitch/2)
+
+    # Create array indices for guide tube locations in lattice
+    template_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8,
+                           11, 14, 2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
+    template_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8,
+                           8, 11, 11, 11, 11, 11, 13, 13, 14, 14, 14])
+
+    # Create 17x17 array of universes
+    assembly.universes = np.tile(fuel_pin_universe, (17, 17))
+    assembly.universes[template_x, template_y] = guide_tube_universe
+
+    # Create root Cell
+    root_cell = openmc.Cell(name='root cell', fill=assembly)
+    root_cell.region = +min_x & -max_x & +min_y & -max_y
+
+    # Create root Universe
+    model.geometry.root_universe = openmc.Universe(name='root universe')
+    model.geometry.root_universe.add_cell(root_cell)
+
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    model.settings.particles = 100
+    model.settings.source = openmc.Source(space=openmc.stats.Box(
+        [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
+
+    plot = openmc.Plot()
+    plot.basis = 'yz'
+    plot.origin = (0, 0, 40000)
+    plot.width = (25, 250)
+    plot.pixels = (300, 3000)
+    plot.color_by = 'material'
+    model.plots.append(plot)
+
+    return model
+
+
+
+def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5'):
+    """Create a 1D slab model.
+
+    Parameters
+    ----------
+    num_regions : int, optional
+        Number of regions in the problem, each with a unique MGXS dataset.
+        Defaults to 1.
+
+    mat_names : Iterable of str, optional
+        List of the material names to use; defaults to ['mat_1', 'mat_2',...].
+
+    mgxslib_name : str, optional
+        MGXS Library file to use; defaults to '2g.h5'.
+
+    Returns
+    -------
+    model : openmc.model.Model
+        One-group, 1D slab model
+
+    """
+
+    openmc.check_type('num_regions', num_regions, Integral)
+    openmc.check_greater_than('num_regions', num_regions, 0)
+    if mat_names is not None:
+        openmc.check_length('mat_names', mat_names, num_regions)
+        openmc.check_iterable_type('mat_names', mat_names, str)
+    else:
+        mat_names = []
+        for i in range(num_regions):
+            mat_names.append('mat_' + str(i + 1))
+
+    # # Make Materials
+    materials_file = openmc.Materials()
+    macros = []
+    mats = []
+    for i in range(len(mat_names)):
+        macros.append(openmc.Macroscopic('mat_' + str(i + 1)))
+        mats.append(openmc.Material(name=mat_names[i]))
+        mats[-1].set_density('macro', 1.0)
+        mats[-1].add_macroscopic(macros[-1])
+
+    materials_file += mats
+
+    materials_file.cross_sections = mgxslib_name
+
+    # # Make Geometry
+    rad_outer = 929.45
+    # Set a cell boundary to exist for every material above (exclude the 0)
+    rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:]
+
+    # Instantiate Universe
+    root = openmc.Universe(universe_id=0, name='root universe')
+    cells = []
+
+    surfs = []
+    surfs.append(openmc.XPlane(x0=0., boundary_type='reflective'))
+    for r, rad in enumerate(rads):
+        if r == len(rads) - 1:
+            surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum'))
+        else:
+            surfs.append(openmc.XPlane(x0=rad))
+
+    # Instantiate Cells
+    cells = []
+    for c in range(len(surfs) - 1):
+        cells.append(openmc.Cell())
+        cells[-1].region = (+surfs[c] & -surfs[c + 1])
+        cells[-1].fill = mats[c]
+
+    # Register Cells with Universe
+    root.add_cells(cells)
+
+    # Instantiate a Geometry, register the root Universe, and export to XML
+    geometry_file = openmc.Geometry(root)
+
+    # # Make Settings
+    # Instantiate a Settings object, set all runtime parameters
+    settings_file = openmc.Settings()
+    settings_file.energy_mode = 'multi-group'
+    settings_file.tabular_legendre = {'enable': False}
+    settings_file.batches = 10
+    settings_file.inactive = 5
+    settings_file.particles = 1000
+
+    # Build source distribution
+    INF = 1000.
+    bounds = [0., -INF, -INF, rads[0], INF, INF]
+    uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
+    settings_file.source = openmc.Source(space=uniform_dist)
+
+    settings_file.output = {'summary': False}
+
+    model = openmc.model.Model()
+    model.geometry = geometry_file
+    model.materials = materials_file
+    model.settings = settings_file
+    model.xs_data = macros
+
+    return model
+
+if __name__ == "__main__":
+    rtn = pwr_core()
+    rtn.export_to_xml()
+    openmc.plot_geometry()
\ No newline at end of file