Merge remote-tracking branch 'upstream/develop' into valgrind

2026-07-27 13:45:36 -04:00 · 2019-12-04 18:48:18 -05:00 · 2019-12-04 18:48:18 -05:00 · eabf3e2a86
commit eabf3e2a86
parent 57cab6b346 8fcd74b14b
241 changed files with 28672 additions and 11069 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -114,7 +114,11 @@ endif()
 #===============================================================================

 add_library(pugixml vendor/pugixml/pugixml.cpp)
-target_include_directories(pugixml PUBLIC vendor/pugixml/)
+target_include_directories(pugixml
+  PUBLIC
+    $<INSTALL_INTERFACE:include/pugixml>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/vendor/pugixml>
+)

 #===============================================================================
 # xtensor header-only library
@ -126,6 +130,7 @@ if (NOT (CMAKE_VERSION VERSION_LESS 3.13))
 endif()

 add_subdirectory(vendor/xtl)
+set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl)
 add_subdirectory(vendor/xtensor)
 target_link_libraries(xtensor INTERFACE xtl)

@ -134,7 +139,11 @@ target_link_libraries(xtensor INTERFACE xtl)
 #===============================================================================

 add_library(gsl INTERFACE)
-target_include_directories(gsl INTERFACE vendor/gsl/include)
+target_include_directories(gsl
+  INTERFACE
+    $<INSTALL_INTERFACE:include/gsl/include>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/vendor/gsl/include>
+)

 # Make sure contract violations throw exceptions
 target_compile_definitions(gsl INTERFACE GSL_THROW_ON_CONTRACT_VIOLATION)
@ -143,6 +152,9 @@ target_compile_definitions(gsl INTERFACE GSL_THROW_ON_CONTRACT_VIOLATION)
 # RPATH information
 #===============================================================================

+# Provide install directory variables as defined by GNU coding standards
+include(GNUInstallDirs)
+
 # This block of code ensures that dynamic libraries can be found via the RPATH
 # whether the executable is the original one from the build directory or the
 # installed one in CMAKE_INSTALL_PREFIX. Ref:
@ -155,16 +167,14 @@ set(CMAKE_SKIP_BUILD_RPATH  FALSE)
 # (but later on when installing)
 set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)

-set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
-
 # add the automatically determined parts of the RPATH
 # which point to directories outside the build tree to the install RPATH
 set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)

 # the RPATH to be used when installing, but only if it's not a system directory
-list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir)
+list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_FULL_LIBDIR}" isSystemDir)
 if("${isSystemDir}" STREQUAL "-1")
-  set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+  set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_FULL_LIBDIR}")
 endif()

 #===============================================================================
@ -172,7 +182,11 @@ endif()
 #===============================================================================

 add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.cc)
-target_include_directories(faddeeva PUBLIC vendor/faddeeva/)
+target_include_directories(faddeeva
+  PUBLIC
+    $<INSTALL_INTERFACE:include/faddeeva>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/vendor/faddeeva>
+)
 target_compile_options(faddeeva PRIVATE ${cxxflags})

 #===============================================================================
@ -283,7 +297,11 @@ set_target_properties(libopenmc PROPERTIES
  OUTPUT_NAME openmc)

 target_include_directories(libopenmc
-  PUBLIC include ${HDF5_INCLUDE_DIRS})
+  PUBLIC
+    $<INSTALL_INTERFACE:include>
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+    ${HDF5_INCLUDE_DIRS}
+)

 # Set compile flags
 target_compile_options(libopenmc PRIVATE ${cxxflags})
@ -343,12 +361,27 @@ add_custom_command(TARGET libopenmc POST_BUILD
 # Install executable, scripts, manpage, license
 #===============================================================================

-install(TARGETS openmc libopenmc
-  RUNTIME DESTINATION bin
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  )
-install(DIRECTORY src/relaxng DESTINATION share/openmc)
-install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
-install(FILES LICENSE DESTINATION "share/doc/openmc" RENAME copyright)
-install(DIRECTORY include/ DESTINATION include)
+set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
+install(TARGETS openmc libopenmc pugixml faddeeva gsl
+  EXPORT openmc-targets
+  RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+  LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+  ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+)
+install(EXPORT openmc-targets
+  FILE OpenMCTargets.cmake
+  NAMESPACE OpenMC::
+  DESTINATION ${INSTALL_CONFIGDIR})
+
+install(DIRECTORY src/relaxng DESTINATION ${CMAKE_INSTALL_DATADIR}/openmc)
+install(FILES cmake/OpenMCConfig.cmake DESTINATION ${INSTALL_CONFIGDIR})
+install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
+install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright)
+install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+
+# Copy headers for vendored dependencies (note that xtensor/xtl are handled
+# separately since they are managed by CMake)
+install(DIRECTORY vendor/pugixml DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+  FILES_MATCHING PATTERN "*.hpp")
+install(DIRECTORY vendor/gsl DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(DIRECTORY vendor/faddeeva DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
--- a/48
+++ b/48
@ -0,0 +1,48 @@
+# Data interface
+openmc/data/ @paulromano
+
+# Python bindings to C/C++ API
+openmc/lib/ @paulromano
+
+# Depletion
+openmc/deplete/ @drewejohnson
+tests/regression_tests/deplete/ @drewejohnson
+tests/unit_tests/test_deplete_*.py @drewejohnson
+
+# MG-related functionality
+openmc/mgxs_library.py @nelsonag
+src/mgxs.cpp @nelsonag
+src/mgxs_interface.cpp @nelsonag
+src/physics_mg.cpp @nelsonag
+src/scattdata.cpp @nelsonag
+src/xsdata.cpp @nelsonag
+
+# CMFD
+openmc/cmfd.py @shikhar413
+src/cmfd_solver.cpp @shikhar413
+
+# DAGMC
+src/dagmc.cpp @pshriwise
+tests/regression_tests/dagmc/ @pshriwise
+tests/unit_tests/dagmc/ @pshriwise
+
+# Photon transport
+openmc/data/BREMX.DAT @amandalund
+openmc/data/compton_profiles.h5 @amandalund
+openmc/data/photon.py @amandalund
+src/photon.cpp @amandalund
+src/bremsstrahlung.cpp @amandalund
+tests/regression_tests/photon_production/ @amandalund
+tests/regression_tests/photon_source/ @amandalund
+
+# RCP and TRISOs
+openmc/model/triso.py @amandalund
+tests/regression_tests/triso/ @amandalund
+tests/unit_tests/test_model_triso.py @amandalund
+
+# Geometry plotting
+src/plot.cpp @pshriwise
+openmc/lib/plot.py @pshriwise
+
+# Resonance covariance
+openmc/data/resonance_covariance.py @icmeyer
--- a/2
+++ b/2
@ -1,4 +1,4 @@
-Copyright (c) 2011-2018 Massachusetts Institute of Technology and OpenMC contributors
+Copyright (c) 2011-2019 Massachusetts Institute of Technology and OpenMC contributors

 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in
--- a/cmake/OpenMCConfig.cmake
+++ b/cmake/OpenMCConfig.cmake
@ -0,0 +1,8 @@
+get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
+
+find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
+find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
+
+if(NOT TARGET OpenMC::libopenmc)
+  include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
+endif()
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@ -27,7 +27,7 @@ MOCK_MODULES = [
    'scipy.interpolate', 'scipy.integrate', 'scipy.optimize', 'scipy.special',
    'scipy.stats', 'scipy.spatial', 'h5py', 'pandas', 'uncertainties',
    'matplotlib', 'matplotlib.pyplot', 'openmoc',
-    'openmc.data.reconstruct'
+    'openmc.data.reconstruct', 'openmc.checkvalue'
 ]
 sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)

@ -79,9 +79,9 @@ copyright = '2011-2019, Massachusetts Institute of Technology and OpenMC contrib
 # built documents.
 #
 # The short X.Y version.
-version = "0.11"
+version = "0.12"
 # The full version, including alpha/beta/rc tags.
-release = "0.11.0-dev"
+release = "0.12.0-dev"

 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
--- a/docs/source/devguide/docbuild.rst
+++ b/docs/source/devguide/docbuild.rst
@ -5,27 +5,28 @@ Building Sphinx Documentation
 =============================

 In order to build the documentation in the ``docs`` directory, you will need to
-have the `Sphinx <http://openmc.readthedocs.io/en/latest/>`_ third-party Python
+have the `Sphinx <https://www.sphinx-doc.org/en/master/>`_ third-party Python
 package. The easiest way to install Sphinx is via pip:

 .. code-block:: sh

-    sudo pip install sphinx
+    pip install sphinx

-Additionally, you will also need a Sphinx extension for numbering figures. The
-`Numfig <http://openmc.readthedocs.io/en/latest/>`_ package can be installed
+Additionally, you will need several Sphinx extensions that can be installed
 directly with pip:

 .. code-block:: sh

-   sudo pip install sphinx-numfig
+    pip install sphinx-numfig
+    pip install sphinxcontrib-katex
+    pip install sphinxcontrib-svg2pdfconverter

 -----------------------------------
 Building Documentation as a Webpage
 -----------------------------------

 To build the documentation as a webpage (what appears at
-http://openmc.readthedocs.io), simply go to the ``docs`` directory and run:
+https://docs.openmc.org), simply go to the ``docs`` directory and run:

 .. code-block:: sh

--- a/docs/source/devguide/docker.rst
+++ b/docs/source/devguide/docker.rst
@ -19,7 +19,7 @@ build a Docker image with OpenMC installed. The image includes OpenMC with
 MPICH and parallel HDF5 in the ``/opt/openmc`` directory, and
 `Miniconda3 <https://conda.io/miniconda.html>`_ with all of the Python
 pre-requisites (NumPy, SciPy, Pandas, etc.) installed. The
-`NJOY2016 <http://www.njoy21.io/NJOY2016/>`_ codebase is installed in
+`NJOY2016 <https://www.njoy21.io/NJOY2016/>`_ codebase is installed in
 ``/opt/NJOY2016`` to support full functionality and testing of the
 ``openmc.data`` Python module. The publicly available nuclear data libraries
 necessary to run OpenMC's test suite -- including NNDC and WMP cross sections
@ -54,4 +54,3 @@ Docker container where you have access to use OpenMC.
 .. _Docker container: https://www.docker.com/resources/what-container
 .. _options: https://docs.docker.com/engine/reference/commandline/run/
 .. _mounting volumes: https://docs.docker.com/storage/volumes/
-
--- a/docs/source/devguide/styleguide.rst
+++ b/docs/source/devguide/styleguide.rst
@ -207,8 +207,8 @@ Documentation
 -------------

 Classes, structs, and functions are to be annotated for the `Doxygen
-<http://www.stack.nl/~dimitri/doxygen/>`_ documentation generation tool. Use the
-``\`` form of Doxygen commands, e.g., ``\brief`` instead of ``@brief``.
+<http://www.doxygen.nl/>`_ documentation generation tool. Use the ``\`` form of
+Doxygen commands, e.g., ``\brief`` instead of ``@brief``.

 ------
 Python
@ -231,7 +231,7 @@ represent a filesystem path should work with both strings and Path_ objects.
 .. _C++ Core Guidelines: http://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines
 .. _PEP8: https://www.python.org/dev/peps/pep-0008/
 .. _numpydoc: https://numpydoc.readthedocs.io/en/latest/format.html
-.. _numpy: http://www.numpy.org/
+.. _numpy: https://numpy.org/
 .. _scipy: https://www.scipy.org/
 .. _matplotlib: https://matplotlib.org/
 .. _pandas: https://pandas.pydata.org/
--- a/docs/source/devguide/workflow.rst
+++ b/docs/source/devguide/workflow.rst
@ -124,7 +124,7 @@ can interfere with virtual environments.
 .. _GitHub: https://github.com/
 .. _git flow: http://nvie.com/git-model
 .. _valgrind: http://valgrind.org/
-.. _style guide: http://openmc.readthedocs.io/en/latest/devguide/styleguide.html
+.. _style guide: https://docs.openmc.org/en/latest/devguide/styleguide.html
 .. _pull request: https://help.github.com/articles/using-pull-requests
 .. _openmc-dev/openmc: https://github.com/openmc-dev/openmc
 .. _paid plan: https://github.com/plans
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@ -2,12 +2,14 @@
 The OpenMC Monte Carlo Code
 ===========================

-OpenMC is a Monte Carlo particle transport simulation code focused on neutron
-criticality calculations. It is capable of simulating 3D models based on
-constructive solid geometry with second-order surfaces. OpenMC supports either
-continuous-energy or multi-group transport.  The continuous-energy particle
+OpenMC is a community-developed Monte Carlo neutron and photon transport
+simulation code. It is capable of performing fixed source, k-eigenvalue, and
+subcritical multiplication calculations on models built using either a
+constructive solid geometry or CAD representation. OpenMC supports both
+continuous-energy and multigroup transport. The continuous-energy particle
 interaction data is based on a native HDF5 format that can be generated from ACE
-files used by the MCNP and Serpent Monte Carlo codes.
+files produced by NJOY. Parallelism is enabled via a hybrid MPI and OpenMP
+programming model.

 OpenMC was originally developed by members of the `Computational Reactor Physics
 Group <http://crpg.mit.edu>`_ at the `Massachusetts Institute of Technology
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@ -333,7 +333,7 @@ or sub-elements:
    velocity sampling) or "dbrc" (Doppler broadening rejection correction).
    Descriptions of each of these methods are documented here_.

-    .. _here: http://dx.doi.org/10.1016/j.anucene.2017.12.044
+    .. _here: https://doi.org/10.1016/j.anucene.2017.12.044

    *Default*: "rvs"

@ -423,12 +423,16 @@ attributes/sub-elements:

    :type:
      The type of spatial distribution. Valid options are "box", "fission",
-      "point", and "cartesian". A "box" spatial distribution has coordinates
-      sampled uniformly in a parallelepiped. A "fission" spatial distribution
-      samples locations from a "box" distribution but only locations in
-      fissionable materials are accepted. A "point" spatial distribution has
-      coordinates specified by a triplet. An "cartesian" spatial distribution
-      specifies independent distributions of x-, y-, and z-coordinates.
+      "point", "cartesian", and "spherical". A "box" spatial distribution has 
+      coordinates sampled uniformly in a parallelepiped. A "fission" spatial
+      distribution samples locations from a "box" distribution but only 
+      locations in fissionable materials are accepted. A "point" spatial 
+      distribution has coordinates specified by a triplet. An "cartesian" 
+      spatial distribution specifies independent distributions of x-, y-, and
+      z-coordinates. A "spherical" spatial distribution specifies independent
+      distributions of r-, theta-, and phi-coordinates where theta is the angle
+      with respect to the z-axis, phi is the azimuthal angle, and the sphere is
+      centered on the coordinate (x0,y0,z0).

      *Default*: None

@ -446,6 +450,9 @@ attributes/sub-elements:
      For an "cartesian" distribution, no parameters are specified. Instead,
      the ``x``, ``y``, and ``z`` elements must be specified.

+      For a "spherical" distribution, no parameters are specified. Instead,
+      the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified.
+
      *Default*: None

    :x:
@ -466,6 +473,28 @@ attributes/sub-elements:
      univariate probability distribution (see the description in
      :ref:`univariate`).

+    :r:
+      For a "spherical" distribution, this element specifies the distribution
+      of r-coordinates. The necessary sub-elements/attributes are those of a
+      univariate probability distribution (see the description in
+      :ref:`univariate`).
+
+    :theta:
+      For a "spherical" distribution, this element specifies the distribution
+      of theta-coordinates. The necessary sub-elements/attributes are those of a
+      univariate probability distribution (see the description in
+      :ref:`univariate`).
+
+    :phi:
+      For a "spherical" distribution, this element specifies the distribution
+      of phi-coordinates. The necessary sub-elements/attributes are those of a
+      univariate probability distribution (see the description in
+      :ref:`univariate`).
+
+    :origin:
+      For a "spherical" distribution, this element specifies the coordinates of
+      the center of the sphere.
+
  :angle:
    An element specifying the angular distribution of source sites. This element
    has the following attributes:
--- a/docs/source/io_formats/statepoint.rst
+++ b/docs/source/io_formats/statepoint.rst
@ -109,6 +109,11 @@ The current version of the statepoint file format is 17.0.

 **/tallies/tally <uid>/**

+:Attributes:
+             - **internal** (*int*) -- Flag indicating the presence of tally
+               data (0) or absence of tally data (1). All user defined
+               tallies will have a value of 0 unless otherwise instructed.
+
 :Datasets: - **n_realizations** (*int*) -- Number of realizations.
           - **n_filters** (*int*) -- Number of filters used.
           - **filters** (*int[]*) -- User-defined unique IDs of the filters on
--- a/docs/source/io_formats/volume.rst
+++ b/docs/source/io_formats/volume.rst
@ -23,6 +23,8 @@ The current version of the volume file format is 1.0.
               bounding box
             - **upper_right** (*double[3]*) -- Upper-right coordinates of
               bounding box
+             - **threshold** (*double*) -- Threshold used for volume uncertainty
+             - **trigger_type** (*char[]*) -- Trigger type used for volume uncertainty

 **/domain_<id>/**

--- a/docs/source/methods/cmfd.rst
+++ b/docs/source/methods/cmfd.rst
@ -538,7 +538,7 @@ Examples of CMFD simulations using OpenMC can be found in [HermanThesis]_.
   .. rubric:: References

 .. [BEAVRS] Nick Horelik, Bryan Herman. *Benchmark for Evaluation And Verification of Reactor
-            Simulations*. Massachusetts Institute of Technology, http://crpg.mit.edu/pub/beavrs
+            Simulations*. Massachusetts Institute of Technology, https://crpg.mit.edu/research/beavrs
            , 2013.

 .. [Gill] Daniel F. Gill. *Newton-Krylov methods for the solution of the k-eigenvalue problem in
--- a/docs/source/methods/cross_sections.rst
+++ b/docs/source/methods/cross_sections.rst
@ -8,13 +8,14 @@ Cross Section Representations
 Continuous-Energy Data
 ----------------------

-The data governing the interaction of neutrons with
-various nuclei for continous-energy problems are represented using the ACE
-format which is used by MCNP_ and Serpent_. ACE-format data can be generated
-with the NJOY_ nuclear data processing system which converts raw
-`ENDF/B data`_ into linearly-interpolable data as required by most Monte Carlo
-codes. The use of a standard cross section format allows for a direct comparison
-of OpenMC with other codes since the same cross section libraries can be used.
+In OpenMC, the data governing the interaction of neutrons with various nuclei
+for continous-energy problems are represented using an HDF5 format that can be
+produced by converting files in the ACE format, which is used by MCNP_ and
+Serpent_. ACE-format data can be generated with the NJOY_ nuclear data
+processing system, which converts raw `ENDF/B data`_ into linearly-interpolable
+data as required by most Monte Carlo codes. Since ACE-format data can be
+converted into OpenMC's HDF5 format, it is possible to perform direct comparison
+of OpenMC with other codes using the same underlying nuclear data library.

 The ACE format contains continuous-energy cross sections for the following types
 of reactions: elastic scattering, fission (or first-chance fission,
@ -31,7 +32,7 @@ data can be used.
 Energy Grid Methods
 -------------------

-The method by which continuous energy cross sections for each nuclide in a
+The method by which continuous-energy cross sections for each nuclide in a
 problem are stored as a function of energy can have a substantial effect on the
 performance of a Monte Carlo simulation. Since the ACE format is based on
 linearly-interpolable cross sections, each nuclide has cross sections tabulated
@ -72,9 +73,9 @@ Windowed Multipole Representation
 ---------------------------------

 In addition to the usual pointwise representation of cross sections, OpenMC
-offers support for an experimental data format called windowed multipole (WMP).
-This data format requires less memory than pointwise cross sections, and it
-allows on-the-fly Doppler broadening to arbitrary temperature.
+offers support for a data format called windowed multipole (WMP). This data
+format requires less memory than pointwise cross sections, and it allows
+on-the-fly Doppler broadening to arbitrary temperature.

 The multipole method was introduced by Hwang_ and the faster windowed multipole
 method by Josey_. In the multipole format, cross section resonances are
@ -258,7 +259,7 @@ where a material has a very large cross sections relative to the other material
 used to minimize this error.

 Finally, the above options for representing the physics do not have to be
-consistent across the problem.  The number of groups and the structure, however,
+consistent across the problem. The number of groups and the structure, however,
 does have to be consistent across the data sets. That is to say that each
 microscopic or macroscopic data set does not have to apply the same scattering
 expansion, treatment of multiplicity or angular representation of the cross
@ -269,11 +270,11 @@ or even isotropic scattering.
 .. _logarithmic mapping technique:
   https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
 .. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381
-.. _Josey: http://dx.doi.org/10.1016/j.jcp.2015.08.013
+.. _Josey: https://doi.org/10.1016/j.jcp.2015.08.013
 .. _WMP Library: https://github.com/mit-crpg/WMP_Library
 .. _MCNP: http://mcnp.lanl.gov
 .. _Serpent: http://montecarlo.vtt.fi
 .. _NJOY: http://t2.lanl.gov/codes.shtml
 .. _ENDF/B data: http://www.nndc.bnl.gov/endf
-.. _Leppanen: http://dx.doi.org/10.1016/j.anucene.2009.03.019
+.. _Leppanen: https://doi.org/10.1016/j.anucene.2009.03.019
 .. _algorithms: http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package
--- a/docs/source/methods/energy_deposition.rst
+++ b/docs/source/methods/energy_deposition.rst
@ -4,11 +4,11 @@
 Heating and Energy Deposition
 =============================

-As particles traverse a problem, some portion of their energy is deposited at 
+As particles traverse a problem, some portion of their energy is deposited at
 collision sites. This energy is deposited when charged particles, including
 electrons and recoil nuclei, undergo electromagnetic interactions with
 surrounding electons and ions. The information describing how much energy
-is deposited for a specific reaction is referred to as 
+is deposited for a specific reaction is referred to as
 "heating numbers" and can be computed using a program like NJOY with the
 ``heatr`` module.

@ -108,7 +108,7 @@ Neutron Transport
 For this case, OpenMC instructs ``heatr`` to produce heating coefficients
 assuming that energy from photons, :math:`E_{\gamma, p}` and
 :math:`E_{\gamma, d}`, is deposited at the fission site.
-Let :math:`N901` represent the total heating number returned from this ``heatr`` 
+Let :math:`N901` represent the total heating number returned from this ``heatr``
 run with :math:`N918` reflecting fission heating computed from NJOY.
 :math:`M901` represent the following modification

@ -119,7 +119,7 @@ run with :math:`N918` reflecting fission heating computed from NJOY.
      + E_{i, \gamma, d}\right]\sigma_{i, f}(E).

 This modified heating data is stored as the MT=901 reaction and will be scored
-if ``901`` is included in :attr:`openmc.Tally.scores`.
+if ``heating-local`` is included in :attr:`openmc.Tally.scores`.

 Coupled neutron-photon transport
 --------------------------------
@ -146,4 +146,4 @@ References
 .. [Mack97] Abdou, M.A., Maynard, C.W., and Wright, R.Q. MACK: computer
   program to calculate neutron energy release parameters (fluence-to-kerma
   factors) and multigroup neutron reaction cross sections from nuclear data
-   in ENDF Format. Oak Ridge National Laboratory report ORNL-TM-3994.
+   in ENDF Format. Oak Ridge National Laboratory report ORNL-TM-3994.
--- a/docs/source/methods/geometry.rst
+++ b/docs/source/methods/geometry.rst
@ -960,8 +960,8 @@ normal using the equations from :ref:`transform-coordinates`. The white boundary
 condition can be applied to any kind of surface, as long as the normal to the
 surface is known as in :ref:`reflection`.

-.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
-.. _surfaces: http://en.wikipedia.org/wiki/Surface
+.. _constructive solid geometry: https://en.wikipedia.org/wiki/Constructive_solid_geometry
+.. _surfaces: https://en.wikipedia.org/wiki/Surface
 .. _MCNP: http://mcnp.lanl.gov
 .. _Serpent: http://montecarlo.vtt.fi
 .. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc
--- a/docs/source/methods/introduction.rst
+++ b/docs/source/methods/introduction.rst
@ -139,9 +139,9 @@ be performed before the run is finished. This include the following:

  - If requested, a source file is written to disk.

-  - All allocatable arrays are deallocated.
+  - Dynamically-allocated memory should be freed.

-.. _probability distributions: http://en.wikipedia.org/wiki/Probability_distribution
-.. _Monte Carlo: http://en.wikipedia.org/wiki/Monte_Carlo_method
-.. _central limit theorem: http://en.wikipedia.org/wiki/Central_limit_theorem
-.. _pseudorandom number: http://en.wikipedia.org/wiki/Pseudorandom_number_generator
+.. _probability distributions: https://en.wikipedia.org/wiki/Probability_distribution
+.. _Monte Carlo: https://en.wikipedia.org/wiki/Monte_Carlo_method
+.. _central limit theorem: https://en.wikipedia.org/wiki/Central_limit_theorem
+.. _pseudorandom number: https://en.wikipedia.org/wiki/Pseudorandom_number_generator
--- a/docs/source/methods/neutron_physics.rst
+++ b/docs/source/methods/neutron_physics.rst
@ -1298,11 +1298,10 @@ section over the range of velocities considered:
 where it should be noted that the maximum is taken over the range :math:`[v_n -
 4/\beta, 4_n + 4\beta]`. This method is known as Doppler broadening rejection
 correction (DBRC) and was first introduced by `Becker et al.`_. OpenMC has an
-implementation of DBRC as well as an accelerated sampling method that are
-described fully in `Walsh et al.`_
+implementation of DBRC as well as an accelerated sampling method that samples the `relative velocity`_ directly.

-.. _Becker et al.: http://dx.doi.org/10.1016/j.anucene.2008.12.001
-.. _Walsh et al.: http://dx.doi.org/10.1016/j.anucene.2014.01.017
+.. _Becker et al.: https://doi.org/10.1016/j.anucene.2008.12.001
+.. _relative velocity: https://doi.org/10.1016/j.anucene.2017.12.044

 .. _sab_tables:

@ -1645,23 +1644,23 @@ another.

 .. |sab| replace:: S(:math:`\alpha,\beta,T`)

-.. _SIGMA1 method: http://dx.doi.org/10.13182/NSE76-1
+.. _SIGMA1 method: https://doi.org/10.13182/NSE76-1

 .. _scaled interpolation: http://www.ans.org/pubs/journals/nse/a_26575

-.. _probability table method: http://dx.doi.org/10.13182/NSE72-3
+.. _probability table method: https://doi.org/10.13182/NSE72-3

-.. _Watt fission spectrum: http://dx.doi.org/10.1103/PhysRev.87.1037
+.. _Watt fission spectrum: https://doi.org/10.1103/PhysRev.87.1037

 .. _Foderaro: http://hdl.handle.net/1721.1/1716

-.. _OECD: http://www.oecd-nea.org/dbprog/MMRW-BOOKS.html
+.. _OECD: http://www.oecd-nea.org/tools/abstract/detail/NEA-1792

-.. _NJOY: https://njoy.github.io/NJOY2016/
+.. _NJOY: https://www.njoy21.io/NJOY2016/

 .. _PREPRO: http://www-nds.iaea.org/ndspub/endf/prepro/

-.. _ENDF-6 Format: http://www-nds.iaea.org/ndspub/documents/endf/endf102/endf102.pdf
+.. _ENDF-6 Format: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf

 .. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf

@ -1669,7 +1668,7 @@ another.

 .. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf

-.. _Romano: http://dx.doi.org/10.1016/j.cpc.2014.11.001
+.. _Romano: https://doi.org/10.1016/j.cpc.2014.11.001

 .. _Sutton and Brown: http://www.osti.gov/bridge/product.biblio.jsp?osti_id=307911

--- a/docs/source/methods/parallelization.rst
+++ b/docs/source/methods/parallelization.rst
@ -607,33 +607,33 @@ is actually independent of the number of nodes:
   Radiation Penetration Calculations on a Parallel Computer,"
   *Trans. Am. Nucl. Soc.*, **17**, 260 (1973).

-.. _first paper: http://www.jstor.org/stable/2280232
+.. _first paper: https://doi.org/10.2307/2280232

 .. _work of Forrest Brown: http://hdl.handle.net/2027.42/24996

-.. _Brissenden and Garlick: http://dx.doi.org/10.1016/0306-4549(86)90095-2
+.. _Brissenden and Garlick: https://doi.org/10.1016/0306-4549(86)90095-2

 .. _MPICH2: http://www.mcs.anl.gov/mpi/mpich

-.. _binomial tree: http://www.cs.auckland.ac.nz/~jmor159/PLDS210/trees.html
+.. _binomial tree: https://www.mcs.anl.gov/~thakur/papers/ijhpca-coll.pdf

-.. _Geary: http://www.jstor.org/stable/10.2307/2342070
+.. _Geary: https://doi.org/10.2307/2342070

 .. _Barnett: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.51.7772

-.. _single-instruction multiple-data: http://en.wikipedia.org/wiki/SIMD
+.. _single-instruction multiple-data: https://en.wikipedia.org/wiki/SIMD

-.. _vector computers: http://en.wikipedia.org/wiki/Vector_processor
+.. _vector computers: https://en.wikipedia.org/wiki/Vector_processor

-.. _single program multiple data: http://en.wikipedia.org/wiki/SPMD
+.. _single program multiple data: https://en.wikipedia.org/wiki/SPMD

-.. _message-passing interface: http://en.wikipedia.org/wiki/Message_Passing_Interface
+.. _message-passing interface: https://en.wikipedia.org/wiki/Message_Passing_Interface

 .. _PVM: http://www.csm.ornl.gov/pvm/pvm_home.html

 .. _MPI: http://www.mcs.anl.gov/research/projects/mpi/

-.. _embarrassingly parallel: http://en.wikipedia.org/wiki/Embarrassingly_parallel
+.. _embarrassingly parallel: https://en.wikipedia.org/wiki/Embarrassingly_parallel

 .. _sends: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Send.html

@ -643,8 +643,8 @@ is actually independent of the number of nodes:

 .. _allgather: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Allgather.html

-.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution
+.. _Cauchy distribution: https://en.wikipedia.org/wiki/Cauchy_distribution

-.. _latency: http://en.wikipedia.org/wiki/Latency_(engineering)#Packet-switched_networks
+.. _latency: https://en.wikipedia.org/wiki/Latency_(engineering)#Packet-switched_networks

-.. _bandwidth: http://en.wikipedia.org/wiki/Bandwidth_(computing)
+.. _bandwidth: https://en.wikipedia.org/wiki/Bandwidth_(computing)
--- a/docs/source/methods/random_numbers.rst
+++ b/docs/source/methods/random_numbers.rst
@ -67,6 +67,6 @@ the idea is to determine the new multiplicative and additive constants in
   .. rubric:: References


-.. _L'Ecuyer: http://dx.doi.org/10.1090/S0025-5718-99-00996-5
+.. _L'Ecuyer: https://doi.org/10.1090/S0025-5718-99-00996-5
 .. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf
-.. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator
+.. _linear congruential generator: https://en.wikipedia.org/wiki/Linear_congruential_generator
--- a/docs/source/methods/tallies.rst
+++ b/docs/source/methods/tallies.rst
@ -484,31 +484,31 @@ improve the estimate of the percentile.

   .. rubric:: References

-.. _following approximation: http://dx.doi.org/10.1080/03610918708812641
+.. _following approximation: https://doi.org/10.1080/03610918708812641

-.. _Bessel's correction: http://en.wikipedia.org/wiki/Bessel's_correction
+.. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction

-.. _random variable: http://en.wikipedia.org/wiki/Random_variable
+.. _random variable: https://en.wikipedia.org/wiki/Random_variable

-.. _stochastic process: http://en.wikipedia.org/wiki/Stochastic_process
+.. _stochastic process: https://en.wikipedia.org/wiki/Stochastic_process

-.. _independent, identically-distributed random variables: http://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables
+.. _independent, identically-distributed random variables: https://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables

-.. _law of large numbers: http://en.wikipedia.org/wiki/Law_of_large_numbers
+.. _law of large numbers: https://en.wikipedia.org/wiki/Law_of_large_numbers

-.. _expected value: http://en.wikipedia.org/wiki/Expected_value
+.. _expected value: https://en.wikipedia.org/wiki/Expected_value

-.. _converges in probability: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability
+.. _converges in probability: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability

-.. _normal distribution: http://en.wikipedia.org/wiki/Normal_distribution
+.. _normal distribution: https://en.wikipedia.org/wiki/Normal_distribution

-.. _converges in distribution: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
+.. _converges in distribution: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution

-.. _confidence intervals: http://en.wikipedia.org/wiki/Confidence_interval
+.. _confidence intervals: https://en.wikipedia.org/wiki/Confidence_interval

-.. _Student's t-distribution: http://en.wikipedia.org/wiki/Student%27s_t-distribution
+.. _Student's t-distribution: https://en.wikipedia.org/wiki/Student%27s_t-distribution

-.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution
+.. _Cauchy distribution: https://en.wikipedia.org/wiki/Cauchy_distribution

 .. _unpublished rational approximation: https://web.archive.org/web/20150926021742/http://home.online.no/~pjacklam/notes/invnorm/

--- a/docs/source/publications.rst
+++ b/docs/source/publications.rst
@ -142,6 +142,11 @@ Geometry and Visualization
 Miscellaneous
 -------------

+- Govatsa Acharya, "`Investigating the Application of Self-Actuated Passive
+  Shutdown System in a Small Lead-Cooled Reactor
+  <https://doi.org/10.13140/RG.2.2.26088.01281>`_," M.S. Thesis, KTH Royal
+  Institute of Technology (2019).
+
 - Shikhar Kumar, Benoit Forget, and Kord Smith, "Analysis of fission source
  convergence for a 3-D SMR core using functional expansion tallies," *Proc.
  M&C*, 937-947, Portland, Oregon, Aug. 25-29 (2019).
@ -366,7 +371,7 @@ Nuclear Data
 - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown,
  "`Uncertainty in Fast Reactor-Relevant Critical Benchmark Simulations Due to
  Unresolved Resonance Structure
-  <https://www.kns.org/paper_file/paper/MC2017_2017_3/P197S03-09WalshJ.pdf>`_,"
+  <https://www.kns.org/files/int_paper/paper/MC2017_2017_3/P197S03-09WalshJ.pdf>`_,"
  *Proc. Int. Conf. Mathematics & Computational Methods Applied to Nuclear
  Science and Engineering*, Jeju, Korea, Apr. 16-20, 2017.

@ -409,7 +414,7 @@ Parallelism

 - Paul K. Romano and Andrew R. Siegel, "`Limits on the efficiency of event-based
  algorithms for Monte Carlo neutron transport
-  <https://www.kns.org/paper_file/paper/MC2017_2017_2/P099S02-02RomanoP.pdf>`_,"
+  <https://www.kns.org/files/int_paper/paper/MC2017_2017_2/P099S02-02RomanoP.pdf>`_,"
  *Proc. Int. Conf. Mathematics & Computational Methods Applied to Nuclear
  Science and Engineering*, Jeju, Korea, Apr. 16-20, 2017.

--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@ -2,6 +2,8 @@
 :mod:`openmc.data` -- Nuclear Data Interface
 --------------------------------------------

+.. module:: openmc.data
+
 Core Classes
 ------------

@ -13,15 +15,15 @@ and product yields.
    :nosignatures:
    :template: myclass.rst

-    openmc.data.IncidentNeutron
-    openmc.data.Reaction
-    openmc.data.Product
-    openmc.data.FissionEnergyRelease
-    openmc.data.DataLibrary
-    openmc.data.Decay
-    openmc.data.FissionProductYields
-    openmc.data.WindowedMultipole
-    openmc.data.ProbabilityTables
+    IncidentNeutron
+    Reaction
+    Product
+    FissionEnergyRelease
+    DataLibrary
+    Decay
+    FissionProductYields
+    WindowedMultipole
+    ProbabilityTables

 The following classes are used for storing atomic data (incident photon cross
 sections, atomic relaxation):
@ -31,9 +33,9 @@ sections, atomic relaxation):
    :nosignatures:
    :template: myclass.rst

-    openmc.data.IncidentPhoton
-    openmc.data.PhotonReaction
-    openmc.data.AtomicRelaxation
+    IncidentPhoton
+    PhotonReaction
+    AtomicRelaxation


 The following classes are used for storing thermal neutron scattering data:
@ -43,10 +45,10 @@ The following classes are used for storing thermal neutron scattering data:
    :nosignatures:
    :template: myclass.rst

-    openmc.data.ThermalScattering
-    openmc.data.ThermalScatteringReaction
-    openmc.data.CoherentElastic
-    openmc.data.IncoherentElastic
+    ThermalScattering
+    ThermalScatteringReaction
+    CoherentElastic
+    IncoherentElastic


 Core Functions
@ -57,12 +59,12 @@ Core Functions
    :nosignatures:
    :template: myfunction.rst

-    openmc.data.atomic_mass
-    openmc.data.gnd_name
-    openmc.data.linearize
-    openmc.data.thin
-    openmc.data.water_density
-    openmc.data.zam
+    atomic_mass
+    gnd_name
+    linearize
+    thin
+    water_density
+    zam

 One-dimensional Functions
 -------------------------
@ -72,13 +74,13 @@ One-dimensional Functions
    :nosignatures:
    :template: myclass.rst

-    openmc.data.Function1D
-    openmc.data.Tabulated1D
-    openmc.data.Polynomial
-    openmc.data.Combination
-    openmc.data.Sum
-    openmc.data.Regions1D
-    openmc.data.ResonancesWithBackground
+    Function1D
+    Tabulated1D
+    Polynomial
+    Combination
+    Sum
+    Regions1D
+    ResonancesWithBackground

 Angle-Energy Distributions
 --------------------------
@ -88,27 +90,27 @@ Angle-Energy Distributions
    :nosignatures:
    :template: myclass.rst

-    openmc.data.AngleEnergy
-    openmc.data.KalbachMann
-    openmc.data.CorrelatedAngleEnergy
-    openmc.data.UncorrelatedAngleEnergy
-    openmc.data.NBodyPhaseSpace
-    openmc.data.LaboratoryAngleEnergy
-    openmc.data.AngleDistribution
-    openmc.data.EnergyDistribution
-    openmc.data.ArbitraryTabulated
-    openmc.data.GeneralEvaporation
-    openmc.data.MaxwellEnergy
-    openmc.data.Evaporation
-    openmc.data.WattEnergy
-    openmc.data.MadlandNix
-    openmc.data.DiscretePhoton
-    openmc.data.LevelInelastic
-    openmc.data.ContinuousTabular
-    openmc.data.CoherentElasticAE
-    openmc.data.IncoherentElasticAE
-    openmc.data.IncoherentElasticAEDiscrete
-    openmc.data.IncoherentInelasticAEDiscrete
+    AngleEnergy
+    KalbachMann
+    CorrelatedAngleEnergy
+    UncorrelatedAngleEnergy
+    NBodyPhaseSpace
+    LaboratoryAngleEnergy
+    AngleDistribution
+    EnergyDistribution
+    ArbitraryTabulated
+    GeneralEvaporation
+    MaxwellEnergy
+    Evaporation
+    WattEnergy
+    MadlandNix
+    DiscretePhoton
+    LevelInelastic
+    ContinuousTabular
+    CoherentElasticAE
+    IncoherentElasticAE
+    IncoherentElasticAEDiscrete
+    IncoherentInelasticAEDiscrete

 Resonance Data
 --------------
@ -118,20 +120,20 @@ Resonance Data
    :nosignatures:
    :template: myclass.rst

-    openmc.data.Resonances
-    openmc.data.ResonanceRange
-    openmc.data.SingleLevelBreitWigner
-    openmc.data.MultiLevelBreitWigner
-    openmc.data.ReichMoore
-    openmc.data.RMatrixLimited
-    openmc.data.ResonanceCovariances
-    openmc.data.ResonanceCovarianceRange
-    openmc.data.SingleLevelBreitWignerCovariance
-    openmc.data.MultiLevelBreitWignerCovariance
-    openmc.data.ReichMooreCovariance
-    openmc.data.ParticlePair
-    openmc.data.SpinGroup
-    openmc.data.Unresolved
+    Resonances
+    ResonanceRange
+    SingleLevelBreitWigner
+    MultiLevelBreitWigner
+    ReichMoore
+    RMatrixLimited
+    ResonanceCovariances
+    ResonanceCovarianceRange
+    SingleLevelBreitWignerCovariance
+    MultiLevelBreitWignerCovariance
+    ReichMooreCovariance
+    ParticlePair
+    SpinGroup
+    Unresolved

 ACE Format
 ----------
@ -144,8 +146,8 @@ Classes
    :nosignatures:
    :template: myclass.rst

-    openmc.data.ace.Library
-    openmc.data.ace.Table
+    ace.Library
+    ace.Table

 Functions
 +++++++++
@ -155,7 +157,7 @@ Functions
    :nosignatures:
    :template: myfunction.rst

-    openmc.data.ace.ascii_to_binary
+    ace.ascii_to_binary

 ENDF Format
 -----------
@ -168,7 +170,7 @@ Classes
    :nosignatures:
    :template: myclass.rst

-    openmc.data.endf.Evaluation
+    endf.Evaluation

 Functions
 +++++++++
@ -178,13 +180,13 @@ Functions
    :nosignatures:
    :template: myfunction.rst

-    openmc.data.endf.float_endf
-    openmc.data.endf.get_cont_record
-    openmc.data.endf.get_evaluations
-    openmc.data.endf.get_head_record
-    openmc.data.endf.get_tab1_record
-    openmc.data.endf.get_tab2_record
-    openmc.data.endf.get_text_record
+    endf.float_endf
+    endf.get_cont_record
+    endf.get_evaluations
+    endf.get_head_record
+    endf.get_tab1_record
+    endf.get_tab2_record
+    endf.get_text_record

 NJOY Interface
 --------------
@ -194,7 +196,7 @@ NJOY Interface
    :nosignatures:
    :template: myfunction.rst

-    openmc.data.njoy.run
-    openmc.data.njoy.make_pendf
-    openmc.data.njoy.make_ace
-    openmc.data.njoy.make_ace_thermal
+    njoy.run
+    njoy.make_pendf
+    njoy.make_ace
+    njoy.make_ace_thermal
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@ -1,11 +1,11 @@
 .. _pythonapi_deplete:

+.. module:: openmc.deplete
+
 ----------------------------------
 :mod:`openmc.deplete` -- Depletion
 ----------------------------------

-.. module:: openmc.deplete
-
 Primary API
 -----------

@ -82,10 +82,10 @@ A minimal example for performing depletion would be:
 Internal Classes and Functions
 ------------------------------

-When running in parallel using `mpi4py <http://mpi4py.scipy.org>`_, the MPI
-intercommunicator used can be changed by modifying the following module
-variable. If it is not explicitly modified, it defaults to
-``mpi4py.MPI.COMM_WORLD``.
+When running in parallel using `mpi4py
+<https://mpi4py.readthedocs.io/en/stable/>`_, the MPI intercommunicator used can
+be changed by modifying the following module variable. If it is not explicitly
+modified, it defaults to ``mpi4py.MPI.COMM_WORLD``.

 .. data:: comm

@ -126,8 +126,15 @@ data, such as number densities and reaction rates for each material.
   Results
   ResultsList

-The following functions are used to solve the depletion equations, with
-:func:`cram.CRAM48` being the default.
+The following class and functions are used to solve the depletion equations,
+with :func:`cram.CRAM48` being the default.
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myintegrator.rst
+
+   cram.IPFCramSolver

 .. autosummary::
   :toctree: generated
@ -161,7 +168,7 @@ Abstract Base Classes

 A good starting point for extending capabilities in :mod:`openmc.deplete` is
 to examine the following abstract base classes. Custom classes can
-inherit from :class:`abc.TransportOperator` to implement alternative 
+inherit from :class:`abc.TransportOperator` to implement alternative
 schemes for collecting reaction rates and other data from a transport code
 prior to depleting materials

@ -185,8 +192,8 @@ OpenMC simulations back on to the :class:`abc.TransportOperator`
   abc.ReactionRateHelper
   abc.TalliedFissionYieldHelper

-Custom integrators can be developed by subclassing from the following abstract
-base classes:
+Custom integrators or depletion solvers can be developed by subclassing from
+the following abstract base classes:

 .. autosummary::
   :toctree: generated
@ -195,3 +202,4 @@ base classes:

   abc.Integrator
   abc.SIIntegrator
+   abc.DepSystemSolver
--- a/docs/source/pythonapi/index.rst
+++ b/docs/source/pythonapi/index.rst
@ -27,7 +27,7 @@ there are many substantial benefits to using the Python API, including:

 For those new to Python, there are many good tutorials available online. We
 recommend going through the modules from `Codecademy
-<https://www.codecademy.com/tracks/python>`_ and/or the `Scipy lectures
+<https://www.codecademy.com/learn/learn-python-3>`_ and/or the `Scipy lectures
 <https://scipy-lectures.github.io/>`_.

 The full API documentation serves to provide more information on a given module
--- a/docs/source/pythonapi/stats.rst
+++ b/docs/source/pythonapi/stats.rst
@ -46,5 +46,6 @@ Spatial Distributions

   openmc.stats.Spatial
   openmc.stats.CartesianIndependent
+   openmc.stats.SphericalIndependent
   openmc.stats.Box
   openmc.stats.Point
--- a/docs/source/releasenotes/0.11.0.rst
+++ b/docs/source/releasenotes/0.11.0.rst
@ -2,11 +2,6 @@
 What's New in 0.11.0
 ====================

-.. note::
-    These release notes are for a future release of OpenMC and are still subject
-    to change. The new features and bug fixes documented here reflect the
-    current status of the ``develop`` branch of OpenMC.
-
 .. currentmodule:: openmc

 -------
@ -40,6 +35,7 @@ random, non-overlapping configuration of spheres within the region.
 New Features
 ------------

+- White boundary conditions can be applied to surfaces
 - Support for rectilinear meshes through :class:`openmc.RectilinearMesh`.
 - The :class:`Geometry`, :class:`Materials`, and :class:`Settings` classes now
  have a ``from_xml`` method that will build an instance from an existing XML
@ -47,7 +43,8 @@ New Features
 - Predefined energy group structures can be found in
  :data:`openmc.mgxs.GROUP_STRUCTURES`.
 - New tally scores: ``H1-production``, ``H2-production``, ``H3-production``,
-  ``He3-production``, ``He4-production``, ``heating``, and ``damage-energy``
+  ``He3-production``, ``He4-production``, ``heating``, ``heating-local``, and
+  ``damage-energy``.
 - Switched to cell-based neighor lists (`PR 1140
  <https://github.com/openmc-dev/openmc/pull/1140>`_)
 - Two new probability distributions that can be used for source distributions:
@ -88,6 +85,7 @@ Python API Changes
 Bug Fixes
 ---------

+- `Rotate azimuthal distributions correctly for source sampling <https://github.com/openmc-dev/openmc/pull/1363>`_
 - `Fix reading ASCII ACE tables in Python 3 <https://github.com/openmc-dev/openmc/pull/1176>`_
 - `Fix bug for distributed temperatures <https://github.com/openmc-dev/openmc/pull/1178>`_
 - `Fix bug for distance to boundary in complex cells <https://github.com/openmc-dev/openmc/pull/1172>`_
@ -106,6 +104,7 @@ This release contains new contributions from the following people:
 - `Brody Bassett <https://github.com/brbass>`_
 - `Will Boyd <https://github.com/wbinventor>`_
 - `Andrew Davis <https://github.com/makeclean>`_
+- `Iurii Drobyshev <https://github.com/dryuri92>`_
 - `Guillaume Giudicelli <https://github.com/GiudGiud>`_
 - `Brittany Grayson <https://github.com/graybri3>`_
 - `Zhuoran Han <https://github.com/hanzhuoran>`_
@ -123,9 +122,12 @@ This release contains new contributions from the following people:
 - `Isaac Meyer <https://github.com/icmeyer>`_
 - `April Novak <https://github.com/aprilnovak>`_
 - `Adam Nelson <https://github.com/nelsonag>`_
+- `Gavin Ridley <https://github.com/gridley>`_
 - `Jose Salcedo Perez <https://github.com/salcedop>`_
 - `Paul Romano <https://github.com/paulromano>`_
 - `Sam Shaner <https://github.com/samuelshaner>`_
 - `Jonathan Shimwell <https://github.com/Shimwell>`_
 - `Patrick Shriwise <https://github.com/pshriwise>`_
 - `John Tramm <https://github.com/jtramm>`_
+- `Jiankai Yu <https://github.com/rockfool>`_
+- `Xiaokang Zhang <https://github.com/zxkjack123>`_
--- a/docs/source/releasenotes/0.4.1.rst
+++ b/docs/source/releasenotes/0.4.1.rst
@ -45,9 +45,9 @@ Bug Fixes
 - `95cfac`_: Fixed error in cell neighbor searches.
 - `83a803`_: Fixed bug related to probability tables.

-.. _b206a8: https://github.com/mit-crpg/openmc/commit/b206a8
-.. _800742: https://github.com/mit-crpg/openmc/commit/800742
-.. _a07c08: https://github.com/mit-crpg/openmc/commit/a07c08
-.. _a75283: https://github.com/mit-crpg/openmc/commit/a75283
-.. _95cfac: https://github.com/mit-crpg/openmc/commit/95cfac
-.. _83a803: https://github.com/mit-crpg/openmc/commit/83a803
+.. _b206a8: https://github.com/openmc-dev/openmc/commit/b206a8
+.. _800742: https://github.com/openmc-dev/openmc/commit/800742
+.. _a07c08: https://github.com/openmc-dev/openmc/commit/a07c08
+.. _a75283: https://github.com/openmc-dev/openmc/commit/a75283
+.. _95cfac: https://github.com/openmc-dev/openmc/commit/95cfac
+.. _83a803: https://github.com/openmc-dev/openmc/commit/83a803
--- a/docs/source/releasenotes/0.4.2.rst
+++ b/docs/source/releasenotes/0.4.2.rst
@ -42,13 +42,13 @@ Bug Fixes
 - d050c7_: Added Bessel's correction to make estimate of variance unbiased.
 - 2a5b9c_: Fixed regression in plotting.

-.. _a27f8f: https://github.com/mit-crpg/openmc/commit/a27f8f
-.. _afe121: https://github.com/mit-crpg/openmc/commit/afe121
-.. _e0968e: https://github.com/mit-crpg/openmc/commit/e0968e
-.. _298db8: https://github.com/mit-crpg/openmc/commit/298db8
-.. _2f3bbe: https://github.com/mit-crpg/openmc/commit/2f3bbe
-.. _671f30: https://github.com/mit-crpg/openmc/commit/671f30
-.. _b2c40e: https://github.com/mit-crpg/openmc/commit/b2c40e
-.. _5524fd: https://github.com/mit-crpg/openmc/commit/5524fd
-.. _d050c7: https://github.com/mit-crpg/openmc/commit/d050c7
-.. _2a5b9c: https://github.com/mit-crpg/openmc/commit/2a5b9c
+.. _a27f8f: https://github.com/openmc-dev/openmc/commit/a27f8f
+.. _afe121: https://github.com/openmc-dev/openmc/commit/afe121
+.. _e0968e: https://github.com/openmc-dev/openmc/commit/e0968e
+.. _298db8: https://github.com/openmc-dev/openmc/commit/298db8
+.. _2f3bbe: https://github.com/openmc-dev/openmc/commit/2f3bbe
+.. _671f30: https://github.com/openmc-dev/openmc/commit/671f30
+.. _b2c40e: https://github.com/openmc-dev/openmc/commit/b2c40e
+.. _5524fd: https://github.com/openmc-dev/openmc/commit/5524fd
+.. _d050c7: https://github.com/openmc-dev/openmc/commit/d050c7
+.. _2a5b9c: https://github.com/openmc-dev/openmc/commit/2a5b9c
--- a/docs/source/releasenotes/0.4.3.rst
+++ b/docs/source/releasenotes/0.4.3.rst
@ -40,12 +40,12 @@ Bug Fixes
 - 3212f5_: Fixed issue with blank line at beginning of XML files.

 .. _nelsonag: https://github.com/nelsonag
-.. _33f29a: https://github.com/mit-crpg/openmc/commit/33f29a
-.. _1c472d: https://github.com/mit-crpg/openmc/commit/1c472d
-.. _3c6e80: https://github.com/mit-crpg/openmc/commit/3c6e80
-.. _3bd35b: https://github.com/mit-crpg/openmc/commit/3bd35b
-.. _0069d5: https://github.com/mit-crpg/openmc/commit/0069d5
-.. _7af2cf: https://github.com/mit-crpg/openmc/commit/7af2cf
-.. _460ef1: https://github.com/mit-crpg/openmc/commit/460ef1
-.. _85a60e: https://github.com/mit-crpg/openmc/commit/85a60e
-.. _3212f5: https://github.com/mit-crpg/openmc/commit/3212f5
+.. _33f29a: https://github.com/openmc-dev/openmc/commit/33f29a
+.. _1c472d: https://github.com/openmc-dev/openmc/commit/1c472d
+.. _3c6e80: https://github.com/openmc-dev/openmc/commit/3c6e80
+.. _3bd35b: https://github.com/openmc-dev/openmc/commit/3bd35b
+.. _0069d5: https://github.com/openmc-dev/openmc/commit/0069d5
+.. _7af2cf: https://github.com/openmc-dev/openmc/commit/7af2cf
+.. _460ef1: https://github.com/openmc-dev/openmc/commit/460ef1
+.. _85a60e: https://github.com/openmc-dev/openmc/commit/85a60e
+.. _3212f5: https://github.com/openmc-dev/openmc/commit/3212f5
--- a/docs/source/releasenotes/0.4.4.rst
+++ b/docs/source/releasenotes/0.4.4.rst
@ -36,8 +36,8 @@ Bug Fixes
 - 7fd617_: Fixed bug with restart runs in parallel.
 - dc4a8f_: Fixed bug with fixed source restart runs.

-.. _4654ee: https://github.com/mit-crpg/openmc/commit/4654ee
-.. _7ee461: https://github.com/mit-crpg/openmc/commit/7ee461
-.. _792eb3: https://github.com/mit-crpg/openmc/commit/792eb3
-.. _7fd617: https://github.com/mit-crpg/openmc/commit/7fd617
-.. _dc4a8f: https://github.com/mit-crpg/openmc/commit/dc4a8f
+.. _4654ee: https://github.com/openmc-dev/openmc/commit/4654ee
+.. _7ee461: https://github.com/openmc-dev/openmc/commit/7ee461
+.. _792eb3: https://github.com/openmc-dev/openmc/commit/792eb3
+.. _7fd617: https://github.com/openmc-dev/openmc/commit/7fd617
+.. _dc4a8f: https://github.com/openmc-dev/openmc/commit/dc4a8f
--- a/docs/source/releasenotes/0.5.0.rst
+++ b/docs/source/releasenotes/0.5.0.rst
@ -39,11 +39,11 @@ Bug Fixes
 - 6f8d9d_: Set default tally labels.
 - 6a3a5e_: Fix problem with corner-crossing in lattices.

-.. _737b90: https://github.com/mit-crpg/openmc/commit/737b90
-.. _a819b4: https://github.com/mit-crpg/openmc/commit/a819b4
-.. _b11696: https://github.com/mit-crpg/openmc/commit/b11696
-.. _2bd46a: https://github.com/mit-crpg/openmc/commit/2bd46a
-.. _7a1f08: https://github.com/mit-crpg/openmc/commit/7a1f08
-.. _c0e3ec: https://github.com/mit-crpg/openmc/commit/c0e3ec
-.. _6f8d9d: https://github.com/mit-crpg/openmc/commit/6f8d9d
-.. _6a3a5e: https://github.com/mit-crpg/openmc/commit/6a3a5e
+.. _737b90: https://github.com/openmc-dev/openmc/commit/737b90
+.. _a819b4: https://github.com/openmc-dev/openmc/commit/a819b4
+.. _b11696: https://github.com/openmc-dev/openmc/commit/b11696
+.. _2bd46a: https://github.com/openmc-dev/openmc/commit/2bd46a
+.. _7a1f08: https://github.com/openmc-dev/openmc/commit/7a1f08
+.. _c0e3ec: https://github.com/openmc-dev/openmc/commit/c0e3ec
+.. _6f8d9d: https://github.com/openmc-dev/openmc/commit/6f8d9d
+.. _6a3a5e: https://github.com/openmc-dev/openmc/commit/6a3a5e
--- a/docs/source/releasenotes/0.5.1.rst
+++ b/docs/source/releasenotes/0.5.1.rst
@ -36,8 +36,8 @@ Bug Fixes
 - 63bfd2_: Fix tracklength tallies with cell filter and universes.
 - 88daf7_: Fix analog tallies with survival biasing.

-.. _94103e: https://github.com/mit-crpg/openmc/commit/94103e
-.. _e77059: https://github.com/mit-crpg/openmc/commit/e77059
-.. _b0fe88: https://github.com/mit-crpg/openmc/commit/b0fe88
-.. _63bfd2: https://github.com/mit-crpg/openmc/commit/63bfd2
-.. _88daf7: https://github.com/mit-crpg/openmc/commit/88daf7
+.. _94103e: https://github.com/openmc-dev/openmc/commit/94103e
+.. _e77059: https://github.com/openmc-dev/openmc/commit/e77059
+.. _b0fe88: https://github.com/openmc-dev/openmc/commit/b0fe88
+.. _63bfd2: https://github.com/openmc-dev/openmc/commit/63bfd2
+.. _88daf7: https://github.com/openmc-dev/openmc/commit/88daf7
--- a/docs/source/releasenotes/0.5.2.rst
+++ b/docs/source/releasenotes/0.5.2.rst
@ -41,15 +41,15 @@ Bug Fixes
 - ab0793_: Corrected PETSC_NULL references to their correct types.
 - 182ebd_: Use analog estimator with energyout filter.

-.. _7632f3: https://github.com/mit-crpg/openmc/commit/7632f3
-.. _f85ac4: https://github.com/mit-crpg/openmc/commit/f85ac4
-.. _49c36b: https://github.com/mit-crpg/openmc/commit/49c36b
-.. _5ccc78: https://github.com/mit-crpg/openmc/commit/5ccc78
-.. _b1f52f: https://github.com/mit-crpg/openmc/commit/b1f52f
-.. _eae7e5: https://github.com/mit-crpg/openmc/commit/eae7e5
-.. _10c1cc: https://github.com/mit-crpg/openmc/commit/10c1cc
-.. _afdb50: https://github.com/mit-crpg/openmc/commit/afdb50
-.. _a3c593: https://github.com/mit-crpg/openmc/commit/a3c593
-.. _3a66e3: https://github.com/mit-crpg/openmc/commit/3a66e3
-.. _ab0793: https://github.com/mit-crpg/openmc/commit/ab0793
-.. _182ebd: https://github.com/mit-crpg/openmc/commit/182ebd
+.. _7632f3: https://github.com/openmc-dev/openmc/commit/7632f3
+.. _f85ac4: https://github.com/openmc-dev/openmc/commit/f85ac4
+.. _49c36b: https://github.com/openmc-dev/openmc/commit/49c36b
+.. _5ccc78: https://github.com/openmc-dev/openmc/commit/5ccc78
+.. _b1f52f: https://github.com/openmc-dev/openmc/commit/b1f52f
+.. _eae7e5: https://github.com/openmc-dev/openmc/commit/eae7e5
+.. _10c1cc: https://github.com/openmc-dev/openmc/commit/10c1cc
+.. _afdb50: https://github.com/openmc-dev/openmc/commit/afdb50
+.. _a3c593: https://github.com/openmc-dev/openmc/commit/a3c593
+.. _3a66e3: https://github.com/openmc-dev/openmc/commit/3a66e3
+.. _ab0793: https://github.com/openmc-dev/openmc/commit/ab0793
+.. _182ebd: https://github.com/openmc-dev/openmc/commit/182ebd
--- a/docs/source/releasenotes/0.5.3.rst
+++ b/docs/source/releasenotes/0.5.3.rst
@ -40,8 +40,8 @@ Bug Fixes
 - c18a6e_: Check for valid secondary mode on S(a,b) tables.
 - 82c456_: Fix bug where last process could have zero particles.

-.. _2b1e8a: https://github.com/mit-crpg/openmc/commit/2b1e8a
-.. _5853d2: https://github.com/mit-crpg/openmc/commit/5853d2
-.. _e178c7: https://github.com/mit-crpg/openmc/commit/e178c7
-.. _c18a6e: https://github.com/mit-crpg/openmc/commit/c18a6e
-.. _82c456: https://github.com/mit-crpg/openmc/commit/82c456
+.. _2b1e8a: https://github.com/openmc-dev/openmc/commit/2b1e8a
+.. _5853d2: https://github.com/openmc-dev/openmc/commit/5853d2
+.. _e178c7: https://github.com/openmc-dev/openmc/commit/e178c7
+.. _c18a6e: https://github.com/openmc-dev/openmc/commit/c18a6e
+.. _82c456: https://github.com/openmc-dev/openmc/commit/82c456
--- a/docs/source/releasenotes/0.5.4.rst
+++ b/docs/source/releasenotes/0.5.4.rst
@ -39,13 +39,13 @@ Bug Fixes
 - cf567c_: ENDF/B-VI data checked for compatibility
 - 6b9461_: Fix p_valid sampling inside of sample_energy

-.. _32c03c: https://github.com/mit-crpg/openmc/commit/32c03c
-.. _c71ef5: https://github.com/mit-crpg/openmc/commit/c71ef5
-.. _8884fb: https://github.com/mit-crpg/openmc/commit/8884fb
-.. _b38af0: https://github.com/mit-crpg/openmc/commit/b38af0
-.. _d28750: https://github.com/mit-crpg/openmc/commit/d28750
-.. _cf567c: https://github.com/mit-crpg/openmc/commit/cf567c
-.. _6b9461: https://github.com/mit-crpg/openmc/commit/6b9461
+.. _32c03c: https://github.com/openmc-dev/openmc/commit/32c03c
+.. _c71ef5: https://github.com/openmc-dev/openmc/commit/c71ef5
+.. _8884fb: https://github.com/openmc-dev/openmc/commit/8884fb
+.. _b38af0: https://github.com/openmc-dev/openmc/commit/b38af0
+.. _d28750: https://github.com/openmc-dev/openmc/commit/d28750
+.. _cf567c: https://github.com/openmc-dev/openmc/commit/cf567c
+.. _6b9461: https://github.com/openmc-dev/openmc/commit/6b9461

 ------------
 Contributors
--- a/docs/source/releasenotes/0.6.0.rst
+++ b/docs/source/releasenotes/0.6.0.rst
@ -35,13 +35,13 @@ Bug Fixes
 - d7a7d0_: Fix bug with <element> specifying xs attribute
 - 85b3cb_: Fix out-of-bounds error with OpenMP threading

-.. _41f7ca: https://github.com/mit-crpg/openmc/commit/41f7ca
-.. _038736: https://github.com/mit-crpg/openmc/commit/038736
-.. _46f9e8: https://github.com/mit-crpg/openmc/commit/46f9e8
-.. _d1ca35: https://github.com/mit-crpg/openmc/commit/d1ca35
-.. _0291c0: https://github.com/mit-crpg/openmc/commit/0291c0
-.. _d7a7d0: https://github.com/mit-crpg/openmc/commit/d7a7d0
-.. _85b3cb: https://github.com/mit-crpg/openmc/commit/85b3cb
+.. _41f7ca: https://github.com/openmc-dev/openmc/commit/41f7ca
+.. _038736: https://github.com/openmc-dev/openmc/commit/038736
+.. _46f9e8: https://github.com/openmc-dev/openmc/commit/46f9e8
+.. _d1ca35: https://github.com/openmc-dev/openmc/commit/d1ca35
+.. _0291c0: https://github.com/openmc-dev/openmc/commit/0291c0
+.. _d7a7d0: https://github.com/openmc-dev/openmc/commit/d7a7d0
+.. _85b3cb: https://github.com/openmc-dev/openmc/commit/85b3cb

 ------------
 Contributors
--- a/docs/source/releasenotes/0.6.1.rst
+++ b/docs/source/releasenotes/0.6.1.rst
@ -38,16 +38,16 @@ Bug Fixes
 - 2a95ef_: Prevent segmentation fault on "current" score without mesh filter
 - 93e482_: Check for negative values in probability tables

-.. _03e890: https://github.com/mit-crpg/openmc/commit/03e890
-.. _4439de: https://github.com/mit-crpg/openmc/commit/4439de
-.. _5808ed: https://github.com/mit-crpg/openmc/commit/5808ed
-.. _2e60c0: https://github.com/mit-crpg/openmc/commit/2e60c0
-.. _3e0870: https://github.com/mit-crpg/openmc/commit/3e0870
-.. _dc4776: https://github.com/mit-crpg/openmc/commit/dc4776
-.. _01178b: https://github.com/mit-crpg/openmc/commit/01178b
-.. _62ec43: https://github.com/mit-crpg/openmc/commit/62ec43
-.. _2a95ef: https://github.com/mit-crpg/openmc/commit/2a95ef
-.. _93e482: https://github.com/mit-crpg/openmc/commit/93e482
+.. _03e890: https://github.com/openmc-dev/openmc/commit/03e890
+.. _4439de: https://github.com/openmc-dev/openmc/commit/4439de
+.. _5808ed: https://github.com/openmc-dev/openmc/commit/5808ed
+.. _2e60c0: https://github.com/openmc-dev/openmc/commit/2e60c0
+.. _3e0870: https://github.com/openmc-dev/openmc/commit/3e0870
+.. _dc4776: https://github.com/openmc-dev/openmc/commit/dc4776
+.. _01178b: https://github.com/openmc-dev/openmc/commit/01178b
+.. _62ec43: https://github.com/openmc-dev/openmc/commit/62ec43
+.. _2a95ef: https://github.com/openmc-dev/openmc/commit/2a95ef
+.. _93e482: https://github.com/openmc-dev/openmc/commit/93e482

 ------------
 Contributors
--- a/docs/source/releasenotes/0.6.2.rst
+++ b/docs/source/releasenotes/0.6.2.rst
@ -33,10 +33,10 @@ Bug Fixes
 - e6abb9_: Fix segfault when tallying in a void material
 - 291b45_: Handle metastable nuclides in NNDC data and multiplicities in MT=5 data

-.. _26fb93: https://github.com/mit-crpg/openmc/commit/26fb93
-.. _2f07c0: https://github.com/mit-crpg/openmc/commit/2f07c0
-.. _e6abb9: https://github.com/mit-crpg/openmc/commit/e6abb9
-.. _291b45: https://github.com/mit-crpg/openmc/commit/291b45
+.. _26fb93: https://github.com/openmc-dev/openmc/commit/26fb93
+.. _2f07c0: https://github.com/openmc-dev/openmc/commit/2f07c0
+.. _e6abb9: https://github.com/openmc-dev/openmc/commit/e6abb9
+.. _291b45: https://github.com/openmc-dev/openmc/commit/291b45

 ------------
 Contributors
--- a/docs/source/releasenotes/0.7.0.rst
+++ b/docs/source/releasenotes/0.7.0.rst
@ -40,11 +40,11 @@ Bug Fixes
 - 6121d9_: Fix bugs related to particle track files
 - 2f0e89_: Fixes for nuclide specification in tallies

-.. _b5f712: https://github.com/mit-crpg/openmc/commit/b5f712
-.. _e6675b: https://github.com/mit-crpg/openmc/commit/e6675b
-.. _04e2c1: https://github.com/mit-crpg/openmc/commit/04e2c1
-.. _6121d9: https://github.com/mit-crpg/openmc/commit/6121d9
-.. _2f0e89: https://github.com/mit-crpg/openmc/commit/2f0e89
+.. _b5f712: https://github.com/openmc-dev/openmc/commit/b5f712
+.. _e6675b: https://github.com/openmc-dev/openmc/commit/e6675b
+.. _04e2c1: https://github.com/openmc-dev/openmc/commit/04e2c1
+.. _6121d9: https://github.com/openmc-dev/openmc/commit/6121d9
+.. _2f0e89: https://github.com/openmc-dev/openmc/commit/2f0e89

 ------------
 Contributors
--- a/docs/source/releasenotes/0.7.1.rst
+++ b/docs/source/releasenotes/0.7.1.rst
@ -62,15 +62,15 @@ Bug Fixes
 - 441fd4_: Fix bug in kappa-fission score
 - 7e5974_: Allow fixed source simulations from Python API

-.. _299322: https://github.com/mit-crpg/openmc/commit/299322
-.. _d74840: https://github.com/mit-crpg/openmc/commit/d74840
-.. _c29a81: https://github.com/mit-crpg/openmc/commit/c29a81
-.. _3edc23: https://github.com/mit-crpg/openmc/commit/3edc23
-.. _629e3b: https://github.com/mit-crpg/openmc/commit/629e3b
-.. _5dbe8b: https://github.com/mit-crpg/openmc/commit/5dbe8b
-.. _ff66f4: https://github.com/mit-crpg/openmc/commit/ff66f4
-.. _441fd4: https://github.com/mit-crpg/openmc/commit/441fd4
-.. _7e5974: https://github.com/mit-crpg/openmc/commit/7e5974
+.. _299322: https://github.com/openmc-dev/openmc/commit/299322
+.. _d74840: https://github.com/openmc-dev/openmc/commit/d74840
+.. _c29a81: https://github.com/openmc-dev/openmc/commit/c29a81
+.. _3edc23: https://github.com/openmc-dev/openmc/commit/3edc23
+.. _629e3b: https://github.com/openmc-dev/openmc/commit/629e3b
+.. _5dbe8b: https://github.com/openmc-dev/openmc/commit/5dbe8b
+.. _ff66f4: https://github.com/openmc-dev/openmc/commit/ff66f4
+.. _441fd4: https://github.com/openmc-dev/openmc/commit/441fd4
+.. _7e5974: https://github.com/openmc-dev/openmc/commit/7e5974

 ------------
 Contributors
--- a/docs/source/releasenotes/0.8.0.rst
+++ b/docs/source/releasenotes/0.8.0.rst
@ -65,17 +65,17 @@ Bug Fixes
 - 8467ae_: Better threshold for allowable lost particles
 - 493c6f_: Fix type of return argument for h5pget_driver_f

-.. _70daa7: https://github.com/mit-crpg/openmc/commit/70daa7
-.. _40b05f: https://github.com/mit-crpg/openmc/commit/40b05f
-.. _9586ed: https://github.com/mit-crpg/openmc/commit/9586ed
-.. _a855e8: https://github.com/mit-crpg/openmc/commit/a855e8
-.. _7294a1: https://github.com/mit-crpg/openmc/commit/7294a1
-.. _12f246: https://github.com/mit-crpg/openmc/commit/12f246
-.. _0227f4: https://github.com/mit-crpg/openmc/commit/0227f4
-.. _51deaa: https://github.com/mit-crpg/openmc/commit/51deaa
-.. _fed74b: https://github.com/mit-crpg/openmc/commit/fed74b
-.. _8467ae: https://github.com/mit-crpg/openmc/commit/8467ae
-.. _493c6f: https://github.com/mit-crpg/openmc/commit/493c6f
+.. _70daa7: https://github.com/openmc-dev/openmc/commit/70daa7
+.. _40b05f: https://github.com/openmc-dev/openmc/commit/40b05f
+.. _9586ed: https://github.com/openmc-dev/openmc/commit/9586ed
+.. _a855e8: https://github.com/openmc-dev/openmc/commit/a855e8
+.. _7294a1: https://github.com/openmc-dev/openmc/commit/7294a1
+.. _12f246: https://github.com/openmc-dev/openmc/commit/12f246
+.. _0227f4: https://github.com/openmc-dev/openmc/commit/0227f4
+.. _51deaa: https://github.com/openmc-dev/openmc/commit/51deaa
+.. _fed74b: https://github.com/openmc-dev/openmc/commit/fed74b
+.. _8467ae: https://github.com/openmc-dev/openmc/commit/8467ae
+.. _493c6f: https://github.com/openmc-dev/openmc/commit/493c6f

 ------------
 Contributors
--- a/docs/source/releasenotes/0.9.0.rst
+++ b/docs/source/releasenotes/0.9.0.rst
@ -111,25 +111,25 @@ Bug Fixes
 - 489540_: Check for void materials in tracklength tallies
 - f0214f_: Fixes/improvements to the ARES algorithm

-.. _c5df6c: https://github.com/mit-crpg/openmc/commit/c5df6c
-.. _1cfa39: https://github.com/mit-crpg/openmc/commit/1cfa39
-.. _335359: https://github.com/mit-crpg/openmc/commit/335359
-.. _17c678: https://github.com/mit-crpg/openmc/commit/17c678
-.. _23ec0b: https://github.com/mit-crpg/openmc/commit/23ec0b
-.. _7eefb7: https://github.com/mit-crpg/openmc/commit/7eefb7
-.. _7880d4: https://github.com/mit-crpg/openmc/commit/7880d4
-.. _ad2d9f: https://github.com/mit-crpg/openmc/commit/ad2d9f
-.. _59fdca: https://github.com/mit-crpg/openmc/commit/59fdca
-.. _9eff5b: https://github.com/mit-crpg/openmc/commit/9eff5b
-.. _7848a9: https://github.com/mit-crpg/openmc/commit/7848a9
-.. _f139ce: https://github.com/mit-crpg/openmc/commit/f139ce
-.. _b8ddfa: https://github.com/mit-crpg/openmc/commit/b8ddfa
-.. _ec3cfb: https://github.com/mit-crpg/openmc/commit/ec3cfb
-.. _5e9b06: https://github.com/mit-crpg/openmc/commit/5e9b06
-.. _c39990: https://github.com/mit-crpg/openmc/commit/c39990
-.. _c6b67e: https://github.com/mit-crpg/openmc/commit/c6b67e
-.. _489540: https://github.com/mit-crpg/openmc/commit/489540
-.. _f0214f: https://github.com/mit-crpg/openmc/commit/f0214f
+.. _c5df6c: https://github.com/openmc-dev/openmc/commit/c5df6c
+.. _1cfa39: https://github.com/openmc-dev/openmc/commit/1cfa39
+.. _335359: https://github.com/openmc-dev/openmc/commit/335359
+.. _17c678: https://github.com/openmc-dev/openmc/commit/17c678
+.. _23ec0b: https://github.com/openmc-dev/openmc/commit/23ec0b
+.. _7eefb7: https://github.com/openmc-dev/openmc/commit/7eefb7
+.. _7880d4: https://github.com/openmc-dev/openmc/commit/7880d4
+.. _ad2d9f: https://github.com/openmc-dev/openmc/commit/ad2d9f
+.. _59fdca: https://github.com/openmc-dev/openmc/commit/59fdca
+.. _9eff5b: https://github.com/openmc-dev/openmc/commit/9eff5b
+.. _7848a9: https://github.com/openmc-dev/openmc/commit/7848a9
+.. _f139ce: https://github.com/openmc-dev/openmc/commit/f139ce
+.. _b8ddfa: https://github.com/openmc-dev/openmc/commit/b8ddfa
+.. _ec3cfb: https://github.com/openmc-dev/openmc/commit/ec3cfb
+.. _5e9b06: https://github.com/openmc-dev/openmc/commit/5e9b06
+.. _c39990: https://github.com/openmc-dev/openmc/commit/c39990
+.. _c6b67e: https://github.com/openmc-dev/openmc/commit/c6b67e
+.. _489540: https://github.com/openmc-dev/openmc/commit/489540
+.. _f0214f: https://github.com/openmc-dev/openmc/commit/f0214f

 ------------
 Contributors
--- a/docs/source/usersguide/beginners.rst
+++ b/docs/source/usersguide/beginners.rst
@ -141,13 +141,13 @@ and `Volume II`_. You may also find it helpful to review the following terms:
 - `Effective multiplication factor`_
 - `Flux`_

-.. _nuclear reactor: http://en.wikipedia.org/wiki/Nuclear_reactor
-.. _Monte Carlo: http://en.wikipedia.org/wiki/Monte_Carlo_method
-.. _fission: http://en.wikipedia.org/wiki/Nuclear_fission
-.. _deterministic: http://en.wikipedia.org/wiki/Deterministic_algorithm
-.. _neutron transport: http://en.wikipedia.org/wiki/Neutron_transport
-.. _discretization: http://en.wikipedia.org/wiki/Discretization
-.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
+.. _nuclear reactor: https://en.wikipedia.org/wiki/Nuclear_reactor
+.. _Monte Carlo: https://en.wikipedia.org/wiki/Monte_Carlo_method
+.. _fission: https://en.wikipedia.org/wiki/Nuclear_fission
+.. _deterministic: https://en.wikipedia.org/wiki/Deterministic_algorithm
+.. _neutron transport: https://en.wikipedia.org/wiki/Neutron_transport
+.. _discretization: https://en.wikipedia.org/wiki/Discretization
+.. _constructive solid geometry: https://en.wikipedia.org/wiki/Constructive_solid_geometry
 .. _git: http://git-scm.com/
 .. _git tutorials: http://git-scm.com/documentation
 .. _Reactor Concepts Manual: http://www.tayloredge.com/periodic/trivia/ReactorConcepts.pdf
@ -156,6 +156,6 @@ and `Volume II`_. You may also find it helpful to review the following terms:
 .. _OpenMC source code: https://github.com/openmc-dev/openmc
 .. _GitHub: https://github.com/
 .. _bug reports: https://github.com/openmc-dev/openmc/issues
-.. _Neutron cross section: http://en.wikipedia.org/wiki/Neutron_cross_section
+.. _Neutron cross section: https://en.wikipedia.org/wiki/Neutron_cross_section
 .. _Effective multiplication factor: https://en.wikipedia.org/wiki/Nuclear_chain_reaction#Effective_neutron_multiplication_factor
-.. _Flux: http://en.wikipedia.org/wiki/Neutron_flux
+.. _Flux: https://en.wikipedia.org/wiki/Neutron_flux
--- a/docs/source/usersguide/cross_sections.rst
+++ b/docs/source/usersguide/cross_sections.rst
@ -16,10 +16,10 @@ recommended to use one of the pregenerated libraries. Alternatively, if you have
 ACE format data that was produced with NJOY_, such as that distributed with
 MCNP_ or Serpent_, it can be converted to the HDF5 format using the :ref:`using
 the Python API <create_xs_library>`. Several sources provide openly available
-ACE data including the `ENDF/B`_, JEFF_, and TENDL_
-libraries. In addition to tabulated cross sections in the HDF5 files, OpenMC
-relies on :ref:`windowed multipole <windowed_multipole>` data to perform
-on-the-fly Doppler broadening.
+ACE data including the `ENDF/B`_, JEFF_, and TENDL_ libraries as well as the
+`LANL Nuclear Data Team <https://nucleardata.lanl.gov/>`_. In addition to
+tabulated cross sections in the HDF5 files, OpenMC relies on :ref:`windowed
+multipole <windowed_multipole>` data to perform on-the-fly Doppler broadening.

 In multi-group mode, OpenMC utilizes an HDF5-based library format which can be
 used to describe nuclide- or material-specific quantities.
@ -30,11 +30,11 @@ Environment Variables

 When :ref:`scripts_openmc` is run, it will look for several environment
 variables that indicate where cross sections can be found. While the location of
-cross sections can also be indicated through the :class:`openmc.Materials` class
-(or in the :ref:`materials.xml <io_materials>` file), if you always use the same
-set of cross section data, it is often easier to just set an environment
-variable that will be picked up by default every time OpenMC is run. The
-following environment variables are used:
+cross sections can also be indicated through the
+:attr:`openmc.Materials.cross_setion` attribute (or in the :ref:`materials.xml
+<io_materials>` file), if you always use the same set of cross section data, it
+is often easier to just set an environment variable that will be picked up by
+default every time OpenMC is run. The following environment variables are used:

 :envvar:`OPENMC_CROSS_SECTIONS`
  Indicates the path to the :ref:`cross_sections.xml <io_cross_sections>`
--- a/docs/source/usersguide/depletion.rst
+++ b/docs/source/usersguide/depletion.rst
@ -127,7 +127,9 @@ of the same material as a unique material definition with::
 For our example problem, this would deplete fuel on the outer region of the problem
 with different reaction rates than those in the center. Materials will be depleted
 corresponding to their local neutron spectra, and have unique compositions at each
-transport step.
+transport step.  The volume of the original ``fuel_3`` material must represent
+the volume of **all** the ``fuel_3`` in the problem. When creating the unique
+materials, this volume will be equally distributed across all material instances.


 .. note::
--- a/docs/source/usersguide/geometry.rst
+++ b/docs/source/usersguide/geometry.rst
@ -121,11 +121,11 @@ For many regions, a bounding-box can be determined automatically::
 While a bounding box can be determined for regions involving half-spaces of
 spheres, cylinders, and axis-aligned planes, it generally cannot be determined
 if the region involves cones, non-axis-aligned planes, or other exotic
-second-order surfaces. For example, the :func:`openmc.get_hexagonal_prism`
+second-order surfaces. For example, the :func:`openmc.model.hexagonal_prism`
 function returns the interior region of a hexagonal prism; because it is bounded
 by a :class:`openmc.Plane`, trying to get its bounding box won't work::

-  >>> hex = openmc.get_hexagonal_prism()
+  >>> hex = openmc.model.hexagonal_prism()
  >>> hex.bounding_box
  (array([-0.8660254,       -inf,       -inf]),
   array([ 0.8660254,        inf,        inf]))
@ -374,7 +374,7 @@ code would work::
  hexlat.universes = [outer_ring, middle_ring, inner_ring]

 If you need to create a hexagonal boundary (composed of six planar surfaces) for
-a hexagonal lattice, :func:`openmc.get_hexagonal_prism` can be used.
+a hexagonal lattice, :func:`openmc.model.hexagonal_prism` can be used.

 .. _usersguide_geom_export:

@ -396,8 +396,15 @@ if needed, lattices, the last step is to create an instance of
   geom.root_universe = root_univ
   geom.export_to_xml()

-.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
-.. _quadratic surfaces: http://en.wikipedia.org/wiki/Quadric
+Note that it's not strictly required to manually create a root universe. You can
+also pass a list of cells to the :class:`openmc.Geometry` constructor and it
+will handle creating the unverse::
+
+   geom = openmc.Geometry([cell1, cell2, cell3])
+   geom.export_to_xml()
+
+.. _constructive solid geometry: https://en.wikipedia.org/wiki/Constructive_solid_geometry
+.. _quadratic surfaces: https://en.wikipedia.org/wiki/Quadric

 --------------------------
 Using CAD-based Geometry
--- a/docs/source/usersguide/install.rst
+++ b/docs/source/usersguide/install.rst
@ -17,7 +17,8 @@ system for installing multiple versions of software packages and their
 dependencies and switching easily between them. `conda-forge
 <https://conda-forge.github.io/>`_ is a community-led conda channel of
 installable packages. For instructions on installing conda, please consult their
-`documentation <http://conda.pydata.org/docs/install/quick.html>`_.
+`documentation
+<https://docs.conda.io/projects/conda/en/latest/user-guide/install/>`_.

 Once you have `conda` installed on your system, add the `conda-forge` channel to
 your configuration with:
@ -444,7 +445,7 @@ as for OpenMC.
 .. admonition:: Optional
   :class: note

-   `mpi4py <http://mpi4py.scipy.org/>`_
+   `mpi4py <https://mpi4py.readthedocs.io/en/stable/>`_
      mpi4py provides Python bindings to MPI for running distributed-memory
      parallel runs. This package is needed if you plan on running depletion
      simulations in parallel using MPI.
@ -483,9 +484,9 @@ Make sure to replace the last string on the second line with the path to the
 schemas.xml file in your own OpenMC source directory.

 .. _GNU Emacs: http://www.gnu.org/software/emacs/
-.. _validation: http://en.wikipedia.org/wiki/XML_validation
+.. _validation: https://en.wikipedia.org/wiki/XML_validation
 .. _RELAX NG: http://relaxng.org/
 .. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html
 .. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html
-.. _Conda: https://conda.io/docs/
+.. _Conda: https://docs.conda.io/en/latest/
 .. _pip: https://pip.pypa.io/en/stable/
--- a/docs/source/usersguide/plots.rst
+++ b/docs/source/usersguide/plots.rst
@ -66,7 +66,7 @@ particular cells/materials should be given colors of your choosing::
  }

 Note that colors can be given as RGB tuples or by a string indicating a valid
-`SVG color <https://www.w3.org/TR/SVG/types.html#ColorKeywords>`_.
+`SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.

 When you're done creating your :class:`openmc.Plot` instances, you need to then
 assign them to a :class:`openmc.Plots` collection and export it to XML::
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@ -221,26 +221,11 @@ selected::

  settings.electron_treatment = 'led'

-.. warning::
-   Currently, collision stopping powers used in the TTB approximation come from
-   the `NIST ESTAR database`_, which provides data for each element calculated
-   using by default the material density at standard temperature and pressure.
-   In OpenMC, stopping powers for compounds are calculated from this elemental
-   data using Bragg's additivity rule. However, this is not a good
-   approximation --- the collision stopping power is a function of certain
-   quantities, such as the mean excitation energy and particularly the density
-   effect correction, that depend on material properties. Data for constituent
-   elements in a compound cannot simply be summed together, but rather these
-   quantities should be calculated for the material. This treatment will be
-   especially poor when the density of a material is different from the
-   densities used in the NIST data.
-
 .. note::
   Some features related to photon transport are not currently implemented,
   including:

     * Tallying photon energy deposition.
-     * Properly accounting for energy deposition in coupled n-p calculations.
     * Generating a photon source from a neutron calculation that can be used
       for a later fixed source photon calculation.
     * Photoneutron reactions.
--- a/docs/source/usersguide/volume.rst
+++ b/docs/source/usersguide/volume.rst
@ -37,6 +37,16 @@ arguments are not necessary. For example,
 Of course, the volumes that you *need* this capability for are often the ones
 with complex definitions.

+A threshold can be applied for the calculation's variance, standard deviation,
+or relative error of volume estimates using :meth:`openmc.VolumeCalculation.set_trigger`::
+
+    vol_calc.set_trigger(1e-05, 'std_dev')
+
+If a threshold is provided, calculations will be performed iteratively using the
+number of samples specified on the calculation until all volume uncertainties are below
+the threshold value. If no threshold is provided, the calculation will run the number of
+samples specified once and return the result.
+
 Once you have one or more :class:`openmc.VolumeCalculation` objects created, you
 can then assign then to :attr:`Settings.volume_calculations`::

--- a/examples/jupyter/chain_simple.xml
+++ b/examples/jupyter/chain_simple.xml
@ -0,0 +1 @@
+../../tests/chain_simple.xml
--- a/examples/jupyter/mg-mode-part-iii.ipynb
+++ b/examples/jupyter/mg-mode-part-iii.ipynb
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@ -110,6 +110,7 @@ extern "C" {
  int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
  int openmc_tally_get_scores(int32_t index, int** scores, int* n);
  int openmc_tally_get_type(int32_t index, int32_t* type);
+  int openmc_tally_get_writable(int32_t index, bool* writable);
  int openmc_tally_reset(int32_t index);
  int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]);
  int openmc_tally_set_active(int32_t index, bool active);
@ -119,6 +120,7 @@ extern "C" {
  int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
  int openmc_tally_set_scores(int32_t index, int n, const char** scores);
  int openmc_tally_set_type(int32_t index, const char* type);
+  int openmc_tally_set_writable(int32_t index, bool writable);
  int openmc_zernike_filter_get_order(int32_t index, int* order);
  int openmc_zernike_filter_get_params(int32_t index, double* x, double* y, double* r);
  int openmc_zernike_filter_set_order(int32_t index, int order);
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@ -179,10 +179,10 @@ public:

  //! \brief Rotational tranfsormation of the filled universe.
  //
-  //! The vector is empty if there is no rotation.  Otherwise, the first three
-  //! values are the rotation angles respectively about the x-, y-, and z-, axes
-  //! in degrees.  The next 9 values give the rotation matrix in row-major
-  //! order.
+  //! The vector is empty if there is no rotation. Otherwise, the first 9 values
+  //! give the rotation matrix in row-major order. When the user specifies
+  //! rotation angles about the x-, y- and z- axes in degrees, these values are
+  //! also present at the end of the vector, making it of length 12.
  std::vector<double> rotation_;

  std::vector<int32_t> offset_;  //!< Distribcell offset table
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@ -20,7 +20,7 @@ using double_4dvec = std::vector<std::vector<std::vector<std::vector<double>>>>;

 // OpenMC major, minor, and release numbers
 constexpr int VERSION_MAJOR {0};
-constexpr int VERSION_MINOR {11};
+constexpr int VERSION_MINOR {12};
 constexpr int VERSION_RELEASE {0};
 constexpr bool VERSION_DEV {true};
 constexpr std::array<int, 3> VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE};
--- a/include/openmc/distribution_spatial.h
+++ b/include/openmc/distribution_spatial.h
@ -37,6 +37,24 @@ private:
  UPtrDist z_; //!< Distribution of z coordinates
 };

+//==============================================================================
+//! Distribution of points specified by spherical coordinates r,theta,phi
+//==============================================================================
+
+class SphericalIndependent : public SpatialDistribution {
+public:
+  explicit SphericalIndependent(pugi::xml_node node);
+
+  //! Sample a position from the distribution
+  //! \return Sampled position
+  Position sample() const;
+private:
+  UPtrDist r_; //!< Distribution of r coordinates
+  UPtrDist theta_; //!< Distribution of theta coordinates
+  UPtrDist phi_; //!< Distribution of phi coordinates
+  Position origin_; //!< Cartesian coordinates of the sphere center
+};
+
 //==============================================================================
 //! Uniform distribution of points over a box
 //==============================================================================
--- a/include/openmc/geometry_aux.h
+++ b/include/openmc/geometry_aux.h
@ -7,9 +7,15 @@
 #include <cstdint>
 #include <string>
 #include <vector>
+#include <unordered_map>

 namespace openmc {

+namespace model {
+  extern std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>> universe_cell_counts;
+  extern std::unordered_map<int32_t, int32_t> universe_level_counts;
+} // namespace model
+
 void read_geometry_xml();

 //==============================================================================
--- a/include/openmc/hdf5_interface.h
+++ b/include/openmc/hdf5_interface.h
@ -345,7 +345,7 @@ read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false)
 }

 template <typename T, std::size_t N>
-void read_dataset_as_shape(hid_t obj_id, const char* name,
+inline void read_dataset_as_shape(hid_t obj_id, const char* name,
                           xt::xtensor<T, N>& arr, bool indep=false)
 {
  hid_t dset = open_dataset(obj_id, name);
@ -367,7 +367,7 @@ void read_dataset_as_shape(hid_t obj_id, const char* name,


 template <typename T, std::size_t N>
-void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
+inline void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
     bool must_have=false)
 {
  if (object_exists(obj_id, name)) {
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@ -11,6 +11,7 @@

 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/particle.h"
 #include "openmc/xsdata.h"


@ -110,7 +111,10 @@ class Mgxs {
    //!
    //! @param xs_id HDF5 group id for the cross section data.
    //! @param temperature Temperatures to read.
-    Mgxs(hid_t xs_id, const std::vector<double>& temperature);
+    //! @param num_group number of energy groups
+    //! @param num_delay number of delayed groups
+    Mgxs(hid_t xs_id, const std::vector<double>& temperature,
+        int num_group, int num_delay);

    //! \brief Constructor that initializes and populates all data to build a
    //!   macroscopic cross section from microscopic cross section.
@ -119,8 +123,11 @@ class Mgxs {
    //! @param mat_kTs temperatures (in units of eV) that data is needed.
    //! @param micros Microscopic objects to combine.
    //! @param atom_densities Atom densities of those microscopic quantities.
+    //! @param num_group number of energy groups
+    //! @param num_delay number of delayed groups
    Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
-         const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities);
+         const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
+         int num_group, int num_delay);

    //! \brief Provides a cross section value given certain parameters
    //!
@ -137,6 +144,11 @@ class Mgxs {
    get_xs(int xstype, int gin, const int* gout, const double* mu,
      const int* dg);

+    inline double
+    get_xs(int xstype, int gin)
+    {return get_xs(xstype, gin, nullptr, nullptr, nullptr);}
+
+
    //! \brief Samples the fission neutron energy and if prompt or delayed.
    //!
    //! @param gin Incoming energy group.
@ -156,15 +168,9 @@ class Mgxs {

    //! \brief Calculates cross section quantities needed for tracking.
    //!
-    //! @param gin Incoming energy group.
-    //! @param sqrtkT Temperature of the material.
-    //! @param u Incoming particle direction.
-    //! @param total_xs Resultant total cross section.
-    //! @param abs_xs Resultant absorption cross section.
-    //! @param nu_fiss_xs Resultant nu-fission cross section.
+    //! @param p The particle whose attributes set which MGXS to get.
    void
-    calculate_xs(int gin, double sqrtkT, Direction u,
-         double& total_xs, double& abs_xs, double& nu_fiss_xs);
+    calculate_xs(Particle& p);

    //! \brief Sets the temperature index in cache given a temperature
    //!
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@ -12,70 +12,71 @@
 namespace openmc {

 //==============================================================================
-// Global variables
+// Global MGXS data container structure
 //==============================================================================

+class MgxsInterface {
+public:
+
+  MgxsInterface() = default;
+
+  // Construct from path to cross sections file, as well as a list
+  // of XS to read and the corresponding temperatures for each XS
+  MgxsInterface(const std::string& path_cross_sections,
+                const std::vector<std::string> xs_to_read,
+                const std::vector<std::vector<double>> xs_temps);
+
+  // Does things to construct after the nuclides and temperatures to
+  // read have been specified.
+  void init();
+
+  // Set which nuclides and temperatures are to be read
+  void set_nuclides_and_temperatures(std::vector<std::string> xs_to_read,
+                                     std::vector<std::vector<double>> xs_temps);
+
+  // Add an Mgxs object to be managed
+  void add_mgxs(hid_t file_id, const std::string& name,
+         const std::vector<double>& temperature);
+
+  // Reads just the header of the cross sections file, to find
+  // min & max energies as well as the available XS
+  void read_header(const std::string& path_cross_sections);
+
+  // Calculate microscopic cross sections from nuclide macro XS
+  void create_macro_xs();
+
+  // Get the kT values which are used in the OpenMC model
+  std::vector<std::vector<double>> get_mat_kTs();
+
+  int num_energy_groups_;
+  int num_delayed_groups_;
+  std::vector<std::string> xs_names_; // available names in HDF5 file
+  std::vector<std::string> xs_to_read_; // XS which appear in materials
+  std::vector<std::vector<double>> xs_temps_to_read_; // temperatures used
+  std::string cross_sections_path_; // path to MGXS h5 file
+  std::vector<Mgxs> nuclides_;
+  std::vector<Mgxs> macro_xs_;
+  std::vector<double> energy_bins_;
+  std::vector<double> energy_bin_avg_;
+  std::vector<double> rev_energy_bins_;
+  std::vector<std::vector<double>> nuc_temps_; // all available temperatures
+};
+
 namespace data {
+  extern MgxsInterface mg;
+}

-extern std::vector<Mgxs> nuclides_MG;
-extern std::vector<Mgxs> macro_xs;
-extern int num_energy_groups;
-extern int num_delayed_groups;
-extern std::vector<double> energy_bins;
-extern std::vector<double> energy_bin_avg;
-extern std::vector<double> rev_energy_bins;
+// Puts available XS in MGXS file to globals so that when
+// materials are read, the MGXS specified in a material can
+// be ensured to be present in the available data.
+void put_mgxs_header_data_to_globals();

-} // namespace data
+// Set which nuclides and temperatures are to be read on
+// mg through global data
+void set_mg_interface_nuclides_and_temps();

-//==============================================================================
-// Mgxs data loading interface methods
-//==============================================================================
-
-void read_mgxs();
-
-void
-add_mgxs(hid_t file_id, const std::string& name,
-     const std::vector<double>& temperature);
-
-void create_macro_xs();
-
-std::vector<std::vector<double>> get_mat_kTs();
-
-void read_mg_cross_sections_header();
-
-//==============================================================================
-// Mgxs tracking/transport/tallying interface methods
-//==============================================================================
-
-extern "C" void
-calculate_xs_c(int i_mat, int gin, double sqrtkT, Direction u,
-     double& total_xs, double& abs_xs, double& nu_fiss_xs);
-
-double
-get_nuclide_xs(int index, int xstype, int gin, const int* gout,
-  const double* mu, const int* dg);
-
-inline double
-get_nuclide_xs(int index, int xstype, int gin)
-{return get_nuclide_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
-
-double
-get_macro_xs(int index, int xstype, int gin, const int* gout,
-  const double* mu, const int* dg);
-
-inline double
-get_macro_xs(int index, int xstype, int gin)
-{return get_macro_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
-
-//==============================================================================
-// General Mgxs methods
-//==============================================================================
-
-extern "C" void
-get_name_c(int index, int name_len, char* name);
-
-extern "C" double
-get_awr_c(int index);
+// After macro XS have been read, materials can be marked as fissionable
+void mark_fissionable_mgxs_materials();

 } // namespace openmc
 #endif // OPENMC_MGXS_INTERFACE_H
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@ -31,7 +31,7 @@ extern "C" bool cmfd_run;             //!< is a CMFD run?
 extern "C" bool dagmc;                //!< indicator of DAGMC geometry
 extern "C" bool entropy_on;           //!< calculate Shannon entropy?
 extern bool legendre_to_tabular;      //!< convert Legendre distributions to tabular?
-extern bool output_summary;           //!< write summary.h5?
+extern "C" bool output_summary;       //!< write summary.h5?
 extern bool output_tallies;           //!< write tallies.out?
 extern bool particle_restart_run;     //!< particle restart run?
 extern "C" bool photon_transport;     //!< photon transport turned on?
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@ -37,6 +37,8 @@ public:

  void set_active(bool active) { active_ = active; }

+  void set_writable(bool writable) { writable_ = writable; }
+
  void set_scores(pugi::xml_node node);

  void set_scores(const std::vector<std::string>& scores);
@ -55,6 +57,8 @@ public:

  int32_t n_filter_bins() const {return n_filter_bins_;}

+  bool writable() const { return writable_;}
+
  //----------------------------------------------------------------------------
  // Other methods.

@ -98,6 +102,9 @@ public:
  //! (e.g. specific cell, specific energy group, etc.)
  xt::xtensor<double, 3> results_;

+  //! True if this tally should be written to statepoint files
+  bool writable_ {true};
+
  //----------------------------------------------------------------------------
  // Miscellaneous public members.

--- a/include/openmc/volume_calc.h
+++ b/include/openmc/volume_calc.h
@ -2,12 +2,14 @@
 #define OPENMC_VOLUME_CALC_H

 #include "openmc/position.h"
+#include "openmc/tallies/trigger.h"

 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"

 #include <string>
 #include <vector>
+#include <gsl/gsl>

 namespace openmc {

@ -16,6 +18,7 @@ namespace openmc {
 //==============================================================================

 class VolumeCalculation {
+
 public:
  // Aliases, types
  struct Result {
@ -23,6 +26,7 @@ public:
    std::vector<int> nuclides; //!< Index of nuclides
    std::vector<double> atoms; //!< Number of atoms for each nuclide
    std::vector<double> uncertainty; //!< Uncertainty on number of atoms
+    int iterations; //!< Number of iterations needed to obtain the results
  }; // Results for a single domain

  // Constructors
@ -44,7 +48,9 @@ public:

  // Data members
  int domain_type_; //!< Type of domain (cell, material, etc.)
-  int n_samples_; //!< Number of samples to use
+  size_t n_samples_; //!< Number of samples to use
+  double threshold_ {-1.0}; //!< Error threshold for domain volumes
+  TriggerMetric trigger_type_ {TriggerMetric::not_active}; //!< Trigger metric for the volume calculation
  Position lower_left_; //!< Lower-left position of bounding box
  Position upper_right_; //!< Upper-right position of bounding box
  std::vector<int> domain_ids_; //!< IDs of domains to find volumes of
--- a/include/openmc/xsdata.h
+++ b/include/openmc/xsdata.h
@ -22,6 +22,7 @@ namespace openmc {
 class XsData {

  private:
+
    //! \brief Reads scattering data from the HDF5 file
    void
    scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
@ -61,6 +62,9 @@ class XsData {
    void
    fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang);

+    //! Number of energy and delayed neutron groups
+    size_t n_g_, n_dg_;
+
  public:

    // The following quantities have the following dimensions:
@ -98,7 +102,10 @@ class XsData {
    //! @param scatter_format The scattering representation of the file.
    //! @param n_pol Number of polar angles.
    //! @param n_azi Number of azimuthal angles.
-    XsData(bool fissionable, int scatter_format, int n_pol, int n_azi);
+    //! @param n_groups Number of energy groups.
+    //! @param n_d_groups Number of delayed neutron groups.
+    XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
+           size_t n_groups, size_t n_d_groups);

    //! \brief Loads the XsData object from the HDF5 file
    //!
--- a/openmc/init.py
+++ b/openmc/init.py
@ -34,4 +34,4 @@ from . import examples
 # Import a few convencience functions that used to be here
 from openmc.model import rectangular_prism, hexagonal_prism

-__version__ = '0.11.0-dev'
+__version__ = '0.12.0-dev'
--- a/openmc/cell.py
+++ b/openmc/cell.py
@ -63,6 +63,10 @@ class Cell(IDManagerMixin):
           \sin\phi \sin\theta \sin\psi & -\sin\phi \cos\psi + \cos\phi
           \sin\theta \sin\psi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi
           \cos\theta \end{array} \right ]
+
+        A rotation matrix can also be specified directly by setting this
+        attribute to a nested list (or 2D numpy array) that specifies each
+        element of the matrix.
    rotation_matrix : numpy.ndarray
        The rotation matrix defined by the angles specified in the
        :attr:`Cell.rotation` property.
@ -227,21 +231,24 @@ class Cell(IDManagerMixin):

    @rotation.setter
    def rotation(self, rotation):
-        cv.check_type('cell rotation', rotation, Iterable, Real)
        cv.check_length('cell rotation', rotation, 3)
        self._rotation = np.asarray(rotation)

        # Save rotation matrix -- the reason we do this instead of having it be
        # automatically calculated when the rotation_matrix property is accessed
        # is so that plotting on a rotated geometry can be done faster.
-        phi, theta, psi = self.rotation*(-pi/180.)
-        c3, s3 = cos(phi), sin(phi)
-        c2, s2 = cos(theta), sin(theta)
-        c1, s1 = cos(psi), sin(psi)
-        self._rotation_matrix = np.array([
-            [c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2],
-            [c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3],
-            [-s2, c2*s3, c2*c3]])
+        if self._rotation.ndim == 2:
+            # User specified rotation matrix directly
+            self._rotation_matrix = self._rotation
+        else:
+            phi, theta, psi = self.rotation*(-pi/180.)
+            c3, s3 = cos(phi), sin(phi)
+            c2, s2 = cos(theta), sin(theta)
+            c1, s1 = cos(psi), sin(psi)
+            self._rotation_matrix = np.array([
+                [c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2],
+                [c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3],
+                [-s2, c2*s3, c2*c3]])

    @translation.setter
    def translation(self, translation):
@ -343,7 +350,7 @@ class Cell(IDManagerMixin):

        return nuclides

-    def get_all_cells(self):
+    def get_all_cells(self, memo=None):
        """Return all cells that are contained within this one if it is filled with a
        universe or lattice

@ -357,12 +364,18 @@ class Cell(IDManagerMixin):

        cells = OrderedDict()

+        if memo and self in memo:
+            return cells
+
+        if memo is not None:
+            memo.add(self)
+
        if self.fill_type in ('universe', 'lattice'):
-            cells.update(self.fill.get_all_cells())
+            cells.update(self.fill.get_all_cells(memo))

        return cells

-    def get_all_materials(self):
+    def get_all_materials(self, memo=None):
        """Return all materials that are contained within the cell

        Returns
@ -381,9 +394,9 @@ class Cell(IDManagerMixin):
                    materials[m.id] = m
        else:
            # Append all Cells in each Cell in the Universe to the dictionary
-            cells = self.get_all_cells()
+            cells = self.get_all_cells(memo)
            for cell in cells.values():
-                materials.update(cell.get_all_materials())
+                materials.update(cell.get_all_materials(memo))

        return materials

@ -456,7 +469,24 @@ class Cell(IDManagerMixin):

        return memo[self]

-    def create_xml_subelement(self, xml_element):
+    def create_xml_subelement(self, xml_element, memo=None):
+        """Add the cell's xml representation to an incoming xml element
+
+        Parameters
+        ----------
+        xml_element : xml.etree.ElementTree.Element
+            XML element to be added to
+
+        memo : set or None
+            A set of object IDs representing geometry entities already
+            written to ``xml_element``. This parameter is used internally
+            and should not be specified by users.
+
+        Returns
+        -------
+        None
+
+        """
        element = ET.Element("cell")
        element.set("id", str(self.id))

@ -475,7 +505,7 @@ class Cell(IDManagerMixin):

        elif self.fill_type in ('universe', 'lattice'):
            element.set("fill", str(self.fill.id))
-            self.fill.create_xml_subelement(xml_element)
+            self.fill.create_xml_subelement(xml_element, memo)

        if self.region is not None:
            # Set the region attribute with the region specification
@ -491,19 +521,22 @@ class Cell(IDManagerMixin):
            # tree. When it reaches a leaf (a Halfspace), it creates a <surface>
            # element for the corresponding surface if none has been created
            # thus far.
-            def create_surface_elements(node, element):
+            def create_surface_elements(node, element, memo=None):
                if isinstance(node, Halfspace):
-                    path = "./surface[@id='{}']".format(node.surface.id)
-                    if xml_element.find(path) is None:
-                        xml_element.append(node.surface.to_xml_element())
+                    if memo and node.surface in memo:
+                        return
+                    if memo is not None:
+                        memo.add(node.surface)
+                    xml_element.append(node.surface.to_xml_element())
+
                elif isinstance(node, Complement):
-                    create_surface_elements(node.node, element)
+                    create_surface_elements(node.node, element, memo)
                else:
                    for subnode in node:
-                        create_surface_elements(subnode, element)
+                        create_surface_elements(subnode, element, memo)

            # Call the recursive function from the top node
-            create_surface_elements(self.region, xml_element)
+            create_surface_elements(self.region, xml_element, memo)

        if self.temperature is not None:
            if isinstance(self.temperature, Iterable):
@ -516,7 +549,7 @@ class Cell(IDManagerMixin):
            element.set("translation", ' '.join(map(str, self.translation)))

        if self.rotation is not None:
-            element.set("rotation", ' '.join(map(str, self.rotation)))
+            element.set("rotation", ' '.join(map(str, self.rotation.ravel())))

        return element

--- a/openmc/data/neutron.py
+++ b/openmc/data/neutron.py
@ -118,7 +118,12 @@ class IncidentNeutron(EqualityMixin):
        if mt in self.reactions:
            return self.reactions[mt]
        else:
-            raise KeyError('No reaction with MT={}.'.format(mt))
+            # Try to create a redundant cross section
+            mts = self.get_reaction_components(mt)
+            if len(mts) > 0:
+                return self._get_redundant_reaction(mt, mts)
+            else:
+                raise KeyError('No reaction with MT={}.'.format(mt))

    def __repr__(self):
        return "<IncidentNeutron: {}>".format(self.name)
@ -588,6 +593,9 @@ class IncidentNeutron(EqualityMixin):

        # If mass number hasn't been specified, make an educated guess
        zaid, xs = ace.name.split('.')
+        if not xs.endswith('c'):
+            raise TypeError(
+                "{} is not a continuous-energy neutron ACE table.".format(ace))
        name, element, Z, mass_number, metastable = \
            get_metadata(int(zaid), metastable_scheme)

@ -871,7 +879,7 @@ class IncidentNeutron(EqualityMixin):
                    fission = data.reactions[18].xs[temp]
                    kerma_fission = get_file3_xs(ev, 318, E)
                    kerma.y = kerma.y - kerma_fission + (
-                        f.fragments(E) + f.betas(E)) * fission.y
+                        f.fragments(E) + f.betas(E)) * fission(E)

                # For local KERMA, we first need to get the values from the
                # HEATR run with photon energy deposited locally and put
@ -884,7 +892,7 @@ class IncidentNeutron(EqualityMixin):
                    kerma_fission_local = get_file3_xs(ev_local, 318, E)
                    kerma_local = kerma_local - kerma_fission_local + (
                        f.fragments(E) + f.prompt_photons(E)
-                        + f.delayed_photons(E) + f.betas(E))*fission.y
+                        + f.delayed_photons(E) + f.betas(E))*fission(E)

                heating_local.xs[temp] = Tabulated1D(E, kerma_local)

@ -908,16 +916,17 @@ class IncidentNeutron(EqualityMixin):
            Redundant reaction

        """
-        # Get energy grid
-        strT = self.temperatures[0]
-        energy = self.energy[strT]

        rx = Reaction(mt)
-        xss = [self.reactions[mt_i].xs[strT] for mt_i in mts]
-        idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
-                   else 0 for xs in xss])
-        rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
-        rx.xs[strT]._threshold_idx = idx
+        # Get energy grid
+        for strT in self.temperatures:
+            energy = self.energy[strT]
+            xss = [self.reactions[mt_i].xs[strT] for mt_i in mts]
+            idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
+                       else 0 for xs in xss])
+            rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
+            rx.xs[strT]._threshold_idx = idx
+
        rx.redundant = True

        return rx
--- a/openmc/data/njoy.py
+++ b/openmc/data/njoy.py
@ -10,7 +10,7 @@ from . import endf


 # For a given MAT number, give a name for the ACE table and a list of ZAID
-# identifiers
+# identifiers. This is based on Appendix C in the ENDF manual.
 ThermalTuple = namedtuple('ThermalTuple', ['name', 'zaids', 'nmix'])
 _THERMAL_DATA = {
    1: ThermalTuple('hh2o', [1001], 1),
@ -23,31 +23,77 @@ _THERMAL_DATA = {
    11: ThermalTuple('dd2o', [1002], 1),
    12: ThermalTuple('parad', [1002], 1),
    13: ThermalTuple('orthod', [1002], 1),
+    14: ThermalTuple('dice', [1002], 1),
    26: ThermalTuple('be', [4009], 1),
    27: ThermalTuple('bebeo', [4009], 1),
-    31: ThermalTuple('graph', [6000, 6012, 6013], 1),
+    28: ThermalTuple('bebe2c', [4009], 1),
+    30: ThermalTuple('graph', [6000, 6012, 6013], 1),
+    31: ThermalTuple('grph10', [6000, 6012, 6013], 1),
+    32: ThermalTuple('grph30', [6000, 6012, 6013], 1),
    33: ThermalTuple('lch4', [1001], 1),
    34: ThermalTuple('sch4', [1001], 1),
+    35: ThermalTuple('sch4p2', [1001], 1),
    37: ThermalTuple('hch2', [1001], 1),
+    38: ThermalTuple('mesi00', [1001], 1),
    39: ThermalTuple('lucite', [1001], 1),
    40: ThermalTuple('benz', [1001, 6000, 6012], 2),
-    41: ThermalTuple('od2o', [8016, 8017, 8018], 1),
+    42: ThermalTuple('tol00', [1001], 1),
    43: ThermalTuple('sisic', [14028, 14029, 14030], 1),
    44: ThermalTuple('csic', [6000, 6012, 6013], 1),
+    45: ThermalTuple('ouo2', [8016, 8017, 8018], 1),
    46: ThermalTuple('obeo', [8016, 8017, 8018], 1),
    47: ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3),
-    48: ThermalTuple('uuo2', [92238], 1),
+    48: ThermalTuple('osap00', [92238], 1),
    49: ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3),
    50: ThermalTuple('oice', [8016, 8017, 8018], 1),
+    51: ThermalTuple('od2o', [8016, 8017, 8018], 1),
    52: ThermalTuple('mg24', [12024], 1),
    53: ThermalTuple('al27', [13027], 1),
    55: ThermalTuple('yyh2', [39089], 1),
    56: ThermalTuple('fe56', [26056], 1),
    58: ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1),
-    59: ThermalTuple('cacah2', [20040, 20042, 20043, 20044, 20046, 20048], 1),
-    75: ThermalTuple('ouo2', [8016, 8017, 8018], 1),
+    59: ThermalTuple('si00', [14028], 1),
+    60: ThermalTuple('asap00', [13027], 1),
+    71: ThermalTuple('n-un', [7014, 7015], 1),
+    72: ThermalTuple('u-un', [92238], 1),
+    75: ThermalTuple('uuo2', [8016, 8017, 8018], 1),
 }

+
+def _get_thermal_data(ev, mat):
+    """Return appropriate ThermalTuple, accounting for bugs."""
+
+    # JEFF assigns MAT=59 to Ca in CaH2 (which is supposed to be silicon).
+    if ev.info['library'][0] == 'JEFF':
+        if ev.material == 59:
+            if 'CaH2' in ''.join(ev.info['description']):
+                zaids = [20040, 20042, 20043, 20044, 20046, 20048]
+                return ThermalTuple('cacah2', zaids, 1)
+
+    # Before ENDF/B-VIII.0, crystalline graphite was MAT=31
+    if ev.info['library'] != ('ENDF/B', 8, 0):
+        if ev.material == 31:
+            return _THERMAL_DATA[30]
+
+    # ENDF/B incorrectly assigns MAT numbers for UO2
+    #
+    # Material | ENDF Manual | VII.0 | VII.1 | VIII.0
+    # ---------|-------------|-------|-------|-------
+    # O in UO2 |     45      |  75   |  75   |   75
+    # U in UO2 |     75      |  76   |  48   |   48
+    if ev.info['library'][0] == 'ENDF/B':
+        if ev.material == 75:
+            return _THERMAL_DATA[45]
+        version = ev.info['library'][1:]
+        if version in ((7, 1), (8, 0)) and ev.material == 48:
+            return _THERMAL_DATA[75]
+        if version == (7, 0) and ev.material == 76:
+            return _THERMAL_DATA[75]
+
+    # If not a problematic material, use the dictionary as is
+    return _THERMAL_DATA[mat]
+
+
 _TEMPLATE_RECONR = """
 reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%%
 {nendf} {npendf}
@ -444,7 +490,8 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
    mat_thermal = ev_thermal.material
    zsymam_thermal = ev_thermal.target['zsymam']

-    data = _THERMAL_DATA[mat_thermal]
+    # Determine name, isotopes based on MAT number
+    data = _get_thermal_data(ev_thermal, mat_thermal)
    zaids = ' '.join(str(zaid) for zaid in data.zaids[:3])

    # Determine name of library
--- a/openmc/data/photon.py
+++ b/openmc/data/photon.py
@ -549,7 +549,9 @@ class IncidentPhoton(EqualityMixin):
            ace = get_table(ace_or_filename)

        # Get atomic number based on name of ACE table
-        zaid = ace.name.split('.')[0]
+        zaid, xs = ace.name.split('.')
+        if not xs.endswith('p'):
+            raise TypeError("{} is not a photoatomic transport ACE table.".format(ace))
        Z = get_metadata(int(zaid))[2]

        # Read each reaction
--- a/openmc/data/resonance.py
+++ b/openmc/data/resonance.py
@ -430,7 +430,7 @@ class MultiLevelBreitWigner(ResonanceRange):

            # Determine penetration at modified energy for competitive reaction
            if gx > 0:
-                Ex = E + self.q[l]*(A + 1)/A
+                Ex = E + self.q_value[l]*(A + 1)/A
                rho = k*self.channel_radius[l](Ex)
                rhohat = k*self.scattering_radius[l](Ex)
                px[i], sx[i] = penetration_shift(l, rho)
@ -915,7 +915,9 @@ class Unresolved(ResonanceRange):
        Minimum energy of the unresolved resonance range in eV
    energy_max : float
        Maximum energy of the unresolved resonance range in eV
-    scatter : openmc.data.Function1D
+    channel : openmc.data.Function1D
+        Channel radii as a function of energy
+    scattering : openmc.data.Function1D
        Scattering radii as a function of energy

    Attributes
@ -923,6 +925,10 @@ class Unresolved(ResonanceRange):
    add_to_background : bool
        If True, file 3 contains partial cross sections to be added to the
        average unresolved cross sections calculated from parameters.
+    atomic_weight_ratio : float
+        Atomic weight ratio of the target nuclide
+    channel_radius : openmc.data.Function1D
+        Channel radii as a function of energy
    energies : Iterable of float
        Energies at which parameters are tabulated
    energy_max : float
@ -938,11 +944,13 @@ class Unresolved(ResonanceRange):

    """

-    def __init__(self, target_spin, energy_min, energy_max, scatter):
-        super().__init__(target_spin, energy_min, energy_max, None, scatter)
+    def __init__(self, target_spin, energy_min, energy_max, channel, scattering):
+        super().__init__(target_spin, energy_min, energy_max, channel,
+                         scattering)
        self.energies = None
        self.parameters = None
        self.add_to_background = False
+        self.atomic_weight_ratio = None

    @classmethod
    def from_endf(cls, file_obj, items, fission_widths):
@ -967,9 +975,9 @@ class Unresolved(ResonanceRange):
        """
        # Read energy-dependent scattering radius if present
        energy_min, energy_max = items[0:2]
-        nro = items[4]
+        nro, naps = items[4:6]
        if nro != 0:
-            params, scattering_radius = get_tab1_record(file_obj)
+            params, ape = get_tab1_record(file_obj)

        # Get SPI, AP, and LSSF
        formalism = items[3]
@ -977,22 +985,23 @@ class Unresolved(ResonanceRange):
            items = get_cont_record(file_obj)
            target_spin = items[0]
            if nro == 0:
-                scattering_radius = items[1]
+                ap = Polynomial((items[1],))
            add_to_background = (items[2] == 0)

        if not fission_widths and formalism == 1:
            # Case A -- fission widths not given, all parameters are
            # energy-independent
            NLS = items[4]
-            columns = ['L', 'J', 'd', 'amun', 'gn', 'gg', 'gf']
+            columns = ['L', 'J', 'd', 'amun', 'gn0', 'gg']
            records = []
            for ls in range(NLS):
                items, values = get_list_record(file_obj)
+                awri = items[0]
                l = items[2]
                NJS = items[5]
                for j in range(NJS):
-                    d, j, amun, gn, gg = values[6*j:6*j + 5]
-                    records.append([l, j, d, amun, gn, gg, 0.0])
+                    d, j, amun, gn0, gg = values[6*j:6*j + 5]
+                    records.append([l, j, d, amun, gn0, gg])
            parameters = pd.DataFrame.from_records(records, columns=columns)
            energies = None

@ -1002,14 +1011,14 @@ class Unresolved(ResonanceRange):
            items, energies = get_list_record(file_obj)
            target_spin = items[0]
            if nro == 0:
-                scattering_radius = items[1]
+                ap = Polynomial((items[1],))
            add_to_background = (items[2] == 0)
            NE, NLS = items[4:6]
-            l_values = np.zeros(NLS, int)
-            records = [[] for e in energies]
-            columns = ['L', 'J', 'd', 'amun', 'gn0', 'gg', 'gf']
+            records = []
+            columns = ['L', 'J', 'E', 'd', 'amun', 'amuf', 'gn0', 'gg', 'gf']
            for ls in range(NLS):
                items = get_cont_record(file_obj)
+                awri = items[0]
                l = items[2]
                NJS = items[4]
                for j in range(NJS):
@ -1020,18 +1029,20 @@ class Unresolved(ResonanceRange):
                    amun = values[2]
                    gn0 = values[3]
                    gg = values[4]
-                    gf = values[6:]
-                    for k, gf_i in enumerate(gf):
-                        records[k].append([l, j, d, amun, gn0, gg, gf])
-            parameters = [pd.DataFrame(r, columns=columns) for r in records]
+                    gfs = values[6:]
+                    for E, gf in zip(energies, gfs):
+                        records.append([l, j, E, d, amun, muf, gn0, gg, gf])
+            parameters = pd.DataFrame.from_records(records, columns=columns)

        elif formalism == 2:
            # Case C -- all parameters are energy-dependent
            NLS = items[4]
-            columns = ['L', 'J', 'E', 'd', 'gx', 'gn0', 'gg', 'gf']
+            columns = ['L', 'J', 'E', 'd', 'amux', 'amun', 'amuf', 'gx', 'gn0',
+                       'gg', 'gf']
            records = []
            for ls in range(NLS):
                items = get_cont_record(file_obj)
+                awri = items[0]
                l = items[2]
                NJS = items[4]
                for j in range(NJS):
@ -1040,8 +1051,8 @@ class Unresolved(ResonanceRange):
                    j = items[0]
                    amux = values[2]
                    amun = values[3]
-                    amug = values[4]
                    amuf = values[5]
+                    energies = []
                    for k in range(1, ne + 1):
                        E = values[6*k]
                        d = values[6*k + 1]
@ -1049,13 +1060,35 @@ class Unresolved(ResonanceRange):
                        gn0 = values[6*k + 3]
                        gg = values[6*k + 4]
                        gf = values[6*k + 5]
-                        records.append([l, j, E, d, gx, gn0, gg, gf])
+                        energies.append(E)
+                        records.append([l, j, E, d, amux, amun, amuf, gx, gn0,
+                                        gg, gf])
            parameters = pd.DataFrame.from_records(records, columns=columns)
-            energies = None

-        urr = cls(target_spin, energy_min, energy_max, scattering_radius)
+        # Calculate channel radius from ENDF-102 equation D.14
+        a = Polynomial((0.123 * (NEUTRON_MASS*awri)**(1./3.) + 0.08,))
+
+        # Determine scattering and channel radius
+        if nro == 0:
+            scattering_radius = ap
+            if naps == 0:
+                channel_radius = a
+            elif naps == 1:
+                channel_radius = ap
+        elif nro == 1:
+            scattering_radius = ape
+            if naps == 0:
+                channel_radius = a
+            elif naps == 1:
+                channel_radius = ape
+            elif naps == 2:
+                channel_radius = ap
+
+        urr = cls(target_spin, energy_min, energy_max, channel_radius,
+                  scattering_radius)
        urr.parameters = parameters
        urr.add_to_background = add_to_background
+        urr.atomic_weight_ratio = awri
        urr.energies = energies

        return urr
--- a/openmc/data/thermal.py
+++ b/openmc/data/thermal.py
@ -33,10 +33,12 @@ _THERMAL_NAMES = {
    'c_Be': ('be', 'be-metal', 'be-met', 'be00'),
    'c_BeO': ('beo',),
    'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o', 'bbeo00'),
+    'c_Be_in_Be2C': ('bebe2c',),
    'c_C6H6': ('benz', 'c6h6'),
    'c_C_in_SiC': ('csic', 'c-sic'),
    'c_Ca_in_CaH2': ('cah', 'cah00'),
    'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw', 'dhw00'),
+    'c_D_in_D2O_ice': ('dice',),
    'c_Fe56': ('fe', 'fe56', 'fe-56'),
    'c_Graphite': ('graph', 'grph', 'gr', 'gr00'),
    'c_Graphite_10p': ('grph10',),
@ -45,6 +47,7 @@ _THERMAL_NAMES = {
    'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00'),
    'c_H_in_CH4_liquid': ('lch4', 'lmeth'),
    'c_H_in_CH4_solid': ('sch4', 'smeth'),
+    'c_H_in_CH4_solid_phase_II': ('sch4p2',),
    'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw', 'lw00'),
    'c_H_in_H2O_solid': ('hice', 'h-ice', 'ice00'),
    'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci'),
@ -600,6 +603,8 @@ class ThermalScattering(EqualityMixin):

        # Get new name that is GND-consistent
        ace_name, xs = ace.name.split('.')
+        if not xs.endswith('t'):
+            raise TypeError("{} is not a thermal scattering ACE table.".format(ace))
        name = get_thermal_name(ace_name)

        # Assign temperature to the running list
--- a/openmc/deplete/init.py
+++ b/openmc/deplete/init.py
@ -4,12 +4,11 @@ openmc.deplete

 A depletion front-end tool.
 """
-from sys import exit
+import sys
+from unittest.mock import Mock

 from h5py import get_config

-from unittest.mock import Mock
-
 from .dummy_comm import DummyCommunicator

 try:
@ -22,7 +21,7 @@ try:
    if not get_config().mpi and comm.size > 1:
        # Raise exception only on process 0
        if comm.rank:
-            exit()
+            sys.exit()
        raise RuntimeError(
            "Need parallel HDF5 installed to perform depletion with MPI"
        )
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@ -28,7 +28,7 @@ from .results_list import ResultsList
 __all__ = [
    "OperatorResult", "TransportOperator", "ReactionRateHelper",
    "EnergyHelper", "FissionYieldHelper", "TalliedFissionYieldHelper",
-    "Integrator", "SIIntegrator"]
+    "Integrator", "SIIntegrator", "DepSystemSolver"]


 OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
@ -546,6 +546,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper):

        # Tally group-wise fission reaction rates
        self._fission_rate_tally = Tally()
+        self._fission_rate_tally.writable = False
        self._fission_rate_tally.scores = ['fission']

        self._fission_rate_tally.filters = [MaterialFilter(materials)]
@ -814,11 +815,11 @@ class SIIntegrator(Integrator):
        reset_particles = False
        if step_index == 0 and hasattr(self.operator, "settings"):
            reset_particles = True
-            self.operator.settings.particles *= self.n_stages
+            self.operator.settings.particles *= self.n_steps
        inherited = super()._get_bos_data_from_operator(
            step_index, step_power, bos_conc)
        if reset_particles:
-            self.operator.settings.particles //= self.n_stages
+            self.operator.settings.particles //= self.n_steps
        return inherited

    def integrate(self):
@ -855,3 +856,40 @@ class SIIntegrator(Integrator):
            Results.save(self.operator, [conc], [res_list[-1]], [t, t],
                         p, self._i_res + len(self), proc_time)
            self.operator.write_bos_data(self._i_res + len(self))
+
+
+class DepSystemSolver(ABC):
+    r"""Abstract class for solving depletion equations
+
+    Responsible for solving
+
+    .. math::
+
+        \frac{\partial \vec{N}}{\partial t} = \bar{A}\vec{N}(t),
+
+    for :math:`0< t\leq t +\Delta t`, given :math:`\vec{N}(0) = \vec{N}_0`
+
+    """
+
+    @abstractmethod
+    def __call__(self, A, n0, dt):
+        """Solve the linear system of equations for depletion
+
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved
+
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``. Should be of identical shape
+            to ``n0``.
+
+        """
--- a/openmc/deplete/cram.py
+++ b/openmc/deplete/cram.py
@ -3,6 +3,7 @@
 Implements two different forms of CRAM for use in openmc.deplete.
 """

+import numbers
 from itertools import repeat
 from multiprocessing import Pool
 import time
@ -11,7 +12,12 @@ import numpy as np
 import scipy.sparse as sp
 import scipy.sparse.linalg as sla

-from . import comm
+from openmc.checkvalue import check_type, check_length
+from .abc import DepSystemSolver
+
+__all__ = [
+    "deplete", "timed_deplete", "CRAM16", "CRAM48",
+    "Cram16Solver", "Cram48Solver", "IPFCramSolver"]


 def deplete(chain, x, rates, dt, matrix_func=None):
@ -81,147 +87,169 @@ def timed_deplete(*args, **kwargs):
    return time.time() - start, results


-def CRAM16(A, n0, dt):
-    """Chebyshev Rational Approximation Method, order 16
+class IPFCramSolver(DepSystemSolver):
+    r"""CRAM depletion solver that uses incomplete partial factorization

-    Algorithm is the 16th order Chebyshev Rational Approximation Method,
-    implemented in the more stable `incomplete partial fraction (IPF)
-    <https://doi.org/10.13182/NSE15-26>`_ form.
+    Provides a :meth:`__call__` that utilizes an incomplete
+    partial factorization (IPF) for the Chebyshev Rational Approximation
+    Method (CRAM), as described in the following paper: M. Pusa, "`Higher-Order
+    Chebyshev Rational Approximation Method and Application to Burnup Equations
+    <https://doi.org/10.13182/NSE15-26>`_," Nucl. Sci. Eng., 182:3, 297-318.

    Parameters
    ----------
-    A : scipy.linalg.csr_matrix
-        Matrix to take exponent of.
-    n0 : numpy.array
-        Vector to operate a matrix exponent on.
-    dt : float
-        Time to integrate to.
+    alpha : numpy.ndarray
+        Complex residues of poles used in the factorization. Must be a
+        vector with even number of items.
+    theta : numpy.ndarray
+        Complex poles. Must have an equal size as ``alpha``.
+    alpha0 : float
+        Limit of the approximation at infinity

-    Returns
-    -------
-    numpy.array
-        Results of the matrix exponent.
-
-    """
-
-    alpha = np.array([+2.124853710495224e-16,
-                      +5.464930576870210e+3 - 3.797983575308356e+4j,
-                      +9.045112476907548e+1 - 1.115537522430261e+3j,
-                      +2.344818070467641e+2 - 4.228020157070496e+2j,
-                      +9.453304067358312e+1 - 2.951294291446048e+2j,
-                      +7.283792954673409e+2 - 1.205646080220011e+5j,
-                      +3.648229059594851e+1 - 1.155509621409682e+2j,
-                      +2.547321630156819e+1 - 2.639500283021502e+1j,
-                      +2.394538338734709e+1 - 5.650522971778156e+0j],
-                     dtype=np.complex128)
-    theta = np.array([+0.0,
-                      +3.509103608414918 + 8.436198985884374j,
-                      +5.948152268951177 + 3.587457362018322j,
-                      -5.264971343442647 + 16.22022147316793j,
-                      +1.419375897185666 + 10.92536348449672j,
-                      +6.416177699099435 + 1.194122393370139j,
-                      +4.993174737717997 + 5.996881713603942j,
-                      -1.413928462488886 + 13.49772569889275j,
-                      -10.84391707869699 + 19.27744616718165j],
-                     dtype=np.complex128)
-
-    n = A.shape[0]
-
-    alpha0 = 2.124853710495224e-16
-
-    k = 8
-
-    y = np.array(n0, dtype=np.float64)
-    for l in range(1, k+1):
-        y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
-
-    y *= alpha0
-    return y
-
-
-def CRAM48(A, n0, dt):
-    """Chebyshev Rational Approximation Method, order 48
-
-    Algorithm is the 48th order Chebyshev Rational Approximation Method,
-    implemented in the more stable `incomplete partial fraction (IPF)
-    <https://doi.org/10.13182/NSE15-26>`_ form.
-
-    Parameters
+    Attributes
    ----------
-    A : scipy.linalg.csr_matrix
-        Matrix to take exponent of.
-    n0 : numpy.array
-        Vector to operate a matrix exponent on.
-    dt : float
-        Time to integrate to.
-
-    Returns
-    -------
-    numpy.array
-        Results of the matrix exponent.
+    alpha : numpy.ndarray
+        Complex residues of poles :attr:`theta` in the incomplete partial
+        factorization. Denoted as :math:`\tilde{\alpha}`
+    theta : numpy.ndarray
+        Complex poles :math:`\theta` of the rational approximation
+    alpha0 : float
+        Limit of the approximation at infinity

    """

-    theta_r = np.array([-4.465731934165702e+1, -5.284616241568964e+0,
-                        -8.867715667624458e+0, +3.493013124279215e+0,
-                        +1.564102508858634e+1, +1.742097597385893e+1,
-                        -2.834466755180654e+1, +1.661569367939544e+1,
-                        +8.011836167974721e+0, -2.056267541998229e+0,
-                        +1.449208170441839e+1, +1.853807176907916e+1,
-                        +9.932562704505182e+0, -2.244223871767187e+1,
-                        +8.590014121680897e-1, -1.286192925744479e+1,
-                        +1.164596909542055e+1, +1.806076684783089e+1,
-                        +5.870672154659249e+0, -3.542938819659747e+1,
-                        +1.901323489060250e+1, +1.885508331552577e+1,
-                        -1.734689708174982e+1, +1.316284237125190e+1])
-    theta_i = np.array([+6.233225190695437e+1, +4.057499381311059e+1,
-                        +4.325515754166724e+1, +3.281615453173585e+1,
-                        +1.558061616372237e+1, +1.076629305714420e+1,
-                        +5.492841024648724e+1, +1.316994930024688e+1,
-                        +2.780232111309410e+1, +3.794824788914354e+1,
-                        +1.799988210051809e+1, +5.974332563100539e+0,
-                        +2.532823409972962e+1, +5.179633600312162e+1,
-                        +3.536456194294350e+1, +4.600304902833652e+1,
-                        +2.287153304140217e+1, +8.368200580099821e+0,
-                        +3.029700159040121e+1, +5.834381701800013e+1,
-                        +1.194282058271408e+0, +3.583428564427879e+0,
-                        +4.883941101108207e+1, +2.042951874827759e+1])
-    theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
+    def __init__(self, alpha, theta, alpha0):
+        check_type("alpha", alpha, np.ndarray, numbers.Complex)
+        check_type("theta", theta, np.ndarray, numbers.Complex)
+        check_length("theta", theta, alpha.size)
+        check_type("alpha0", alpha0, numbers.Real)
+        self.alpha = alpha
+        self.theta = theta
+        self.alpha0 = alpha0

-    alpha_r = np.array([+6.387380733878774e+2, +1.909896179065730e+2,
-                        +4.236195226571914e+2, +4.645770595258726e+2,
-                        +7.765163276752433e+2, +1.907115136768522e+3,
-                        +2.909892685603256e+3, +1.944772206620450e+2,
-                        +1.382799786972332e+5, +5.628442079602433e+3,
-                        +2.151681283794220e+2, +1.324720240514420e+3,
-                        +1.617548476343347e+4, +1.112729040439685e+2,
-                        +1.074624783191125e+2, +8.835727765158191e+1,
-                        +9.354078136054179e+1, +9.418142823531573e+1,
-                        +1.040012390717851e+2, +6.861882624343235e+1,
-                        +8.766654491283722e+1, +1.056007619389650e+2,
-                        +7.738987569039419e+1, +1.041366366475571e+2])
-    alpha_i = np.array([-6.743912502859256e+2, -3.973203432721332e+2,
-                        -2.041233768918671e+3, -1.652917287299683e+3,
-                        -1.783617639907328e+4, -5.887068595142284e+4,
-                        -9.953255345514560e+3, -1.427131226068449e+3,
-                        -3.256885197214938e+6, -2.924284515884309e+4,
-                        -1.121774011188224e+3, -6.370088443140973e+4,
-                        -1.008798413156542e+6, -8.837109731680418e+1,
-                        -1.457246116408180e+2, -6.388286188419360e+1,
-                        -2.195424319460237e+2, -6.719055740098035e+2,
-                        -1.693747595553868e+2, -1.177598523430493e+1,
-                        -4.596464999363902e+3, -1.738294585524067e+3,
-                        -4.311715386228984e+1, -2.777743732451969e+2])
-    alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
-    n = A.shape[0]
+    def __call__(self, A, n0, dt):
+        """Solve depletion equations using IPF CRAM

-    alpha0 = 2.258038182743983e-47
+        Parameters
+        ----------
+        A : scipy.sparse.csr_matrix
+            Sparse transmutation matrix ``A[j, i]`` desribing rates at
+            which isotope ``i`` transmutes to isotope ``j``
+        n0 : numpy.ndarray
+            Initial compositions, typically given in number of atoms in some
+            material or an atom density
+        dt : float
+            Time [s] of the specific interval to be solved

-    k = 24
+        Returns
+        -------
+        numpy.ndarray
+            Final compositions after ``dt``

-    y = np.array(n0, dtype=np.float64)
-    for l in range(k):
-        y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
+        """
+        A = sp.csr_matrix(A * dt, dtype=np.float64)
+        y = np.asarray(n0, dtype=np.float64)
+        ident = sp.eye(A.shape[0])
+        for alpha, theta in zip(self.alpha, self.theta):
+            y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
+        return y * self.alpha0

-    y *= alpha0
-    return y
+
+# Coefficients for IPF Cram 16
+c16_alpha = np.array([
+    +5.464930576870210e+3 - 3.797983575308356e+4j,
+    +9.045112476907548e+1 - 1.115537522430261e+3j,
+    +2.344818070467641e+2 - 4.228020157070496e+2j,
+    +9.453304067358312e+1 - 2.951294291446048e+2j,
+    +7.283792954673409e+2 - 1.205646080220011e+5j,
+    +3.648229059594851e+1 - 1.155509621409682e+2j,
+    +2.547321630156819e+1 - 2.639500283021502e+1j,
+    +2.394538338734709e+1 - 5.650522971778156e+0j],
+    dtype=np.complex128)
+
+c16_theta = np.array([
+    +3.509103608414918 + 8.436198985884374j,
+    +5.948152268951177 + 3.587457362018322j,
+    -5.264971343442647 + 16.22022147316793j,
+    +1.419375897185666 + 10.92536348449672j,
+    +6.416177699099435 + 1.194122393370139j,
+    +4.993174737717997 + 5.996881713603942j,
+    -1.413928462488886 + 13.49772569889275j,
+    -10.84391707869699 + 19.27744616718165j],
+    dtype=np.complex128)
+
+c16_alpha0 = 2.124853710495224e-16
+Cram16Solver = IPFCramSolver(c16_alpha, c16_theta, c16_alpha0)
+CRAM16 = Cram16Solver.__call__
+
+del c16_alpha, c16_alpha0, c16_theta
+
+# Coefficients for 48th order IPF Cram
+
+theta_r = np.array([
+    -4.465731934165702e+1, -5.284616241568964e+0,
+    -8.867715667624458e+0, +3.493013124279215e+0,
+    +1.564102508858634e+1, +1.742097597385893e+1,
+    -2.834466755180654e+1, +1.661569367939544e+1,
+    +8.011836167974721e+0, -2.056267541998229e+0,
+    +1.449208170441839e+1, +1.853807176907916e+1,
+    +9.932562704505182e+0, -2.244223871767187e+1,
+    +8.590014121680897e-1, -1.286192925744479e+1,
+    +1.164596909542055e+1, +1.806076684783089e+1,
+    +5.870672154659249e+0, -3.542938819659747e+1,
+    +1.901323489060250e+1, +1.885508331552577e+1,
+    -1.734689708174982e+1, +1.316284237125190e+1])
+
+theta_i = np.array([
+    +6.233225190695437e+1, +4.057499381311059e+1,
+    +4.325515754166724e+1, +3.281615453173585e+1,
+    +1.558061616372237e+1, +1.076629305714420e+1,
+    +5.492841024648724e+1, +1.316994930024688e+1,
+    +2.780232111309410e+1, +3.794824788914354e+1,
+    +1.799988210051809e+1, +5.974332563100539e+0,
+    +2.532823409972962e+1, +5.179633600312162e+1,
+    +3.536456194294350e+1, +4.600304902833652e+1,
+    +2.287153304140217e+1, +8.368200580099821e+0,
+    +3.029700159040121e+1, +5.834381701800013e+1,
+    +1.194282058271408e+0, +3.583428564427879e+0,
+    +4.883941101108207e+1, +2.042951874827759e+1])
+
+c48_theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
+
+alpha_r = np.array([
+    +6.387380733878774e+2, +1.909896179065730e+2,
+    +4.236195226571914e+2, +4.645770595258726e+2,
+    +7.765163276752433e+2, +1.907115136768522e+3,
+    +2.909892685603256e+3, +1.944772206620450e+2,
+    +1.382799786972332e+5, +5.628442079602433e+3,
+    +2.151681283794220e+2, +1.324720240514420e+3,
+    +1.617548476343347e+4, +1.112729040439685e+2,
+    +1.074624783191125e+2, +8.835727765158191e+1,
+    +9.354078136054179e+1, +9.418142823531573e+1,
+    +1.040012390717851e+2, +6.861882624343235e+1,
+    +8.766654491283722e+1, +1.056007619389650e+2,
+    +7.738987569039419e+1, +1.041366366475571e+2])
+
+alpha_i = np.array([
+    -6.743912502859256e+2, -3.973203432721332e+2,
+    -2.041233768918671e+3, -1.652917287299683e+3,
+    -1.783617639907328e+4, -5.887068595142284e+4,
+    -9.953255345514560e+3, -1.427131226068449e+3,
+    -3.256885197214938e+6, -2.924284515884309e+4,
+    -1.121774011188224e+3, -6.370088443140973e+4,
+    -1.008798413156542e+6, -8.837109731680418e+1,
+    -1.457246116408180e+2, -6.388286188419360e+1,
+    -2.195424319460237e+2, -6.719055740098035e+2,
+    -1.693747595553868e+2, -1.177598523430493e+1,
+    -4.596464999363902e+3, -1.738294585524067e+3,
+    -4.311715386228984e+1, -2.777743732451969e+2])
+
+c48_alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
+
+c48_alpha0 = 2.258038182743983e-47
+
+Cram48Solver = IPFCramSolver(c48_alpha, c48_theta, c48_alpha0)
+
+del c48_alpha, c48_alpha0, c48_theta, alpha_r, alpha_i, theta_r, theta_i
+
+CRAM48 = Cram48Solver.__call__
--- a/openmc/deplete/helpers.py
+++ b/openmc/deplete/helpers.py
@ -59,6 +59,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
            ``"(n, gamma)"``, needed for the reaction rate tally.
        """
        self._rate_tally = Tally()
+        self._rate_tally.writable = False
        self._rate_tally.scores = scores
        self._rate_tally.filters = [MaterialFilter(materials)]

@ -194,6 +195,7 @@ class EnergyScoreHelper(EnergyHelper):

        """
        self._tally = Tally()
+        self._tally.writable = False
        self._tally.scores = [self.score]

    def reset(self):
@ -570,6 +572,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
        func_filter = EnergyFunctionFilter()
        func_filter.set_data((0, self._upper_energy), (0, self._upper_energy))
        weighted_tally = Tally()
+        weighted_tally.writable = False
        weighted_tally.scores = ['fission']
        weighted_tally.filters = filters + [func_filter]
        self._weighted_tally = weighted_tally
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@ -82,6 +82,7 @@ class Operator(TransportOperator):
        in the previous results.
    diff_burnable_mats : bool, optional
        Whether to differentiate burnable materials with multiple instances.
+        Volumes are divided equally from the original material volume.
        Default: False.
    energy_mode : {"energy-deposition", "fission-q"}
        Indicator for computing system energy. ``"energy-deposition"`` will
@ -309,7 +310,12 @@ class Operator(TransportOperator):
            # Assign distribmats to cells
            for cell in self.geometry.get_all_material_cells().values():
                if cell.fill in distribmats and cell.num_instances > 1:
-                    cell.fill = [cell.fill.clone()
+                    mat = cell.fill
+                    if mat.volume is None:
+                        raise RuntimeError("Volume not specified for depletable "
+                                           "material with ID={}.".format(mat.id))
+                    mat.volume /= mat.num_instances
+                    cell.fill = [mat.clone()
                                 for i in range(cell.num_instances)]

    def _get_burnable_mats(self):
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@ -88,7 +88,7 @@ class Geometry(object):
        """
        # Create XML representation
        root_element = ET.Element("geometry")
-        self.root_universe.create_xml_subelement(root_element)
+        self.root_universe.create_xml_subelement(root_element, memo=set())

        # Sort the elements in the file
        root_element[:] = sorted(root_element, key=lambda x: (
@ -272,9 +272,9 @@ class Geometry(object):

        """
        if self.root_universe is not None:
-            return self.root_universe.get_all_cells()
+            return self.root_universe.get_all_cells(memo=set())
        else:
-            return []
+            return OrderedDict()

    def get_all_universes(self):
        """Return all universes in the geometry.
@ -301,7 +301,10 @@ class Geometry(object):
            instances

        """
-        return self.root_universe.get_all_materials()
+        if self.root_universe is not None:
+            return self.root_universe.get_all_materials(memo=set())
+        else:
+            return OrderedDict()

    def get_all_material_cells(self):
        """Return all cells filled by a material
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@ -332,7 +332,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):

        return nuclides

-    def get_all_cells(self):
+    def get_all_cells(self, memo=None):
        """Return all cells that are contained within the lattice

        Returns
@ -342,16 +342,22 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
            instances

        """
-
        cells = OrderedDict()
+
+        if memo and self in memo:
+            return cells
+
+        if memo is not None:
+            memo.add(self)
+
        unique_universes = self.get_unique_universes()

        for universe_id, universe in unique_universes.items():
-            cells.update(universe.get_all_cells())
+            cells.update(universe.get_all_cells(memo))

        return cells

-    def get_all_materials(self):
+    def get_all_materials(self, memo=None):
        """Return all materials that are contained within the lattice

        Returns
@ -365,9 +371,9 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
        materials = OrderedDict()

        # Append all Cells in each Cell in the Universe to the dictionary
-        cells = self.get_all_cells()
+        cells = self.get_all_cells(memo)
        for cell_id, cell in cells.items():
-            materials.update(cell.get_all_materials())
+            materials.update(cell.get_all_materials(memo))

        return materials

@ -468,7 +474,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
        if memo is None:
            memo = {}

-        # If no nemoize'd clone exists, instantiate one
+        # If no memoize'd clone exists, instantiate one
        if self not in memo:
            clone = deepcopy(self)
            clone.id = None
@ -753,15 +759,29 @@ class RectLattice(Lattice):
                    0 <= idx[1] < self.shape[1] and
                    0 <= idx[2] < self.shape[2])

-    def create_xml_subelement(self, xml_element):
+    def create_xml_subelement(self, xml_element, memo=None):
+        """Add the lattice xml representation to an incoming xml element

-        # Determine if XML element already contains subelement for this Lattice
-        path = './lattice[@id=\'{0}\']'.format(self._id)
-        test = xml_element.find(path)
+        Parameters
+        ----------
+        xml_element : xml.etree.ElementTree.Element
+            XML element to be added to

-        # If the element does contain the Lattice subelement, then return
-        if test is not None:
+        memo : set or None
+            A set of object id's representing geometry entities already
+            written to the xml_element. This parameter is used internally
+            and should not be specified by users.
+
+        Returns
+        -------
+        None
+
+        """
+        # If the element already contains the Lattice subelement, then return
+        if memo and self in memo:
            return
+        if memo is not None:
+            memo.add(self)

        lattice_subelement = ET.Element("lattice")
        lattice_subelement.set("id", str(self._id))
@ -777,7 +797,7 @@ class RectLattice(Lattice):
        if self._outer is not None:
            outer = ET.SubElement(lattice_subelement, "outer")
            outer.text = '{0}'.format(self._outer._id)
-            self._outer.create_xml_subelement(xml_element)
+            self._outer.create_xml_subelement(xml_element, memo)

        # Export Lattice cell dimensions
        dimension = ET.SubElement(lattice_subelement, "dimension")
@ -801,7 +821,7 @@ class RectLattice(Lattice):
                        universe_ids += '{0} '.format(universe._id)

                        # Create XML subelement for this Universe
-                        universe.create_xml_subelement(xml_element)
+                        universe.create_xml_subelement(xml_element, memo)

                    # Add newline character when we reach end of row of cells
                    universe_ids += '\n'
@ -819,7 +839,7 @@ class RectLattice(Lattice):
                    universe_ids += '{0} '.format(universe._id)

                    # Create XML subelement for this Universe
-                    universe.create_xml_subelement(xml_element)
+                    universe.create_xml_subelement(xml_element, memo)

                # Add newline character when we reach end of row of cells
                universe_ids += '\n'
@ -1274,14 +1294,12 @@ class HexLattice(Lattice):
        else:
            return g < self.num_rings and 0 <= idx[2] < self.num_axial

-    def create_xml_subelement(self, xml_element):
-        # Determine if XML element already contains subelement for this Lattice
-        path = './hex_lattice[@id=\'{0}\']'.format(self._id)
-        test = xml_element.find(path)
-
-        # If the element does contain the Lattice subelement, then return
-        if test is not None:
+    def create_xml_subelement(self, xml_element, memo=None):
+        # If this subelement has already been written, return
+        if memo and self in memo:
            return
+        if memo is not None:
+            memo.add(self)

        lattice_subelement = ET.Element("hex_lattice")
        lattice_subelement.set("id", str(self._id))
@ -1297,7 +1315,7 @@ class HexLattice(Lattice):
        if self._outer is not None:
            outer = ET.SubElement(lattice_subelement, "outer")
            outer.text = '{0}'.format(self._outer._id)
-            self._outer.create_xml_subelement(xml_element)
+            self._outer.create_xml_subelement(xml_element, memo)

        lattice_subelement.set("n_rings", str(self._num_rings))
        # If orientation is "x" export it to XML
@ -1319,13 +1337,13 @@ class HexLattice(Lattice):
            for z in range(self._num_axial):
                # Initialize the center universe.
                universe = self._universes[z][-1][0]
-                universe.create_xml_subelement(xml_element)
+                universe.create_xml_subelement(xml_element, memo)

                # Initialize the remaining universes.
                for r in range(self._num_rings-1):
                    for theta in range(6*(self._num_rings - 1 - r)):
                        universe = self._universes[z][r][theta]
-                        universe.create_xml_subelement(xml_element)
+                        universe.create_xml_subelement(xml_element, memo)

                # Get a string representation of the universe IDs.
                slices.append(self._repr_axial_slice(self._universes[z]))
@ -1337,13 +1355,13 @@ class HexLattice(Lattice):
        else:
            # Initialize the center universe.
            universe = self._universes[-1][0]
-            universe.create_xml_subelement(xml_element)
+            universe.create_xml_subelement(xml_element, memo)

            # Initialize the remaining universes.
            for r in range(self._num_rings - 1):
                for theta in range(6*(self._num_rings - 1 - r)):
                    universe = self._universes[r][theta]
-                    universe.create_xml_subelement(xml_element)
+                    universe.create_xml_subelement(xml_element, memo)

            # Get a string representation of the universe IDs.
            universe_ids = self._repr_axial_slice(self._universes)
--- a/openmc/lib/settings.py
+++ b/openmc/lib/settings.py
@ -26,6 +26,7 @@ class _Settings(object):
    restart_run = _DLLGlobal(c_bool, 'restart_run')
    run_CE = _DLLGlobal(c_bool, 'run_CE')
    verbosity = _DLLGlobal(c_int, 'verbosity')
+    output_summary = _DLLGlobal(c_bool, 'output_summary')

    @property
    def run_mode(self):
--- a/openmc/lib/tally.py
+++ b/openmc/lib/tally.py
@ -53,6 +53,9 @@ _dll.openmc_tally_get_scores.errcheck = _error_handler
 _dll.openmc_tally_get_type.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_tally_get_type.restype = c_int
 _dll.openmc_tally_get_type.errcheck = _error_handler
+_dll.openmc_tally_get_writable.argtypes = [c_int32, POINTER(c_bool)]
+_dll.openmc_tally_get_writable.restype = c_int
+_dll.openmc_tally_get_writable.errcheck = _error_handler
 _dll.openmc_tally_reset.argtypes = [c_int32]
 _dll.openmc_tally_reset.restype = c_int
 _dll.openmc_tally_reset.errcheck = _error_handler
@ -81,6 +84,9 @@ _dll.openmc_tally_set_scores.errcheck = _error_handler
 _dll.openmc_tally_set_type.argtypes = [c_int32, c_char_p]
 _dll.openmc_tally_set_type.restype = c_int
 _dll.openmc_tally_set_type.errcheck = _error_handler
+_dll.openmc_tally_set_writable.argtypes = [c_int32, c_bool]
+_dll.openmc_tally_set_writable.restype = c_int
+_dll.openmc_tally_set_writable.errcheck = _error_handler
 _dll.tallies_size.restype = c_size_t


@ -344,6 +350,16 @@ class Tally(_FortranObjectWithID):

        return std_dev

+    @property
+    def writable(self):
+        writable = c_bool()
+        _dll.openmc_tally_get_writable(self._index, writable)
+        return writable.value
+
+    @writable.setter
+    def writable(self, writable):
+        _dll.openmc_tally_set_writable(self._index, writable)
+
    def reset(self):
        """Reset results and num_realizations of tally"""
        _dll.openmc_tally_reset(self._index)
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@ -53,6 +53,12 @@ class MeshBase(IDManagerMixin, metaclass=ABCMeta):
        else:
            self._name = ''

+    def __repr__(self):
+        string = type(self).__name__ + '\n'
+        string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
+        string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
+        return string
+
    @classmethod
    def from_hdf5(cls, group):
        """Create mesh from HDF5 group
@ -126,7 +132,10 @@ class RegularMesh(MeshBase):

    @property
    def n_dimension(self):
-        return len(self._dimension)
+        if self._dimension is not None:
+            return len(self._dimension)
+        else:
+            return None

    @property
    def lower_left(self):
@ -187,11 +196,9 @@ class RegularMesh(MeshBase):
        self._width = width

    def __repr__(self):
-        string = 'RegularMesh\n'
-        string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
-        string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
-        string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
-        string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._dimension)
+        string = super().__repr__()
+        string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
+        string += '{0: <16}{1}{2}\n'.format('\tMesh Cells', '=\t', self._dimension)
        string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._lower_left)
        string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._upper_right)
        string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width)
@ -522,6 +529,27 @@ class RectilinearMesh(MeshBase):
        cv.check_type('mesh z_grid', grid, Iterable, Real)
        self._z_grid = grid

+    def __repr__(self):
+        fmt = '{0: <16}{1}{2}\n'
+        string = super().__repr__()
+        string += fmt.format('\tDimensions', '=\t', self.n_dimension)
+        x_grid_str = str(self._x_grid) if not self._x_grid else len(self._x_grid)
+        string += fmt.format('\tN X pnts:', '=\t', x_grid_str)
+        if self._x_grid:
+            string += fmt.format('\tX Min:', '=\t', self._x_grid[0])
+            string += fmt.format('\tX Max:', '=\t', self._x_grid[-1])
+        y_grid_str = str(self._y_grid) if not self._y_grid else len(self._y_grid)
+        string += fmt.format('\tN Y pnts:', '=\t', y_grid_str)
+        if self._y_grid:
+            string += fmt.format('\tY Min:', '=\t', self._y_grid[0])
+            string += fmt.format('\tY Max:', '=\t', self._y_grid[-1])
+        z_grid_str = str(self._z_grid) if not self._z_grid else len(self._z_grid)
+        string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
+        if self._z_grid:
+            string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
+            string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
+        return string
+
    @classmethod
    def from_hdf5(cls, group):
        mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
--- a/openmc/mgxs/init.py
+++ b/openmc/mgxs/init.py
@ -18,7 +18,7 @@ GROUP_STRUCTURES = {}

 .. _CASMO: https://www.studsvik.com/SharepointFiles/CASMO-5%20Development%20and%20Applications.pdf
 .. _XMAS-172: https://www-nds.iaea.org/wimsd/energy.htm
-.. _SHEM-361: https://www.polymtl.ca/merlin/libraries.htm
+.. _SHEM-361: https://www.polymtl.ca/merlin/downloads/FP214.pdf
 .. _activation: https://fispact.ukaea.uk/wiki/Keyword:GETXS
 .. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure
 .. _UKAEA-1102: https://fispact.ukaea.uk/wiki/UKAEA-1102_group_structure
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@ -986,7 +986,7 @@ class Library(object):
            xsdata.num_azimuthal = self.num_azimuthal

        if nuclide != 'total':
-            xsdata.atomic_weight_ratio = self._nuclides[nuclide][1]
+            xsdata.atomic_weight_ratio = self._nuclides[nuclide]

        if subdomain is None:
            subdomain = 'all'
--- a/openmc/model/triso.py
+++ b/openmc/model/triso.py
@ -827,13 +827,17 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
    triso_locations = {idx: [] for idx in indices}
    for t in trisos:
        for idx in t.classify(lattice):
-            if idx in sorted(triso_locations):
+            if idx in triso_locations:
                # Create copy of TRISO particle with materials preserved and
                # different cell/surface IDs
-                t_copy = copy.deepcopy(t)
+                t_copy = copy.copy(t)
                t_copy.id = None
                t_copy.fill = t.fill
-                t_copy._surface.id = None
+                t_copy._surface = openmc.Sphere(r=t._surface.r,
+                                                x0=t._surface.x0,
+                                                y0=t._surface.y0,
+                                                z0=t._surface.z0)
+                t_copy.region = -t_copy._surface
                triso_locations[idx].append(t_copy)
            else:
                warnings.warn('TRISO particle is partially or completely '
--- a/openmc/plots.py
+++ b/openmc/plots.py
@ -173,7 +173,7 @@ class Plot(IDManagerMixin):
    OpenMC is capable of generating two-dimensional slice plots and
    three-dimensional voxel plots. Colors that are used in plots can be given as
    RGB tuples, e.g. (255, 255, 255) would be white, or by a string indicating a
-    valid `SVG color <https://www.w3.org/TR/SVG/types.html#ColorKeywords>`_.
+    valid `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.

    Parameters
    ----------
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@ -375,13 +375,18 @@ class StatePoint(object):
                for tally_id in tally_ids:
                    group = tallies_group['tally {}'.format(tally_id)]

-                    # Read the number of realizations
-                    n_realizations = group['n_realizations'][()]
+                    # Check if tally is internal and therefore has no data
+                    if group.attrs.get("internal"):
+                        continue

                    # Create Tally object and assign basic properties
                    tally = openmc.Tally(tally_id)
                    tally._sp_filename = self._f.filename
                    tally.name = group['name'][()].decode() if 'name' in group else ''
+
+                    # Read the number of realizations
+                    n_realizations = group['n_realizations'][()]
+
                    tally.estimator = group['estimator'][()].decode()
                    tally.num_realizations = n_realizations

--- a/openmc/stats/multivariate.py
+++ b/openmc/stats/multivariate.py
@ -272,6 +272,8 @@ class Spatial(metaclass=ABCMeta):
        distribution = get_text(elem, 'type')
        if distribution == 'cartesian':
            return CartesianIndependent.from_xml_element(elem)
+        elif distribution == 'spherical':
+            return SphericalIndependent.from_xml_element(elem)
        elif distribution == 'box' or distribution == 'fission':
            return Box.from_xml_element(elem)
        elif distribution == 'point':
@ -281,7 +283,7 @@ class Spatial(metaclass=ABCMeta):
 class CartesianIndependent(Spatial):
    """Spatial distribution with independent x, y, and z distributions.

-    This distribution allows one to specify a coordinates whose x-, y-, and z-
+    This distribution allows one to specify coordinates whose x-, y-, and z-
    components are sampled independently from one another.

    Parameters
@ -304,7 +306,6 @@ class CartesianIndependent(Spatial):

    """

-
    def __init__(self, x, y, z):
        super().__init__()
        self.x = x
@ -375,6 +376,122 @@ class CartesianIndependent(Spatial):
        return cls(x, y, z)


+class SphericalIndependent(Spatial):
+    """Spatial distribution represented in spherical coordinates.
+
+    This distribution allows one to specify coordinates whose :math:`r`,
+    :math:`\theta`, and :math:`\phi` components are sampled independently from
+    one another and centered on the coordinates (x0, y0, z0).
+
+    Parameters
+    ----------
+    r : openmc.stats.Univariate
+        Distribution of r-coordinates
+    theta : openmc.stats.Univariate
+        Distribution of theta-coordinates (angle relative to the z-axis)
+    phi : openmc.stats.Univariate
+        Distribution of phi-coordinates (azimuthal angle)
+    origin: Iterable of float, optional
+        coordinates (x0, y0, z0) of the center of the sphere. Defaults to 
+        (0.0, 0.0, 0.0)
+
+    Attributes
+    ----------
+    r : openmc.stats.Univariate
+        Distribution of r-coordinates
+    theta : openmc.stats.Univariate
+        Distribution of theta-coordinates (angle relative to the z-axis)
+    phi : openmc.stats.Univariate
+        Distribution of phi-coordinates (azimuthal angle)
+    origin: Iterable of float, optional
+        coordinates (x0, y0, z0) of the center of the sphere. Defaults to 
+        (0.0, 0.0, 0.0)
+
+    """
+
+    def __init__(self, r, theta, phi, origin=(0.0, 0.0, 0.0)):
+        super().__init__()
+        self.r = r
+        self.theta = theta
+        self.phi = phi
+        self.origin = origin
+
+    @property
+    def r(self):
+        return self._r
+
+    @property
+    def theta(self):
+        return self._theta
+
+    @property
+    def phi(self):
+        return self._phi
+
+    @property
+    def origin(self):
+        return self._origin
+
+    @r.setter
+    def r(self, r):
+        cv.check_type('r coordinate', r, Univariate)
+        self._r = r
+
+    @theta.setter
+    def theta(self, theta):
+        cv.check_type('theta coordinate', theta, Univariate)
+        self._theta = theta
+
+    @phi.setter
+    def phi(self, phi):
+        cv.check_type('phi coordinate', phi, Univariate)
+        self._phi = phi
+
+    @origin.setter
+    def origin(self, origin):
+        cv.check_type('origin coordinates', origin, Iterable, Real)
+        origin = np.asarray(origin)
+        self._origin = origin
+
+    def to_xml_element(self):
+        """Return XML representation of the spatial distribution
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing spatial distribution data
+
+        """
+        element = ET.Element('space')
+        element.set('type', 'spherical')
+        element.append(self.r.to_xml_element('r'))
+        element.append(self.theta.to_xml_element('theta'))
+        element.append(self.phi.to_xml_element('phi'))
+        element.set("origin", ' '.join(map(str, self.origin)))
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate spatial distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.SphericalIndependent
+            Spatial distribution generated from XML element
+
+        """
+        r = Univariate.from_xml_element(elem.find('r'))
+        theta = Univariate.from_xml_element(elem.find('theta'))
+        phi = Univariate.from_xml_element(elem.find('phi'))
+        origin = [float(x) for x in elem.get('origin').split()]
+        return cls(r, theta, phi, origin=origin)
+
+
 class Box(Spatial):
    """Uniform distribution of coordinates in a rectangular cuboid.

--- a/openmc/summary.py
+++ b/openmc/summary.py
@ -97,16 +97,17 @@ class Summary(object):
                self._macroscopics = name.decode()

    def _read_geometry(self):
-        if "dagmc" in self._f['geometry'].attrs.keys():
-            return

-        # Read in and initialize the Materials and Geometry
+        # Read in and initialize the Materials
        self._read_materials()
-        self._read_surfaces()
-        cell_fills = self._read_cells()
-        self._read_universes()
-        self._read_lattices()
-        self._finalize_geometry(cell_fills)
+
+        # Read native geometry only
+        if "dagmc" not in self._f['geometry'].attrs.keys():
+            self._read_surfaces()
+            cell_fills = self._read_cells()
+            self._read_universes()
+            self._read_lattices()
+            self._finalize_geometry(cell_fills)

    def _read_materials(self):
        for group in self._f['materials'].values():
@ -153,8 +154,9 @@ class Summary(object):

                if 'rotation' in group:
                    rotation = group['rotation'][()]
-                    rotation = np.asarray(rotation, dtype=np.int)
-                    cell._rotation = rotation
+                    if rotation.size == 9:
+                        rotation.shape = (3, 3)
+                    cell.rotation = rotation

            elif fill_type == 'material':
                cell.temperature = group['temperature'][()]
--- a/openmc/universe.py
+++ b/openmc/universe.py
@ -418,7 +418,7 @@ class Universe(IDManagerMixin):

        return nuclides

-    def get_all_cells(self):
+    def get_all_cells(self, memo=None):
        """Return all cells that are contained within the universe

        Returns
@ -431,16 +431,22 @@ class Universe(IDManagerMixin):

        cells = OrderedDict()

+        if memo and self in memo:
+            return cells
+
+        if memo is not None:
+            memo.add(self)
+
        # Add this Universe's cells to the dictionary
        cells.update(self._cells)

        # Append all Cells in each Cell in the Universe to the dictionary
        for cell in self._cells.values():
-            cells.update(cell.get_all_cells())
+            cells.update(cell.get_all_cells(memo))

        return cells

-    def get_all_materials(self):
+    def get_all_materials(self, memo=None):
        """Return all materials that are contained within the universe

        Returns
@ -454,9 +460,9 @@ class Universe(IDManagerMixin):
        materials = OrderedDict()

        # Append all Cells in each Cell in the Universe to the dictionary
-        cells = self.get_all_cells()
+        cells = self.get_all_cells(memo)
        for cell in cells.values():
-            materials.update(cell.get_all_materials())
+            materials.update(cell.get_all_materials(memo))

        return materials

@ -512,19 +518,40 @@ class Universe(IDManagerMixin):

        return memo[self]

-    def create_xml_subelement(self, xml_element):
+    def create_xml_subelement(self, xml_element, memo=None):
+        """Add the universe xml representation to an incoming xml element
+
+        Parameters
+        ----------
+        xml_element : xml.etree.ElementTree.Element
+            XML element to be added to
+
+        memo : set or None
+            A set of object id's representing geometry entities already
+            written to the xml_element. This parameter is used internally
+            and should not be specified by users.
+
+        Returns
+        -------
+        None
+
+        """
        # Iterate over all Cells
        for cell_id, cell in self._cells.items():
-            path = "./cell[@id='{}']".format(cell_id)

-            # If the cell was not already written, write it
-            if xml_element.find(path) is None:
-                # Create XML subelement for this Cell
-                cell_element = cell.create_xml_subelement(xml_element)
+            # If the cell was already written, move on
+            if memo and cell in memo:
+                continue

-                # Append the Universe ID to the subelement and add to Element
-                cell_element.set("universe", str(self._id))
-                xml_element.append(cell_element)
+            if memo is not None:
+                memo.add(cell)
+
+            # Create XML subelement for this Cell
+            cell_element = cell.create_xml_subelement(xml_element, memo)
+
+            # Append the Universe ID to the subelement and add to Element
+            cell_element.set("universe", str(self._id))
+            xml_element.append(cell_element)

    def _determine_paths(self, path='', instances_only=False):
        """Count the number of instances for each cell in the universe, and
--- a/openmc/volume.py
+++ b/openmc/volume.py
@ -45,6 +45,8 @@ class VolumeCalculation(object):
        Lower-left coordinates of bounding box used to sample points
    upper_right : Iterable of float
        Upper-right coordinates of bounding box used to sample points
+    threshold : float
+        Threshold for the maximum standard deviation of volume in the calculation
    atoms : dict
        Dictionary mapping unique IDs of domains to a mapping of nuclides to
        total number of atoms for each nuclide present in the domain. For
@ -54,12 +56,20 @@ class VolumeCalculation(object):
        in each domain specified.
    volumes : dict
        Dictionary mapping unique IDs of domains to estimated volumes in cm^3.
+    threshold : float
+        Threshold for the maxmimum standard deviation of volumes.
+    trigger_type : {'variance', 'std_dev', 'rel_err'}
+        Value type used to halt volume calculation
+    iterations : int
+        Number of iterations over samples (for calculations with a trigger).

    """
-    def __init__(self, domains, samples, lower_left=None,
-                 upper_right=None):
+    def __init__(self, domains, samples, lower_left=None, upper_right=None):
        self._atoms = {}
        self._volumes = {}
+        self._threshold = None
+        self._trigger_type = None
+        self._iterations = None

        cv.check_type('domains', domains, Iterable,
                      (openmc.Cell, openmc.Material, openmc.Universe))
@ -123,6 +133,18 @@ class VolumeCalculation(object):
    def upper_right(self):
        return self._upper_right

+    @property
+    def threshold(self):
+        return self._threshold
+
+    @property
+    def trigger_type(self):
+        return self._trigger_type
+
+    @property
+    def iterations(self):
+        return self._iterations
+
    @property
    def domain_type(self):
        return self._domain_type
@ -170,6 +192,26 @@ class VolumeCalculation(object):
        cv.check_length(name, upper_right, 3)
        self._upper_right = upper_right

+    @threshold.setter
+    def threshold(self, threshold):
+        name = 'volume std. dev. threshold'
+        cv.check_type(name, threshold, Real)
+        cv.check_greater_than(name, threshold, 0.0)
+        self._threshold = threshold
+
+    @trigger_type.setter
+    def trigger_type(self, trigger_type):
+        cv.check_value('tally trigger type', trigger_type,
+                       ('variance', 'std_dev', 'rel_err'))
+        self._trigger_type = trigger_type
+
+    @iterations.setter
+    def iterations(self, iterations):
+        name = 'volume calculation iterations'
+        cv.check_type(name, iterations, Integral)
+        cv.check_greater_than(name, iterations, 0)
+        self._iterations = iterations
+
    @volumes.setter
    def volumes(self, volumes):
        cv.check_type('volumes', volumes, Mapping)
@ -180,6 +222,19 @@ class VolumeCalculation(object):
        cv.check_type('atoms', atoms, Mapping)
        self._atoms = atoms

+    def set_trigger(self, threshold, trigger_type):
+        """Set a trigger on the voulme calculation
+
+        Parameters
+        ----------
+        threshold : float
+            Threshold for the maxmimum standard deviation of volumes
+        trigger_type : {'variance', 'std_dev', 'rel_err'}
+            Value type used to halt volume calculation
+        """
+        self.trigger_type = trigger_type
+        self.threshold = threshold
+
    @classmethod
    def from_hdf5(cls, filename):
        """Load stochastic volume calculation results from HDF5 file.
@ -203,6 +258,10 @@ class VolumeCalculation(object):
            lower_left = f.attrs['lower_left']
            upper_right = f.attrs['upper_right']

+            threshold = f.attrs.get('threshold')
+            trigger_type = f.attrs.get('trigger_type')
+            iterations = f.attrs.get('iterations', 1)
+
            volumes = {}
            atoms = {}
            ids = []
@ -212,13 +271,12 @@ class VolumeCalculation(object):
                    ids.append(domain_id)
                    group = f[obj_name]
                    volume = ufloat(*group['volume'][()])
+                    volumes[domain_id] = volume
                    nucnames = group['nuclides'][()]
                    atoms_ = group['atoms'][()]
-
                    atom_dict = OrderedDict()
                    for name_i, atoms_i in zip(nucnames, atoms_):
                        atom_dict[name_i.decode()] = ufloat(*atoms_i)
-                    volumes[domain_id] = volume
                    atoms[domain_id] = atom_dict

        # Instantiate some throw-away domains that are used by the constructor
@ -234,6 +292,11 @@ class VolumeCalculation(object):

        # Instantiate the class and assign results
        vol = cls(domains, samples, lower_left, upper_right)
+
+        if trigger_type is not None:
+            vol.set_trigger(threshold, trigger_type.decode())
+
+        vol.iterations = iterations
        vol.volumes = volumes
        vol.atoms = atoms
        return vol
@ -278,4 +341,8 @@ class VolumeCalculation(object):
        ll_elem.text = ' '.join(str(x) for x in self.lower_left)
        ur_elem = ET.SubElement(element, "upper_right")
        ur_elem.text = ' '.join(str(x) for x in self.upper_right)
+        if self.threshold:
+            trigger_elem = ET.SubElement(element, "threshold")
+            trigger_elem.set("type", self.trigger_type)
+            trigger_elem.set("threshold", str(self.threshold))
        return element
--- a/scripts/openmc-plot-mesh-tally
+++ b/scripts/openmc-plot-mesh-tally
@ -176,8 +176,8 @@ class MeshPlotter(tk.Frame):

            # Set combobox items
            if filterType in ['Energy', 'Energyout']:
-                combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2])
-                                      for i in range(len(f.bins) - 1)]
+                combobox['values'] = ['{} to {}'.format(*ebin)
+                                      for ebin in f.bins]
            else:
                combobox['values'] = [str(i) for i in f.bins]

@ -210,7 +210,7 @@ class MeshPlotter(tk.Frame):
                continue
            elif f.short_name in ['Energy', 'Energyout']:
                index = self.filterBoxes[f.short_name].current()
-                ebin = (f.bins[index], f.bins[index + 1])
+                ebin = f.bins[index]
                spec_list.append((type(f), (ebin,)))
            else:
                index = self.filterBoxes[f.short_name].current()
--- a/scripts/openmc-update-inputs
+++ b/scripts/openmc-update-inputs
@ -298,4 +298,4 @@ if __name__ == '__main__':
            move(fname, fname + '.original')

            # Write a new geometry file.
-            tree.write(fname)
+            tree.write(fname, xml_declaration=True)
--- a/setup.py
+++ b/setup.py
@ -69,7 +69,7 @@ kwargs = {
        'pandas', 'lxml', 'uncertainties'
    ],
    'extras_require': {
-        'test': ['pytest', 'pytest-cov'],
+        'test': ['pytest', 'pytest-cov', 'colorama'],
        'vtk': ['vtk'],
    },
 }
--- a/src/cell.cpp
+++ b/src/cell.cpp
@ -1,8 +1,10 @@

 #include "openmc/cell.h"

+#include <algorithm>
 #include <cctype>
 #include <cmath>
+#include <iterator>
 #include <sstream>
 #include <set>
 #include <string>
@ -305,6 +307,10 @@ CSGCell::CSGCell(pugi::xml_node cell_node)

  if (fill_present) {
    fill_ = std::stoi(get_node_value(cell_node, "fill"));
+    if (fill_ == universe_) {
+      fatal_error("Cell " + std::to_string(id_) +
+        " is filled with the same universe that it is contained in.");
+    }
  } else {
    fill_ = C_NONE;
  }
@ -438,35 +444,39 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
    }

    auto rot {get_node_array<double>(cell_node, "rotation")};
-    if (rot.size() != 3) {
+    if (rot.size() != 3 && rot.size() != 9) {
      std::stringstream err_msg;
      err_msg << "Non-3D rotation vector applied to cell " << id_;
      fatal_error(err_msg);
    }

-    // Store the rotation angles.
-    rotation_.reserve(12);
-    rotation_.push_back(rot[0]);
-    rotation_.push_back(rot[1]);
-    rotation_.push_back(rot[2]);
-
    // Compute and store the rotation matrix.
-    auto phi = -rot[0] * PI / 180.0;
-    auto theta = -rot[1] * PI / 180.0;
-    auto psi = -rot[2] * PI / 180.0;
-    rotation_.push_back(std::cos(theta) * std::cos(psi));
-    rotation_.push_back(-std::cos(phi) * std::sin(psi)
-                        + std::sin(phi) * std::sin(theta) * std::cos(psi));
-    rotation_.push_back(std::sin(phi) * std::sin(psi)
-                        + std::cos(phi) * std::sin(theta) * std::cos(psi));
-    rotation_.push_back(std::cos(theta) * std::sin(psi));
-    rotation_.push_back(std::cos(phi) * std::cos(psi)
-                        + std::sin(phi) * std::sin(theta) * std::sin(psi));
-    rotation_.push_back(-std::sin(phi) * std::cos(psi)
-                        + std::cos(phi) * std::sin(theta) * std::sin(psi));
-    rotation_.push_back(-std::sin(theta));
-    rotation_.push_back(std::sin(phi) * std::cos(theta));
-    rotation_.push_back(std::cos(phi) * std::cos(theta));
+    rotation_.reserve(rot.size() == 9 ? 9 : 12);
+    if (rot.size() == 3) {
+      double phi = -rot[0] * PI / 180.0;
+      double theta = -rot[1] * PI / 180.0;
+      double psi = -rot[2] * PI / 180.0;
+      rotation_.push_back(std::cos(theta) * std::cos(psi));
+      rotation_.push_back(-std::cos(phi) * std::sin(psi)
+                          + std::sin(phi) * std::sin(theta) * std::cos(psi));
+      rotation_.push_back(std::sin(phi) * std::sin(psi)
+                          + std::cos(phi) * std::sin(theta) * std::cos(psi));
+      rotation_.push_back(std::cos(theta) * std::sin(psi));
+      rotation_.push_back(std::cos(phi) * std::cos(psi)
+                          + std::sin(phi) * std::sin(theta) * std::sin(psi));
+      rotation_.push_back(-std::sin(phi) * std::cos(psi)
+                          + std::cos(phi) * std::sin(theta) * std::sin(psi));
+      rotation_.push_back(-std::sin(theta));
+      rotation_.push_back(std::sin(phi) * std::cos(theta));
+      rotation_.push_back(std::cos(phi) * std::cos(theta));
+
+      // When user specifies angles, write them at end of vector
+      rotation_.push_back(rot[0]);
+      rotation_.push_back(rot[1]);
+      rotation_.push_back(rot[2]);
+    } else {
+      std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));
+    }
  }
 }

@ -578,8 +588,12 @@ CSGCell::to_hdf5(hid_t cell_group) const
      write_dataset(group, "translation", translation_);
    }
    if (!rotation_.empty()) {
-      std::array<double, 3> rot {rotation_[0], rotation_[1], rotation_[2]};
-      write_dataset(group, "rotation", rot);
+      if (rotation_.size() == 12) {
+        std::array<double, 3> rot {rotation_[9], rotation_[10], rotation_[11]};
+        write_dataset(group, "rotation", rot);
+      } else {
+        write_dataset(group, "rotation", rotation_);
+      }
    }

  } else if (type_ == FILL_LATTICE) {
--- a/src/cross_sections.cpp
+++ b/src/cross_sections.cpp
@ -155,7 +155,8 @@ void read_cross_sections_xml()
  if (settings::run_CE) {
    read_ce_cross_sections_xml();
  } else {
-    read_mg_cross_sections_header();
+    data::mg.read_header(settings::path_cross_sections);
+    put_mgxs_header_data_to_globals();
  }

  // Establish mapping between (type, material) and index in libraries
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@ -268,14 +268,12 @@ void load_dagmc_geometry()
    std::string cmp_str = mat_value;
    to_lower(cmp_str);

-    if (cmp_str.find("graveyard") != std::string::npos) {
+    if (cmp_str == "graveyard") {
      graveyard = vol_handle;
    }

    // material void checks
-    if (cmp_str.find("void") != std::string::npos   ||
-        cmp_str.find("vacuum") != std::string::npos ||
-        cmp_str.find("graveyard") != std::string::npos) {
+    if (cmp_str == "void" || cmp_str == "vacuum" || cmp_str == "graveyard") {
      c->material_.push_back(MATERIAL_VOID);
    } else {
      if (using_uwuw) {
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@ -61,6 +61,12 @@ Direction PolarAzimuthal::sample() const

  // Sample azimuthal angle
  double phi = phi_->sample();
+
+  // If the reference direction is along the z-axis, rotate the aziumthal angle
+  // to match spherical coordinate conventions.
+  // TODO: apply this change directly to rotate_angle
+  if (u_ref_.x == 0 && u_ref_.y == 0) phi += 0.5*PI;
+
  return rotate_angle(u_ref_, mu, &phi);
 }

--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@ -51,6 +51,73 @@ Position CartesianIndependent::sample() const
  return {x_->sample(), y_->sample(), z_->sample()};
 }

+//==============================================================================
+// SphericalIndependent implementation
+//==============================================================================
+
+SphericalIndependent::SphericalIndependent(pugi::xml_node node)
+{
+  // Read distribution for r-coordinate
+  if (check_for_node(node, "r")) {
+    pugi::xml_node node_dist = node.child("r");
+    r_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at r=0
+    double x[] {0.0};
+    double p[] {1.0};
+    r_ = std::make_unique<Discrete>(x, p, 1);
+  }
+
+  // Read distribution for theta-coordinate
+  if (check_for_node(node, "theta")) {
+    pugi::xml_node node_dist = node.child("theta");
+    theta_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at theta=0
+    double x[] {0.0};
+    double p[] {1.0};
+    theta_ = std::make_unique<Discrete>(x, p, 1);
+  }
+
+  // Read distribution for phi-coordinate
+  if (check_for_node(node, "phi")) {
+    pugi::xml_node node_dist = node.child("phi");
+    phi_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at phi=0
+    double x[] {0.0};
+    double p[] {1.0};
+    phi_ = std::make_unique<Discrete>(x, p, 1);
+  }
+
+  // Read sphere center coordinates
+  if (check_for_node(node, "origin")) {
+    auto origin = get_node_array<double>(node, "origin");
+    if (origin.size() == 3) {
+      origin_ = origin;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Origin for spherical source distribution must be length 3";
+      fatal_error(err_msg);
+    }
+  } else {
+    // If no coordinates were specified, default to (0, 0, 0)
+    origin_ = {0.0, 0.0, 0.0};
+  }
+
+}
+
+Position SphericalIndependent::sample() const
+{
+  double r = r_->sample();
+  double theta = theta_->sample();
+  double phi = phi_->sample();
+  double x = r*sin(theta)*cos(phi) + origin_.x;
+  double y =  r*sin(theta)*sin(phi) + origin_.y;
+  double z = r*cos(theta) + origin_.z;
+  return {x, y, z};
+}
+
 //==============================================================================
 // SpatialBox implementation
 //==============================================================================
--- a/Show more
+++ b/Show more