diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py
index b68b1032a..f948666cc 100644
--- a/openmc/opencg_compatible.py
+++ b/openmc/opencg_compatible.py
@@ -975,6 +975,19 @@ def get_openmc_lattice(opencg_lattice):
 
 
 def get_opencg_geometry(openmc_geometry):
+    """Return an OpenCG geometry corresponding to an OpenMC geometry.
+
+    Parameters
+    ----------
+    openmc_geometry : openmc.universe.Geometry
+        OpenMC geometry
+
+    Returns
+    -------
+    opencg_geometry : opencg.Geometry
+        Equivalent OpenCG geometry
+
+    """
 
     if not isinstance(openmc_geometry, openmc.Geometry):
         msg = 'Unable to get OpenCG geometry from {0} which is ' \
diff --git a/openmc/statepoint.py b/openmc/statepoint.py
index 4655186b0..89afef7ec 100644
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@@ -522,7 +522,8 @@ class StatePoint(object):
                 sum_sq = np.reshape(sum_sq, new_shape)
 
                 # Set the data for this Tally
-                tally.set_results(sum=sum, sum_sq=sum_sq)
+                tally.sum = sum
+                tally.sum_sq = sum_sq
 
         # Indicate that Tally results have been read
         self._results = True
diff --git a/openmc/summary.py b/openmc/summary.py
index 61cce8b5c..1a5d03df0 100644
--- a/openmc/summary.py
+++ b/openmc/summary.py
@@ -529,6 +529,11 @@ class Summary(object):
             self.tallies[tally_id] = tally
 
     def make_opencg_geometry(self):
+        """Create OpenCG geometry based on the information contained in the summary
+        file. The geometry is stored as the 'opencg_geometry' attribute.
+
+        """
+
         try:
             from openmc.opencg_compatible import get_opencg_geometry
         except ImportError:
@@ -540,6 +545,22 @@ class Summary(object):
             self.opencg_geometry = get_opencg_geometry(self.openmc_geometry)
 
     def get_nuclide_by_zaid(self, zaid):
+        """Return a Nuclide object given the 'zaid' identifier for the nuclide.
+
+        Parameters
+        ----------
+        zaid : int
+            1000*Z + A, where Z is the atomic number of the nuclide and A is the
+            mass number. For example, the zaid for U-235 is 92235.
+
+        Returns
+        -------
+        nuclide : openmc.nuclide.Nuclide
+            Nuclide matching the specified zaid
+
+        """
+
+
         for index, nuclide in self.nuclides.items():
             if nuclide._zaid == zaid:
                 return nuclide
@@ -547,6 +568,20 @@ class Summary(object):
         return None
 
     def get_material_by_id(self, material_id):
+        """Return a Material object given the material id
+
+        Parameters
+        ----------
+        id : int
+            Unique identifier for the material
+
+        Returns
+        -------
+        material : openmc.material.Material
+            Material with given id
+
+        """
+
         for index, material in self.materials.items():
             if material._id == material_id:
                 return material
@@ -554,6 +589,20 @@ class Summary(object):
         return None
 
     def get_surface_by_id(self, surface_id):
+        """Return a Surface object given the surface id
+
+        Parameters
+        ----------
+        id : int
+            Unique identifier for the surface
+
+        Returns
+        -------
+        surface : openmc.surface.Surface
+            Surface with given id
+
+        """
+
         for index, surface in self.surfaces.items():
             if surface._id == surface_id:
                 return surface
@@ -561,6 +610,20 @@ class Summary(object):
         return None
 
     def get_cell_by_id(self, cell_id):
+        """Return a Cell object given the cell id
+
+        Parameters
+        ----------
+        id : int
+            Unique identifier for the cell
+
+        Returns
+        -------
+        cell : openmc.universe.Cell
+            Cell with given id
+
+        """
+
         for index, cell in self.cells.items():
             if cell._id == cell_id:
                 return cell
@@ -568,6 +631,20 @@ class Summary(object):
         return None
 
     def get_universe_by_id(self, universe_id):
+        """Return a Universe object given the universe id
+
+        Parameters
+        ----------
+        id : int
+            Unique identifier for the universe
+
+        Returns
+        -------
+        universe : openmc.universe.Universe
+            Universe with given id
+
+        """
+
         for index, universe in self.universes.items():
             if universe._id == universe_id:
                 return universe
@@ -575,6 +652,20 @@ class Summary(object):
         return None
 
     def get_lattice_by_id(self, lattice_id):
+        """Return a Lattice object given the lattice id
+
+        Parameters
+        ----------
+        id : int
+            Unique identifier for the lattice
+
+        Returns
+        -------
+        lattice : openmc.universe.Lattice
+            Lattice with given id
+
+        """
+
         for index, lattice in self.lattices.items():
             if lattice._id == lattice_id:
                 return lattice
diff --git a/openmc/surface.py b/openmc/surface.py
index 0e958b456..0f10dae38 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -7,16 +7,44 @@ from openmc.constants import BC_TYPES
 # A static variable for auto-generated Surface IDs
 AUTO_SURFACE_ID = 10000
 
+
 def reset_auto_surface_id():
     global AUTO_SURFACE_ID
     AUTO_SURFACE_ID = 10000
 
 
-
 class Surface(object):
+    """A two-dimensional surface that can be used define regions of space with an
+    associated boundary condition.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    name : str
+        Name of the surface
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier for the surface
+    name : str
+        Name of the surface
+    type : str
+        Type of the surface, e.g. 'x-plane'
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface.
+    coeffs : dict
+        Dictionary of surface coefficients
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission', name=''):
-
         # Initialize class attributes
         self.id = surface_id
         self.name = name
@@ -32,35 +60,28 @@ class Surface(object):
         # proper order
         self._coeff_keys = []
 
-
     @property
     def id(self):
         return self._id
 
-
     @property
     def name(self):
         return self._name
 
-
     @property
     def type(self):
         return self._type
 
-
     @property
     def boundary_type(self):
         return self._boundary_type
 
-
     @property
     def coeffs(self):
         return self._coeffs
 
-
     @id.setter
     def id(self, surface_id):
-
         if surface_id is None:
             global AUTO_SURFACE_ID
             self._id = AUTO_SURFACE_ID
@@ -80,10 +101,8 @@ class Surface(object):
         else:
             self._id = surface_id
 
-
     @name.setter
     def name(self, name):
-
         if not is_string(name):
             msg = 'Unable to set name for Surface ID={0} with a non-string ' \
                   'value {1}'.format(self._id, name)
@@ -92,16 +111,14 @@ class Surface(object):
         else:
             self._name = name
 
-
     @boundary_type.setter
     def boundary_type(self, boundary_type):
-
         if not is_string(boundary_type):
             msg = 'Unable to set boundary type for Surface ID={0} with a ' \
                   'non-string value {1}'.format(self._id, boundary_type)
             raise ValueError(msg)
 
-        elif not boundary_type in BC_TYPES.values():
+        elif boundary_type not in BC_TYPES.values():
             msg = 'Unable to set boundary type for Surface ID={0} to ' \
                   '{1} which is not trasmission, vacuum or ' \
                   'reflective'.format(boundary_type)
@@ -110,9 +127,7 @@ class Surface(object):
         else:
             self._boundary_type = boundary_type
 
-
     def __repr__(self):
-
         string = 'Surface\n'
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
         string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@@ -128,9 +143,7 @@ class Surface(object):
 
         return string
 
-
     def create_xml_subelement(self):
-
         element = ET.Element("surface")
         element.set("id", str(self._id))
 
@@ -139,65 +152,84 @@ class Surface(object):
 
         element.set("type", self._type)
         element.set("boundary", self._boundary_type)
-
-        coeffs = ''
-
-        for coeff in self._coeff_keys:
-            coeffs += '{0} '.format(self._coeffs[coeff])
-
-        element.set("coeffs", coeffs.rstrip(' '))
+        element.set("coeffs", ' '.join(map(str, self._coeffs)))
 
         return element
 
 
-
 class Plane(Surface):
+    """An arbitrary plane of the form :math:`Ax + By + Cz = D`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    A : float
+        The 'A' parameter for the surface
+    B : float
+        The 'B' parameter for the surface
+    C : float
+        The 'C' parameter for the surface
+    D : float
+        The 'D' parameter for the surface
+    name : str
+        Name of the plane
+
+    Attributes
+    ----------
+    a : float
+        The 'A' parameter for the surface
+    b : float
+        The 'B' parameter for the surface
+    c : float
+        The 'C' parameter for the surface
+    d : float
+        The 'D' parameter for the surface
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
-                 A=None, B=None, C=None, D=None, name='',):
-
+                 A=None, B=None, C=None, D=None, name=''):
         # Initialize Plane class attributes
         super(Plane, self).__init__(surface_id, boundary_type, name=name)
 
         self._type = 'plane'
         self._coeff_keys = ['A', 'B', 'C', 'D']
 
-        if not A is None:
+        if A is not None:
             self.a = A
 
-        if not B is None:
+        if B is not None:
             self.b = B
 
-        if not C is None:
+        if C is not None:
             self.c = C
 
-        if not D is None:
+        if D is not None:
             self.d = D
 
-
     @property
     def a(self):
         return self.coeffs['A']
 
-
     @property
     def b(self):
         return self.coeffs['B']
 
-
     @property
     def c(self):
         return self.coeffs['C']
 
-
     @property
     def d(self):
         return self.coeffs['D']
 
-
     @a.setter
     def a(self, A):
-
         if not is_integer(A) and not is_float(A):
             msg = 'Unable to set A coefficient for Plane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, A)
@@ -205,10 +237,8 @@ class Plane(Surface):
 
         self._coeffs['A'] = A
 
-
     @b.setter
     def b(self, B):
-
         if not is_integer(B) and not is_float(B):
             msg = 'Unable to set B coefficient for Plane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, B)
@@ -216,10 +246,8 @@ class Plane(Surface):
 
         self._coeffs['B'] = B
 
-
     @c.setter
     def c(self, C):
-
         if not is_integer(C) and not is_float(C):
             msg = 'Unable to set C coefficient for Plane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, C)
@@ -227,10 +255,8 @@ class Plane(Surface):
 
         self._coeffs['C'] = C
 
-
     @d.setter
     def d(self, D):
-
         if not is_integer(D) and not is_float(D):
             msg = 'Unable to set D coefficient for Plane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, D)
@@ -239,30 +265,47 @@ class Plane(Surface):
         self._coeffs['D'] = D
 
 
-
 class XPlane(Plane):
+    """A plane perpendicular to the x axis, i.e. a surface of the form :math:`x -
+    x_0 = 0`
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        Location of the plane
+    name : str
+        Name of the plane
+
+    Attributes
+    ----------
+    x0 : float
+        Location of the plane
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, name=''):
-
         # Initialize XPlane class attributes
         super(XPlane, self).__init__(surface_id, boundary_type, name=name)
 
         self._type = 'x-plane'
         self._coeff_keys = ['x0']
 
-        if not x0 is None:
+        if x0 is not None:
             self.x0 = x0
 
-
     @property
     def x0(self):
         return self.coeff['x0']
 
-
     @x0.setter
     def x0(self, x0):
-
         if not is_integer(x0) and not is_float(x0):
             msg = 'Unable to set x0 coefficient for XPlane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, x0)
@@ -271,30 +314,47 @@ class XPlane(Plane):
         self._coeffs['x0'] = x0
 
 
-
 class YPlane(Plane):
+    """A plane perpendicular to the y axis, i.e. a surface of the form :math:`y -
+    y_0 = 0`
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    y0 : float
+        Location of the plane
+    name : str
+        Name of the plane
+
+    Attributes
+    ----------
+    y0 : float
+        Location of the plane
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  y0=None, name=''):
-
         # Initialize YPlane class attributes
         super(YPlane, self).__init__(surface_id, boundary_type, name=name)
 
         self._type = 'y-plane'
         self._coeff_keys = ['y0']
 
-        if not y0 is None:
+        if y0 is not None:
             self.y0 = y0
 
-
     @property
     def y0(self):
         return self.coeffs['y0']
 
-
     @y0.setter
     def y0(self, y0):
-
         if not is_integer(y0) and not is_float(y0):
             msg = 'Unable to set y0 coefficient for XPlane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, y0)
@@ -303,30 +363,47 @@ class YPlane(Plane):
         self._coeffs['y0'] = y0
 
 
-
 class ZPlane(Plane):
+    """A plane perpendicular to the z axis, i.e. a surface of the form :math:`z -
+    z_0 = 0`
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    z0 : float
+        Location of the plane
+    name : str
+        Name of the plane
+
+    Attributes
+    ----------
+    z0 : float
+        Location of the plane
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  z0=None, name=''):
-
         # Initialize ZPlane class attributes
         super(ZPlane, self).__init__(surface_id, boundary_type, name=name)
 
         self._type = 'z-plane'
         self._coeff_keys = ['z0']
 
-        if not z0 is None:
+        if z0 is not None:
             self.z0 = z0
 
-
     @property
     def z0(self):
         return self.coeffs['z0']
 
-
     @z0.setter
     def z0(self, z0):
-
         if not is_integer(z0) and not is_float(z0):
             msg = 'Unable to set z0 coefficient for ZPlane ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, z0)
@@ -335,29 +412,45 @@ class ZPlane(Plane):
         self._coeffs['z0'] = z0
 
 
-
 class Cylinder(Surface):
+    """A cylinder whose length is parallel to the x-, y-, or z-axis.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    R : float
+        Radius of the cylinder
+    name : str
+        Name of the cylinder
+
+    Attributes
+    ----------
+    r : float
+        Radius of the cylinder
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  R=None, name=''):
-
         # Initialize Cylinder class attributes
         super(Cylinder, self).__init__(surface_id, boundary_type, name=name)
 
         self._coeff_keys = ['R']
 
-        if not R is None:
+        if R is not None:
             self.r = R
 
-
     @property
     def r(self):
         return self.coeffs['R']
 
-
     @r.setter
     def r(self, R):
-
         if not is_integer(R) and not is_float(R):
             msg = 'Unable to set R coefficient for Cylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, R)
@@ -366,38 +459,60 @@ class Cylinder(Surface):
         self._coeffs['R'] = R
 
 
-
 class XCylinder(Cylinder):
+    """An infinite cylinder whose length is parallel to the x-axis. This is a
+    quadratic surface of the form :math:`(y - y_0)^2 + (z - z_0)^2 = R^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    y0 : float
+        y-coordinate of the center of the cylinder
+    z0 : float
+        z-coordinate of the center of the cylinder
+    R : float
+        Radius of the cylinder
+    name : str
+        Name of the cylinder
+
+    Attributes
+    ----------
+    y0 : float
+        y-coordinate of the center of the cylinder
+    z0 : float
+        z-coordinate of the center of the cylinder
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  y0=None, z0=None, R=None, name=''):
-
         # Initialize XCylinder class attributes
         super(XCylinder, self).__init__(surface_id, boundary_type, R, name=name)
 
         self._type = 'x-cylinder'
         self._coeff_keys = ['y0', 'z0', 'R']
 
-        if not y0 is None:
+        if y0 is not None:
             self.y0 = y0
 
-        if not z0 is None:
+        if z0 is not None:
             self.z0 = z0
 
-
     @property
     def y0(self):
         return self.coeffs['y0']
 
-
     @property
     def z0(self):
         return self.coeffs['z0']
 
-
     @y0.setter
     def y0(self, y0):
-
         if not is_integer(y0) and not is_float(y0):
             msg = 'Unable to set y0 coefficient for XCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, y0)
@@ -405,10 +520,8 @@ class XCylinder(Cylinder):
 
         self._coeffs['y0'] = y0
 
-
     @z0.setter
     def z0(self, z0):
-
         if not is_integer(z0) and not is_float(z0):
             msg = 'Unable to set z0 coefficient for XCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, z0)
@@ -417,38 +530,60 @@ class XCylinder(Cylinder):
         self._coeffs['z0'] = z0
 
 
-
 class YCylinder(Cylinder):
+    """An infinite cylinder whose length is parallel to the y-axis. This is a
+    quadratic surface of the form :math:`(x - x_0)^2 + (z - z_0)^2 = R^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the center of the cylinder
+    z0 : float
+        z-coordinate of the center of the cylinder
+    R : float
+        Radius of the cylinder
+    name : str
+        Name of the cylinder
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the center of the cylinder
+    z0 : float
+        z-coordinate of the center of the cylinder
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, z0=None, R=None, name=''):
-
         # Initialize YCylinder class attributes
         super(YCylinder, self).__init__(surface_id, boundary_type, R, name=name)
 
         self._type = 'y-cylinder'
         self._coeff_keys = ['x0', 'z0', 'R']
 
-        if not x0 is None:
+        if x0 is not None:
             self.x0 = x0
 
-        if not z0 is None:
+        if z0 is not None:
             self.z0 = z0
 
-
     @property
     def x0(self):
         return self.coeffs['x0']
 
-
     @property
     def z0(self):
         return self.coeffs['z0']
 
-
     @x0.setter
     def x0(self, x0):
-
         if not is_integer(x0) and not is_float(x0):
             msg = 'Unable to set x0 coefficient for YCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, x0)
@@ -456,10 +591,8 @@ class YCylinder(Cylinder):
 
         self._coeffs['x0'] = x0
 
-
     @z0.setter
     def z0(self, z0):
-
         if not is_integer(z0) and not is_float(z0):
             msg = 'Unable to set z0 coefficient for YCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, z0)
@@ -469,37 +602,59 @@ class YCylinder(Cylinder):
 
 
 class ZCylinder(Cylinder):
+    """An infinite cylinder whose length is parallel to the z-axis. This is a
+    quadratic surface of the form :math:`(x - x_0)^2 + (y - y_0)^2 = R^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the center of the cylinder
+    y0 : float
+        y-coordinate of the center of the cylinder
+    R : float
+        Radius of the cylinder
+    name : str
+        Name of the cylinder
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the center of the cylinder
+    y0 : float
+        y-coordinate of the center of the cylinder
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, R=None, name=''):
-
         # Initialize ZCylinder class attributes
-        # Initialize YPlane class attributes
         super(ZCylinder, self).__init__(surface_id, boundary_type, R, name=name)
 
         self._type = 'z-cylinder'
         self._coeff_keys = ['x0', 'y0', 'R']
 
-        if not x0 is None:
+        if x0 is not None:
             self.x0 = x0
 
-        if not y0 is None:
+        if y0 is not None:
             self.y0 = y0
 
-
     @property
     def x0(self):
         return self.coeffs['x0']
 
-
     @property
     def y0(self):
         return self.coeffs['y0']
 
-
     @x0.setter
     def x0(self, x0):
-
         if not is_integer(x0) and not is_float(x0):
             msg = 'Unable to set x0 coefficient for ZCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, x0)
@@ -507,10 +662,8 @@ class ZCylinder(Cylinder):
 
         self._coeffs['x0'] = x0
 
-
     @y0.setter
     def y0(self, y0):
-
         if not is_integer(y0) and not is_float(y0):
             msg = 'Unable to set y0 coefficient for ZCylinder ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, y0)
@@ -519,54 +672,79 @@ class ZCylinder(Cylinder):
         self._coeffs['y0'] = y0
 
 
-
 class Sphere(Surface):
+    """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the center of the sphere
+    y0 : float
+        y-coordinate of the center of the sphere
+    z0 : float
+        z-coordinate of the center of the sphere
+    R : float
+        Radius of the sphere
+    name : str
+        Name of the sphere
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the center of the sphere
+    y0 : float
+        y-coordinate of the center of the sphere
+    z0 : float
+        z-coordinate of the center of the sphere
+    R : float
+        Radius of the sphere
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, z0=None, R=None, name=''):
-
         # Initialize Sphere class attributes
         super(Sphere, self).__init__(surface_id, boundary_type, name=name)
 
         self._type = 'sphere'
         self._coeff_keys = ['x0', 'y0', 'z0', 'R']
 
-        if not x0 is None:
+        if x0 is not None:
             self.x0 = x0
 
-        if not y0 is None:
+        if y0 is not None:
             self.y0 = y0
 
-        if not z0 is None:
+        if z0 is not None:
             self.z0 = z0
 
-        if not R is None:
+        if R is not None:
             self.r = R
 
-
     @property
     def x0(self):
         return self.coeffs['x0']
 
-
     @property
     def y0(self):
         return self.coeffs['y0']
 
-
     @property
     def z0(self):
         return self.coeffs['z0']
 
-
     @property
     def r(self):
         return self.coeffs['R']
 
-
     @x0.setter
     def x0(self, x0):
-
         if not is_integer(x0) and not is_float(x0):
             msg = 'Unable to set x0 coefficient for Sphere ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, x0)
@@ -574,10 +752,8 @@ class Sphere(Surface):
 
         self._coeffs['x0'] = x0
 
-
     @y0.setter
     def y0(self, y0):
-
         if not is_integer(y0) and not is_float(y0):
             msg = 'Unable to set y0 coefficient for Sphere ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, y0)
@@ -585,10 +761,8 @@ class Sphere(Surface):
 
         self._coeffs['y0'] = y0
 
-
     @z0.setter
     def z0(self, z0):
-
         if not is_integer(z0) and not is_float(z0):
             msg = 'Unable to set z0 coefficient for Sphere ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, z0)
@@ -596,10 +770,8 @@ class Sphere(Surface):
 
         self._coeffs['z0'] = z0
 
-
     @r.setter
     def r(self, R):
-
         if not is_integer(R) and not is_float(R):
             msg = 'Unable to set R coefficient for Sphere ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, R)
@@ -608,53 +780,78 @@ class Sphere(Surface):
         self._coeffs['R'] = R
 
 
-
 class Cone(Surface):
+    """A conical surface parallel to the x-, y-, or z-axis.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+    name : str
+        Name of the cone
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, z0=None, R2=None, name=''):
-
         # Initialize Cone class attributes
         super(Cone, self).__init__(surface_id, boundary_type, name=name)
 
         self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
 
-        if not x0 is None:
+        if x0 is not None:
             self.x0 = x0
 
-        if not y0 is None:
+        if y0 is not None:
             self.y0 = y0
 
-        if not z0 is None:
+        if z0 is not None:
             self.z0 = z0
 
-        if not R2 is None:
+        if R2 is not None:
             self.r2 = R2
 
-
     @property
     def x0(self):
         return self.coeffs['x0']
 
-
     @property
     def y0(self):
         return self.coeffs['y0']
 
-
     @property
     def z0(self):
         return self.coeffs['z0']
 
-
     @property
     def r2(self):
         return self.coeffs['r2']
 
-
     @x0.setter
     def x0(self, x0):
-
         if not is_integer(x0) and not is_float(x0):
             msg = 'Unable to set x0 coefficient for Cone ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, x0)
@@ -662,10 +859,8 @@ class Cone(Surface):
 
         self._coeffs['x0'] = x0
 
-
     @y0.setter
     def y0(self, y0):
-
         if not is_integer(y0) and not is_float(y0):
             msg = 'Unable to set y0 coefficient for Cone ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, y0)
@@ -673,10 +868,8 @@ class Cone(Surface):
 
         self._coeffs['y0'] = y0
 
-
     @z0.setter
     def z0(self, z0):
-
         if not is_integer(z0) and not is_float(z0):
             msg = 'Unable to set z0 coefficient for Cone ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, z0)
@@ -684,10 +877,8 @@ class Cone(Surface):
 
         self._coeffs['z0'] = z0
 
-
     @r2.setter
     def r2(self, R2):
-
         if not is_integer(R2) and not is_float(R2):
             msg = 'Unable to set R^2 coefficient for Cone ID={0} to a ' \
                   'non-integer value {1}'.format(self._id, R2)
@@ -696,12 +887,44 @@ class Cone(Surface):
         self._coeffs['R2'] = R2
 
 
-
 class XCone(Cone):
+    """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
+    R^2 (x - x_0)^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+    name : str
+        Name of the cone
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, z0=None, R2=None, name=''):
-
         # Initialize XCone class attributes
         super(XCone, self).__init__(surface_id, boundary_type, x0, y0,
                                     z0, R2, name=name)
@@ -709,12 +932,44 @@ class XCone(Cone):
         self._type = 'x-cone'
 
 
-
 class YCone(Cone):
+    """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 =
+    R^2 (y - y_0)^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+    name : str
+        Name of the cone
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, z0=None, R2=None, name=''):
-
         # Initialize YCone class attributes
         super(YCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
                                     R2, name=name)
@@ -722,12 +977,44 @@ class YCone(Cone):
         self._type = 'y-cone'
 
 
-
 class ZCone(Cone):
+    """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 =
+    R^2 (z - z_0)^2`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+    name : str
+        Name of the cone
+
+    Attributes
+    ----------
+    x0 : float
+        x-coordinate of the apex
+    y0 : float
+        y-coordinate of the apex
+    z0 : float
+        z-coordinate of the apex
+    R2 : float
+        Parameter related to the aperature
+
+    """
 
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=None, y0=None, z0=None, R2=None, name=''):
-
         # Initialize ZCone class attributes
         super(ZCone, self).__init__(surface_id, boundary_type, x0, y0, z0,
                                     R2, name=name)
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 512e93c9a..f87a1a79c 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -1,5 +1,6 @@
 import copy
 import os
+import pickle
 import itertools
 from xml.etree import ElementTree as ET
 
@@ -21,9 +22,58 @@ def reset_auto_tally_id():
 
 
 class Tally(object):
+    """A tally defined by a set of scores that are accumulated for a list of
+    nuclides given a set of filters.
+
+    Parameters
+    ----------
+    tally_id : int
+        Unique identifier for the tally
+    name : str
+        Name of the taly
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier for the tally
+    name : str
+        Name of the tally
+    filters : list of openmc.filter.Filter
+        List of specified filters for the tally
+    nuclides : list of openmc.nuclide.Nuclide
+        List of nuclides to score results for
+    scores : list of str
+        List of defined scores, e.g. 'flux', 'fission', etc.
+    estimator : {'analog', 'tracklength'}
+        Type of estimator for the tally
+    triggers : list of openmc.trigger.Trigger
+        List of tally triggers
+    num_score_bins : int
+        Total number of scores, accounting for the fact that a single
+        user-specified score might have multiple bins
+    num_scores : int
+        Total number of user-specified scores
+    num_filter_bins : int
+        Total number of filter bins accounting for all filters
+    num_bins : int
+        Total number of bins for the tally
+    num_realizations : int
+        Total number of realizations
+    with_summary : bool
+        Whether or not a Summary has been linked
+    sum : ndarray
+        An array containing the sum of each independent realization for each bin
+    sum_sq : ndarray
+        An array containing the sum of each independent realization squared for
+        each bin
+    mean : ndarray
+        An array containing the sample mean for each bin
+    std_dev : ndarray
+        An array containing the sample standard deviation for each bin
+
+    """
 
     def __init__(self, tally_id=None, name=''):
-
         # Initialize Tally class attributes
         self.id = tally_id
         self.name = name
@@ -42,14 +92,11 @@ class Tally(object):
         self._mean = None
         self._std_dev = None
 
-
     def __deepcopy__(self, memo):
-
         existing = memo.get(id(self))
 
         # If this is the first time we have tried to copy this object, create a copy
         if existing is None:
-
             clone = type(self).__new__(type(self))
             clone.id = self.id
             clone.name = self.name
@@ -63,19 +110,19 @@ class Tally(object):
 
             clone.filters = []
             for filter in self.filters:
-              clone.add_filter(copy.deepcopy(filter, memo))
+                clone.add_filter(copy.deepcopy(filter, memo))
 
             clone.nuclides = []
             for nuclide in self.nuclides:
-              clone.add_nuclide(copy.deepcopy(nuclide, memo))
+                clone.add_nuclide(copy.deepcopy(nuclide, memo))
 
             clone.scores = []
             for score in self.scores:
-              clone.add_score(score)
+                clone.add_score(score)
 
             clone.triggers = []
             for trigger in self.triggers:
-              clone.add_trigger(trigger)
+                clone.add_trigger(trigger)
 
             memo[id(self)] = clone
 
@@ -85,9 +132,7 @@ class Tally(object):
         else:
             return existing
 
-
     def __eq__(self, tally2):
-
         # Check all filters
         if len(self.filters) != len(tally2.filters):
             return False
@@ -117,7 +162,6 @@ class Tally(object):
 
         return True
 
-
     def __hash__(self):
         hashable = []
 
@@ -135,9 +179,7 @@ class Tally(object):
 
         return hash(tuple(hashable))
 
-
     def __add__(self, other):
-
         # FIXME: Error checking: must check that results has been
         # set and that # bins is the same
 
@@ -145,55 +187,40 @@ class Tally(object):
         new_tally._mean = self._mean + other._mean
         new_tally._std_dev = np.sqrt(self.std_dev**2 + other.std_dev**2)
 
-
     @property
     def id(self):
         return self._id
 
-
     @property
     def name(self):
         return self._name
 
-
     @property
     def filters(self):
         return self._filters
 
-
-    @property
-    def num_filters(self):
-        return len(self._filters)
-
-
     @property
     def nuclides(self):
         return self._nuclides
 
-
     @property
     def num_nuclides(self):
         return len(self._nuclides)
 
-
     @property
     def scores(self):
         return self._scores
 
-
     @property
     def num_scores(self):
         return len(self._scores)
 
-
     @property
     def num_score_bins(self):
         return self._num_score_bins
 
-
     @property
     def num_filter_bins(self):
-
         num_bins = 1
 
         for filter in self._filters:
@@ -201,7 +228,6 @@ class Tally(object):
 
         return num_bins
 
-
     @property
     def num_bins(self):
         num_bins = self.num_filter_bins
@@ -209,50 +235,40 @@ class Tally(object):
         num_bins *= self.num_score_bins
         return num_bins
 
-
     @property
     def estimator(self):
         return self._estimator
 
-
     @property
     def triggers(self):
         return self._triggers
 
-
     @property
     def num_realizations(self):
         return self._num_realizations
 
-
     @property
     def with_summary(self):
         return self._with_summary
 
-
     @property
     def sum(self):
         return self._sum
 
-
     @property
     def sum_sq(self):
         return self._sum_sq
 
-
     @property
     def mean(self):
         return self._mean
 
-
     @property
     def std_dev(self):
         return self._std_dev
 
-
     @estimator.setter
     def estimator(self, estimator):
-
         if not estimator in ['analog', 'tracklength']:
             msg = 'Unable to set the estimator for Tally ID={0} to {1} since ' \
                   'it is not a valid estimator type'.format(self.id, estimator)
@@ -260,8 +276,15 @@ class Tally(object):
 
         self._estimator = estimator
 
-
     def add_trigger(self, trigger):
+        """Add a tally trigger to the tally
+
+        Parameters
+        ----------
+        trigger : openmc.trigger.Trigger
+            Trigger to add
+
+        """
 
         if not isinstance(trigger, Trigger):
             msg = 'Unable to add a tally trigger for Tally ID={0} to ' \
@@ -270,10 +293,8 @@ class Tally(object):
 
         self._triggers.append(trigger)
 
-
     @id.setter
     def id(self, tally_id):
-
         if tally_id is None:
             global AUTO_TALLY_ID
             self._id = AUTO_TALLY_ID
@@ -292,10 +313,8 @@ class Tally(object):
         else:
             self._id = tally_id
 
-
     @name.setter
     def name(self, name):
-
         if not is_string(name):
             msg = 'Unable to set name for Tally ID={0} with a non-string ' \
                   'value "{1}"'.format(self.id, name)
@@ -304,8 +323,15 @@ class Tally(object):
         else:
             self._name = name
 
-
     def add_filter(self, filter):
+        """Add a filter to the tally
+
+        Parameters
+        ----------
+        filter : openmc.filter.Filter
+            Filter to add
+
+        """
 
         global filters
 
@@ -316,12 +342,27 @@ class Tally(object):
 
         self._filters.append(filter)
 
-
     def add_nuclide(self, nuclide):
+        """Specify that scores for a particular nuclide should be accumulated
+
+        Parameters
+        ----------
+        nuclide : openmc.nuclide.Nuclide
+            Nuclide to add
+
+        """
+
         self._nuclides.append(nuclide)
 
-
     def add_score(self, score):
+        """Specify a quantity to be scored
+
+        Parameters
+        ----------
+        score : str
+            Score to be accumulated, e.g. 'flux'
+
+        """
 
         if not is_string(score):
             msg = 'Unable to add score "{0}" to Tally ID={1} since it is ' \
@@ -334,15 +375,12 @@ class Tally(object):
         else:
             self._scores.append(score)
 
-
     @num_score_bins.setter
     def num_score_bins(self, num_score_bins):
         self._num_score_bins = num_score_bins
 
-
     @num_realizations.setter
     def num_realizations(self, num_realizations):
-
         if not is_integer(num_realizations):
             msg = 'Unable to set the number of realizations to "{0}" for ' \
                   'Tally ID={1} since it is not an ' \
@@ -357,10 +395,8 @@ class Tally(object):
 
         self._num_realizations = num_realizations
 
-
     @with_summary.setter
     def with_summary(self, with_summary):
-
         if not isinstance(with_summary, bool):
             msg = 'Unable to set with_summary to a non-boolean ' \
                   'value "{0}"'.format(with_summary)
@@ -368,56 +404,85 @@ class Tally(object):
 
         self._with_summary = with_summary
 
-
-    def set_results(self, sum, sum_sq):
-
+    @sum.setter
+    def sum(self, sum):
         if not isinstance(sum, (tuple, list, np.ndarray)):
             msg = 'Unable to set the sum to "{0}" for Tally ID={1} since ' \
                   'it is not a Python tuple/list or NumPy ' \
                   'array'.format(sum, self.id)
             raise ValueError(msg)
+        self._sum = sum
 
+    @sum_sq.setter
+    def sum_sq(self, sum_sq):
         if not isinstance(sum_sq, (tuple, list, np.ndarray)):
             msg = 'Unable to set the sum to "{0}" for Tally ID={1} since ' \
                   'it is not a Python tuple/list or NumPy ' \
                   'array'.format(sum_sq, self.id)
             raise ValueError(msg)
-
-        self._sum = sum
         self._sum_sq = sum_sq
 
-
     def remove_score(self, score):
+        """Remove a score from the tally
 
-        if not score in self.scores:
+        Parameters
+        ----------
+        score : str
+            Score to remove
+
+        """
+
+        if score not in self.scores:
             msg = 'Unable to remove score "{0}" from Tally ID={1} since the ' \
                   'Tally does not contain this score'.format(score, self.id)
             ValueError(msg)
 
         self._scores.remove(score)
 
-
     def remove_filter(self, filter):
+        """Remove a filter from the tally
 
-        if not filter in self.filters:
+        Parameters
+        ----------
+        filter : openmc.filter.Filter
+            Filter to remove
+
+        """
+
+        if filter not in self.filters:
             msg = 'Unable to remove filter "{0}" from Tally ID={1} since the ' \
                   'Tally does not contain this filter'.format(filter, self.id)
             ValueError(msg)
 
         self._filters.remove(filter)
 
-
     def remove_nuclide(self, nuclide):
+        """Remove a nuclide from the tally
 
-        if not nuclide in self.nuclides:
+        Parameters
+        ----------
+        nuclide : openmc.nuclide.Nuclide
+            Nuclide to remove
+
+        """
+
+        if nuclide not in self.nuclides:
             msg = 'Unable to remove nuclide "{0}" from Tally ID={1} since the ' \
                   'Tally does not contain this nuclide'.format(nuclide, self.id)
             ValueError(msg)
 
         self._nuclides.remove(nuclide)
 
-
     def compute_std_dev(self, t_value=1.0):
+        """Compute the sample mean and standard deviation for each bin in the tally
+
+        Parameters
+        ----------
+        t_value : float, optional
+            Student's t-value applied to the uncertainty. Defaults to 1.0,
+            meaning the reported value is the sample standard deviation.
+
+        """
 
         # Calculate sample mean and standard deviation
         self._mean = self.sum / self.num_realizations
@@ -425,9 +490,7 @@ class Tally(object):
                                  self.mean**2) / (self.num_realizations - 1))
         self._std_dev *= t_value
 
-
     def __repr__(self):
-
         string = 'Tally\n'
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
         string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name)
@@ -453,8 +516,15 @@ class Tally(object):
 
         return string
 
-
     def can_merge(self, tally):
+        """Determine if another tally can be merged with this one
+
+        Parameters
+        ----------
+        tally : Tally
+            Tally to check for merging
+
+        """
 
         if not isinstance(tally, Tally):
             return False
@@ -491,8 +561,20 @@ class Tally(object):
         # Tallies are mergeable if all conditional checks passed
         return True
 
-
     def merge(self, tally):
+        """Merge another tally with this one
+
+        Parameters
+        ----------
+        tally : Tally
+            Tally to merge with this one
+
+        Returns
+        -------
+        merged_tally : Tally
+            Merged tallies
+
+        """
 
         if not self.can_merge(tally):
             msg = 'Unable to merge tally ID={0} with {1}'.format(tally.id, self.id)
@@ -522,8 +604,15 @@ class Tally(object):
 
         return merged_tally
 
-
     def get_tally_xml(self):
+        """Return XML representation of the tally
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing tally data
+
+        """
 
         element = ET.Element("tally")
 
@@ -536,12 +625,10 @@ class Tally(object):
 
         # Optional Tally filters
         for filter in self.filters:
-
             subelement = ET.SubElement(element, "filter")
             subelement.set("type", str(filter.type))
 
             if not filter.bins is None:
-
                 bins = ''
                 for bin in filter.bins:
                     bins += '{0} '.format(bin)
@@ -550,7 +637,6 @@ class Tally(object):
 
         # Optional Nuclides
         if len(self.nuclides) > 0:
-
             nuclides = ''
             for nuclide in self.nuclides:
                 if isinstance(nuclide, Nuclide):
@@ -568,7 +654,6 @@ class Tally(object):
             raise ValueError(msg)
 
         else:
-
             scores = ''
             for score in self.scores:
                 scores += '{0} '.format(score)
@@ -587,8 +672,20 @@ class Tally(object):
 
         return element
 
-
     def find_filter(self, filter_type):
+        """Return a filter in the tally that matches a specified type
+
+        Parameters
+        ----------
+        filter_type : str
+            Type of the filter, e.g. 'mesh'
+
+        Returns
+        -------
+        filter : openmc.filter.Filter
+            Filter from this tally with matching type
+
+        """
 
         filter = None
 
@@ -606,27 +703,27 @@ class Tally(object):
 
         return filter
 
-
     def get_filter_index(self, filter_type, filter_bin):
         """Returns the index in the Tally's results array for a Filter bin
 
         Parameters
         ----------
         filter_type : str
-             The type of Filter (e.g., 'cell', 'energy', etc.)
+            The type of Filter (e.g., 'cell', 'energy', etc.)
 
         filter_bin : int, list
-              The bin is an integer ID for 'material', 'surface', 'cell',
-              'cellborn', and 'universe' Filters. The bin is an integer for
-              the cell instance ID for 'distribcell' Filters. The bin is
-              a 2-tuple of floats for 'energy' and 'energyout' filters
-              corresponding to the energy boundaries of the bin of interest.
-              The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to
-              the mesh cell of interest.
+            The bin is an integer ID for 'material', 'surface', 'cell',
+            'cellborn', and 'universe' Filters. The bin is an integer for the
+            cell instance ID for 'distribcell' Filters. The bin is a 2-tuple of
+            floats for 'energy' and 'energyout' filters corresponding to the
+            energy boundaries of the bin of interest.  The bin is a (x,y,z)
+            3-tuple for 'mesh' filters corresponding to the mesh cell of
+            interest.
 
         Returns
         -------
              The index in the Tally data array for this filter bin.
+
         """
 
         # Find the equivalent Filter in this Tally's list of Filters
@@ -636,23 +733,25 @@ class Tally(object):
         filter_index = filter.get_bin_index(filter_bin)
         return filter_index
 
-
     def get_nuclide_index(self, nuclide):
         """Returns the index in the Tally's results array for a Nuclide bin
 
         Parameters
         ----------
         nuclide : str
-             The name of the Nuclide (e.g., 'H-1', 'U-238')
+            The name of the Nuclide (e.g., 'H-1', 'U-238')
 
         Returns
         -------
-             The index in the Tally data array for this nuclide.
+        nuclide_index : int
+            The index in the Tally data array for this nuclide.
 
         Raises
         ------
-             KeyError : An error when the argument passed to the 'nuclide'
-             parameter cannot be found in the Tally.
+        KeyError
+            When the argument passed to the 'nuclide' parameter cannot be found
+            in the Tally.
+
         """
 
         nuclide_index = -1
@@ -679,23 +778,25 @@ class Tally(object):
         else:
             return nuclide_index
 
-
     def get_score_index(self, score):
         """Returns the index in the Tally's results array for a score bin
 
         Parameters
         ----------
         score : str
-             The score string (e.g., 'absorption', 'nu-fission')
+            The score string (e.g., 'absorption', 'nu-fission')
 
         Returns
         -------
-             The index in the Tally data array for this score.
+        score_index : int
+            The index in the Tally data array for this score.
 
         Raises
         ------
-             ValueError: An error when the argument passed to the 'score'
-             parameter cannot be found in the Tally.
+        ValueError
+            When the argument passed to the 'score' parameter cannot be found in
+            the Tally.
+
         """
 
         try:
@@ -708,7 +809,6 @@ class Tally(object):
 
         return score_index
 
-
     def get_values(self, scores=[], filters=[], filter_bins=[],
                    nuclides=[], value='mean'):
         """Returns a tally score value given a list of filters to satisfy.
@@ -720,42 +820,45 @@ class Tally(object):
         Parameters
         ----------
         scores : list
-             A list of one or more score strings
-             (e.g., ['absorption', 'nu-fission']; default is [])
+            A list of one or more score strings
+            (e.g., ['absorption', 'nu-fission']; default is [])
 
         filters : list
-             A list of filter type strings
-             (e.g., ['mesh', 'energy']; default is [])
+            A list of filter type strings
+            (e.g., ['mesh', 'energy']; default is [])
 
         filter_bins : list
-             A list of the filter bins corresponding to the filter_types
-             parameter (e.g., [1, (0., 0.625e-6)]; default is []). Each bin
-             in the list is the integer ID for 'material', 'surface', 'cell',
-             'cellborn', and 'universe' Filters. Each bin is an integer for
-             the cell instance ID for 'distribcell Filters. Each bin is
-             a 2-tuple of floats for 'energy' and 'energyout' filters
-             corresponding to the energy boundaries of the bin of interest.
-             The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding
-             to the mesh cell of interest. The order of the bins in the list
-             must correspond of the filter_types parameter.
+            A list of the filter bins corresponding to the filter_types
+            parameter (e.g., [1, (0., 0.625e-6)]; default is []). Each bin in
+            the list is the integer ID for 'material', 'surface', 'cell',
+            'cellborn', and 'universe' Filters. Each bin is an integer for the
+            cell instance ID for 'distribcell Filters. Each bin is a 2-tuple of
+            floats for 'energy' and 'energyout' filters corresponding to the
+            energy boundaries of the bin of interest.  The bin is a (x,y,z)
+            3-tuple for 'mesh' filters corresponding to the mesh cell of
+            interest. The order of the bins in the list must correspond of the
+            filter_types parameter.
 
         nuclides : list
-             A list of nuclide name strings
-             (e.g., ['U-235', 'U-238']; default is [])
+            A list of nuclide name strings
+            (e.g., ['U-235', 'U-238']; default is [])
 
         value : str
-             A string for the type of value to return  - 'mean' (default),
-             'std_dev', 'rel_err', 'sum', or 'sum_sq' are accepted
+            A string for the type of value to return  - 'mean' (default),
+            'std_dev', 'rel_err', 'sum', or 'sum_sq' are accepted
 
         Returns
         -------
-             A scalar or NumPy array of the Tally data indexed in the order
-             each filter, nuclide and score is listed in the parameters.
+        float or ndarray
+            A scalar or NumPy array of the Tally data indexed in the order
+            each filter, nuclide and score is listed in the parameters.
 
         Raises
         ------
-             ValueError : An error when this routine is called before the Tally
-             is populated with data by the StatePoint.read_results() routine.
+        ValueError
+            When this routine is called before the Tally is populated with data
+            by the StatePoint.read_results() routine.
+
         """
 
         # Ensure that StatePoint.read_results() was called first
@@ -772,13 +875,11 @@ class Tally(object):
         ############################      FILTERS      #########################
         # Determine the score indices from any of the requested scores
         if filters:
-
             # Initialize empty list of indices for each bin in each Filter
             filter_indices = []
 
             # Loop over all of the Tally's Filters
             for i, filter in enumerate(self.filters):
-
                 # Initialize empty list of indices for this Filter's bins
                 filter_indices.append([])
 
@@ -793,7 +894,6 @@ class Tally(object):
 
                 # If not a user-requested Filter, get all bins
                 if not user_filter:
-
                     # Create list of 2- or 3-tuples tuples for mesh cell bins
                     if filter.type == 'mesh':
                         dimension = filter.mesh.dimension
@@ -866,7 +966,6 @@ class Tally(object):
 
         return data.squeeze()
 
-
     def get_pandas_dataframe(self, filters=True, nuclides=True,
                              scores=True, summary=None):
         """Build a Pandas DataFrame for the Tally data.
@@ -880,31 +979,34 @@ class Tally(object):
         Parameters
         ----------
         filters : bool
-             Include columns with filter bin information (default is True).
+            Include columns with filter bin information (default is True).
 
         nuclides : bool
-             Include columns with nuclide bin information (default is True).
+            Include columns with nuclide bin information (default is True).
 
         scores : bool
-             Include columns with score bin information (default is True).
+            Include columns with score bin information (default is True).
 
         summary : None or Summary
-             An optional Summary object to be used to construct columns for
-             for distribcell tally filters (default is None). The geometric
-             information in the Summary object is embedded into a Multi-index
-             column with a geometric "path" to each distribcell intance.
-             NOTE: This option requires the OpenCG Python package.
+            An optional Summary object to be used to construct columns for for
+            distribcell tally filters (default is None). The geometric
+            information in the Summary object is embedded into a Multi-index
+            column with a geometric "path" to each distribcell intance.  NOTE:
+            This option requires the OpenCG Python package.
 
         Returns
         -------
-             A Pandas DataFrame with each column annotated by filter, nuclide
-             and score bin information (if these parameters are True), and the
-             mean and standard deviation of the Tally's data.
+        pandas.DataFrame
+            A Pandas DataFrame with each column annotated by filter, nuclide and
+            score bin information (if these parameters are True), and the mean
+            and standard deviation of the Tally's data.
 
         Raises
         ------
-             KeyError : An error when this routine is called before the Tally
-             is populated with data by the StatePoint.read_results() routine.
+        KeyError
+            When this routine is called before the Tally is populated with data
+            by the StatePoint.read_results() routine.
+
         """
 
         # Ensure that StatePoint.read_results() was called first
@@ -940,9 +1042,7 @@ class Tally(object):
 
         # Build DataFrame columns for filters if user requested them
         if filters:
-
             for filter in self.filters:
-
                 # mesh filters
                 if filter.type == 'mesh':
 
@@ -989,9 +1089,7 @@ class Tally(object):
 
                 # distribcell filters
                 elif filter.type == 'distribcell':
-
                     if isinstance(summary, Summary):
-
                         # Attempt to import the OpenCG package
                         try:
                             import opencg
@@ -1180,7 +1278,6 @@ class Tally(object):
         df = df.dropna(axis=1)
         return df
 
-
     def export_results(self, filename='tally-results', directory='.',
                       format='hdf5', append=True):
         """Exports tallly results to an HDF5 or Python pickle binary file.
@@ -1188,22 +1285,24 @@ class Tally(object):
         Parameters
         ----------
         filename : str
-             The name of the file for the results (default is 'tally-results')
+            The name of the file for the results (default is 'tally-results')
 
         directory : str
-             The name of the directory for the results (default is '.')
+            The name of the directory for the results (default is '.')
 
         format : str
-             The format for the exported file - HDF5 ('hdf5', default) and
-             Python pickle ('pkl') files are supported
+            The format for the exported file - HDF5 ('hdf5', default) and
+            Python pickle ('pkl') files are supported
 
         append : bool
-             Whether or not to append the results to the file (default is True)
+            Whether or not to append the results to the file (default is True)
 
         Raises
         ------
-             KeyError : An error when this routine is called before the Tally
-             is populated with data by the StatePoint.read_results() routine.
+        KeyError
+            When this routine is called before the Tally is populated with data
+            by the StatePoint.read_results() routine.
+
         """
 
         # Ensure that StatePoint.read_results() was called first
@@ -1239,10 +1338,8 @@ class Tally(object):
         if not os.path.exists(directory):
             os.makedirs(directory)
 
-
         # HDF5 binary file
         if format == 'hdf5':
-
             import h5py
 
             filename = directory + '/' + filename + '.h5'
@@ -1267,7 +1364,6 @@ class Tally(object):
             for nuclide in self.nuclides:
                 nuclides.append(nuclide.name)
 
-
             tally_group.create_dataset('nuclides', data=np.array(nuclides))
 
             # Create an HDF5 sub-group for the Filters
@@ -1285,12 +1381,8 @@ class Tally(object):
             # Close the Tally results HDF5 file
             tally_results.close()
 
-
         # Python pickle binary file
         elif format == 'pkl':
-
-            import pickle
-
             # Load the dictionary from the Pickle file
             filename = directory + '/' + filename + '.pkl'
 
@@ -1335,16 +1427,29 @@ class Tally(object):
 
 
 class TalliesFile(object):
+    """Tallies file used for an OpenMC simulation. Corresponds directly to the
+    tallies.xml input file.
+
+    """
 
     def __init__(self):
-
         # Initialize TalliesFile class attributes
         self._tallies = []
         self._meshes = []
         self._tallies_file = ET.Element("tallies")
 
-
     def add_tally(self, tally, merge=False):
+        """Add a tally to the file
+
+        Parameters
+        ----------
+        tally : Tally
+            Tally to add to file
+        merge : bool
+            Indicate whether the tally should be merged with an existing tally,
+            if possible. Defaults to False.
+
+        """
 
         if not isinstance(tally, Tally):
             msg = 'Unable to add a non-Tally {0} to the TalliesFile'.format(tally)
@@ -1358,7 +1463,6 @@ class TalliesFile(object):
 
                 # If a mergeable tally is found
                 if tally2.can_merge(tally):
-
                     # Replace tally 2 with the merged tally
                     merged_tally = tally2.merge(tally)
                     self._tallies[i] = merged_tally
@@ -1372,23 +1476,32 @@ class TalliesFile(object):
         else:
             self._tallies.append(tally)
 
-
     def remove_tally(self, tally):
+        """Remove a tally from the file
+
+        Parameters
+        ----------
+        tally : Tally
+            Tally to remove
+
+        """
+
         self._tallies.remove(tally)
 
-
     def merge_tallies(self):
+        """Merge any mergable tallies together. Note that n-way merges are
+        possible.
+
+        """
 
         for i, tally1 in enumerate(self._tallies):
             for j, tally2 in enumerate(self._tallies):
-
                 # Do not merge the same tally with itself
                 if i == j:
                     continue
 
                 # If the two tallies are mergeable
                 if tally1.can_merge(tally2):
-
                     # Replace tally 1 with the merged tally
                     merged_tally = tally1.merge(tally2)
                     self._tallies[i] = merged_tally
@@ -1399,8 +1512,15 @@ class TalliesFile(object):
                     # Continue iterating from the first loop
                     break
 
-
     def add_mesh(self, mesh):
+        """Add a mesh to the file
+
+        Parameters
+        ----------
+        mesh : openmc.mesh.Mesh
+            Mesh to add to the file
+
+        """
 
         if not isinstance(mesh, Mesh):
             msg = 'Unable to add a non-Mesh {0} to the TalliesFile'.format(mesh)
@@ -1408,33 +1528,38 @@ class TalliesFile(object):
 
         self._meshes.append(mesh)
 
-
     def remove_mesh(self, mesh):
+        """Remove a mesh from the file
+
+        Parameters
+        ----------
+        mesh : openmc.mesh.Mesh
+            Mesh to remove from the file
+
+        """
+
         self._meshes.remove(mesh)
 
-
-    def create_tally_subelements(self):
-
+    def _create_tally_subelements(self):
         for tally in self._tallies:
             xml_element = tally.get_tally_xml()
             self._tallies_file.append(xml_element)
 
-
-    def create_mesh_subelements(self):
-
+    def _create_mesh_subelements(self):
         for mesh in self._meshes:
-
             if len(mesh._name) > 0:
                 self._tallies_file.append(ET.Comment(mesh._name))
 
             xml_element = mesh.get_mesh_xml()
             self._tallies_file.append(xml_element)
 
-
     def export_to_xml(self):
+        """Create a tallies.xml file that can be used for a simulation.
 
-        self.create_mesh_subelements()
-        self.create_tally_subelements()
+        """
+
+        self._create_mesh_subelements()
+        self._create_tally_subelements()
 
         # Clean the indentation in the file to be user-readable
         clean_xml_indentation(self._tallies_file)