Merge pull request #746 from smharper/diff_tally6

Differential Tallies

docs/source/io_formats/statepoint.rst

@@ -162,6 +162,18 @@ if run_mode == 'k-eigenvalue':
 
     Width of each mesh cell in each dimension.
 
+**/tallies/derivatives/derivative <id>/independent variable** (*char[]*)
+
+    Independent variable of tally derivative
+
+**/tallies/derivatives/derivative <id>/material** (*int*)
+
+    ID of the perturbed material
+
+**/tallies/derivatives/derivative <id>/nuclide** (*char[]*)
+
+    Alias of the perturbed nuclide
+
 **/tallies/n_tallies** (*int*)
 
     Number of user-defined tallies.
@@ -204,6 +216,10 @@ if run_mode == 'k-eigenvalue':
     Array of nuclides to tally. Note that if no nuclide is specified in the user
     input, a single 'total' nuclide appears here.
 
+**/tallies/tally <uid>/derivative** (*int*)
+
+    ID of the derivative applied to the tally.
+
 **/tallies/tally <uid>/n_score_bins** (*int*)
 
     Number of scoring bins for a single nuclide. In general, this can be greater

docs/source/usersguide/input.rst

@@ -1916,6 +1916,13 @@ The ``<tally>`` element accepts the following sub-elements:
 
      *Default*: "all"
 
+  :derivative:
+    The id of a ``derivative`` element. This derivative will be applied to all
+    scores in the tally. Differential tallies are currently only implemented
+    for collision and analog estimators.
+
+     *Default*: None
+
 ``<mesh>`` Element
 ------------------
 
@@ -1944,6 +1951,39 @@ attributes/sub-elements:
       One of ``<upper_right>`` or ``<width>`` must be specified, but not both
       (even if they are consistent with one another).
 
+``<derivative>`` Element
+------------------------
+
+OpenMC can take the first-order derivative of many tallies with respect to
+material perturbations. It works by propagating a derivative through the
+transport equation. Essentially, OpenMC keeps track of how each particle's
+weight would change as materials are perturbed, and then accounts for that
+weight change in the tallies. Note that this assumes material perturbations are
+small enough not to change the distribution of fission sites. This element has
+the following attributes/sub-elements:
+
+  :id:
+    A unique integer that can be used to identify the derivative.
+
+  :variable:
+    The independent variable of the derivative. Accepted options are "density",
+    "nuclide_density", and "temperature". A "density" derivative will give the
+    derivative with respect to the density of the material in [g / cm^3]. A
+    "nuclide_density" derivative will give the derivative with respect to the
+    density of a particular nuclide in units of [atom / b / cm].  A
+    "temperature" derivative is with respect to a material temperature in units
+    of [K].  The temperature derivative requires windowed multipole to be
+    turned on.  Note also that the temperature derivative only accounts for
+    resolved resonance Doppler broadening.  It does not account for thermal
+    expansion, S(a, b) scattering, resonance scattering, or unresolved Doppler
+    broadening.
+
+  :material:
+    The perturbed material. (Necessary for all derivative types)
+
+  :nuclide:
+    The perturbed nuclide. (Necessary only for "nuclide_density")
+
 ``<assume_separate>`` Element
 -----------------------------
 

openmc/__init__.py

@@ -16,6 +16,7 @@ from openmc.universe import *
 from openmc.mesh import *
 from openmc.filter import *
 from openmc.trigger import *
+from openmc.tally_derivative import *
 from openmc.tallies import *
 from openmc.mgxs_library import *
 from openmc.cmfd import *

openmc/statepoint.py

@@ -97,6 +97,9 @@ class StatePoint(object):
         Dictionary whose keys are tally IDs and whose values are Tally objects
     tallies_present : bool
         Indicate whether user-defined tallies are present
+    tally_derivatives : dict
+        Dictionary whose keys are tally derivative IDs and whose values are
+        TallyDerivative objects
     version: tuple of Integral
         Version of OpenMC
     summary : None or openmc.Summary
@@ -129,6 +132,7 @@ class StatePoint(object):
         self._summary = None
         self._global_tallies = None
         self._sparse = False
+        self._derivs_read = False
 
         # Automatically link in a summary file if one exists
         if autolink:
@@ -391,6 +395,12 @@ class StatePoint(object):
                     base, tally_key)].value.decode()
                 tally.num_realizations = n_realizations
 
+                # Read derivative information.
+                if 'derivative' in self._f['{0}{1}'.format(base, tally_key)]:
+                    deriv_id = self._f['{0}{1}/derivative'.format(
+                                                        base, tally_key)].value
+                    tally.derivative = self.tally_derivatives[deriv_id]
+
                 # Read the number of Filters
                 n_filters = \
                     self._f['{0}{1}/n_filters'.format(base, tally_key)].value
@@ -452,6 +462,43 @@ class StatePoint(object):
     def tallies_present(self):
         return self._f['tallies/tallies_present'].value
 
+    @property
+    def tally_derivatives(self):
+        if not self._derivs_read:
+            # Initialize dictionaries for the Meshes
+            # Keys   - Derivative IDs
+            # Values - TallyDerivative objects
+            self._derivs = {}
+
+            # Populate the dictionary if any derivatives are present.
+            if 'derivatives' in self._f['tallies']:
+                # Read the derivative ids.
+                base = 'tallies/derivatives'
+                deriv_ids = [int(k.split(' ')[1]) for k in self._f[base]]
+
+                # Create each derivative object and add it to the dictionary.
+                for d_id in deriv_ids:
+                    base = 'tallies/derivatives/derivative {:d}'.format(d_id)
+                    deriv = openmc.TallyDerivative(derivative_id=d_id)
+                    deriv.variable = \
+                         self._f[base + '/independent variable'].value.decode()
+                    if deriv.variable == 'density':
+                        deriv.material = self._f[base + '/material'].value
+                    elif deriv.variable == 'nuclide_density':
+                        deriv.material = self._f[base + '/material'].value
+                        deriv.nuclide = \
+                             self._f[base + '/nuclide'].value.decode()
+                    elif deriv.variable == 'temperature':
+                        deriv.material = self._f[base + '/material'].value
+                    else:
+                        raise RuntimeError('Unrecognized tally differential '
+                                           'variable')
+                    self._derivs[d_id] = deriv
+
+            self._derivs_read = True
+
+        return self._derivs
+
     @property
     def version(self):
         return (self._f['version_major'].value,

openmc/tallies.py

@@ -104,6 +104,8 @@ class Tally(object):
     sparse : bool
         Whether or not the tally uses SciPy's LIL sparse matrix format for
         compressed data storage
+    derivative : openmc.TallyDerivative
+        A material perturbation derivative to apply to all scores in the tally.
 
     """
 
@@ -116,6 +118,7 @@ class Tally(object):
         self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores')
         self._estimator = None
         self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers')
+        self._derivative = None
 
         self._num_realizations = 0
         self._with_summary = False
@@ -151,6 +154,10 @@ class Tally(object):
             if nuclide not in other.nuclides:
                 return False
 
+        # Check derivatives
+        if self.derivative != other.derivative:
+            return False
+
         # Check all scores
         if len(self.scores) != len(other.scores):
             return False
@@ -172,27 +179,31 @@ class Tally(object):
 
     def __repr__(self):
         string = 'Tally\n'
-        string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
-        string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name)
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        string += '{: <16}=\t{}\n'.format('\tName', self.name)
 
-        string += '{0: <16}{1}\n'.format('\tFilters', '=\t')
+        if self.derivative is not None:
+            string += '{: <16}=\t{}\n'.format('\tDerivative ID',
+                                              str(self.derivative.id))
+
+        string += '{: <16}=\n'.format('\tFilters')
 
         for self_filter in self.filters:
-            string += '{0: <16}\t\t{1}\t{2}\n'.format('',
+            string += '{: <16}\t\t{}\t{}\n'.format('',
                 type(self_filter).__name__, self_filter.bins)
 
-        string += '{0: <16}{1}'.format('\tNuclides', '=\t')
+        string += '{: <16}=\t'.format('\tNuclides')
 
         for nuclide in self.nuclides:
             if isinstance(nuclide, openmc.Nuclide):
-                string += '{0} '.format(nuclide.name)
+                string += nuclide.name + ' '
             else:
-                string += '{0} '.format(nuclide)
+                string += nuclide + ' '
 
         string += '\n'
 
-        string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self.scores)
-        string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self.estimator)
+        string += '{: <16}=\t{}\n'.format('\tScores', self.scores)
+        string += '{: <16}=\t{}\n'.format('\tEstimator', self.estimator)
 
         return string
 
@@ -375,6 +386,10 @@ class Tally(object):
     def derived(self):
         return self._derived
 
+    @property
+    def derivative(self):
+        return self._derivative
+
     @property
     def sparse(self):
         return self._sparse
@@ -429,6 +444,12 @@ class Tally(object):
         else:
             self._name = ''
 
+    @derivative.setter
+    def derivative(self, deriv):
+        if deriv is not None:
+            cv.check_type('tally derivative', deriv, openmc.TallyDerivative)
+        self._derivative = deriv
+
     @filters.setter
     def filters(self, filters):
         cv.check_type('tally filters', filters, MutableSequence)
@@ -1055,6 +1076,11 @@ class Tally(object):
         for trigger in self.triggers:
             trigger.get_trigger_xml(element)
 
+        # Optional derivatives
+        if self.derivative is not None:
+            subelement = ET.SubElement(element, "derivative")
+            subelement.text = str(self.derivative.id)
+
         return element
 
     def contains_filter(self, filter_type):
@@ -1489,7 +1515,8 @@ class Tally(object):
         return data
 
     def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True,
-                             distribcell_paths=True, float_format='{:.2e}'):
+                             derivative=True, distribcell_paths=True,
+                             float_format='{:.2e}'):
         """Build a Pandas DataFrame for the Tally data.
 
         This method constructs a Pandas DataFrame object for the Tally data
@@ -1507,6 +1534,8 @@ class Tally(object):
             Include columns with nuclide bin information (default is True).
         scores : bool
             Include columns with score bin information (default is True).
+        derivative : bool
+            Include columns with differential tally info (default is True).
         distribcell_paths : bool, optional
             Construct columns for distribcell tally filters (default is True).
             The geometric information in the Summary object is embedded into a
@@ -1587,6 +1616,14 @@ class Tally(object):
             tile_factor = data_size / len(self.scores)
             df[column_name] = np.tile(scores, int(tile_factor))
 
+        # Include columns for derivatives if user requested it
+        if derivative and (self.derivative is not None):
+            df['d_variable'] = self.derivative.variable
+            if self.derivative.material is not None:
+                df['d_material'] = self.derivative.material
+            if self.derivative.nuclide is not None:
+                df['d_nuclide'] = self.derivative.nuclide
+
         # Append columns with mean, std. dev. for each tally bin
         df['mean'] = self.mean.ravel()
         df['std. dev.'] = self.std_dev.ravel()
@@ -3554,6 +3591,18 @@ class Tallies(cv.CheckedList):
                         self._tallies_file.append(xml_element)
                         already_written.add(f.mesh)
 
+    def _create_derivative_subelements(self):
+        # Get a list of all derivatives referenced in a tally.
+        derivs = []
+        for tally in self:
+            deriv = tally.derivative
+            if deriv is not None and deriv not in derivs:
+                derivs.append(deriv)
+
+        # Add the derivatives to the XML tree.
+        for d in derivs:
+            self._tallies_file.append(d.to_xml_element())
+
     def export_to_xml(self):
         """Create a tallies.xml file that can be used for a simulation.
 
@@ -3564,6 +3613,7 @@ class Tallies(cv.CheckedList):
 
         self._create_mesh_subelements()
         self._create_tally_subelements()
+        self._create_derivative_subelements()
 
         # Clean the indentation in the file to be user-readable
         clean_xml_indentation(self._tallies_file)

openmc/tally_derivative.py

@@ -0,0 +1,141 @@
+from __future__ import division
+
+import sys
+from numbers import Integral
+from xml.etree import ElementTree as ET
+
+from six import string_types
+
+import openmc.checkvalue as cv
+from openmc.mixin import EqualityMixin
+
+
+# "Static" variable for auto-generated TallyDerivative IDs
+AUTO_TALLY_DERIV_ID = 10000
+
+def reset_auto_tally_deriv_id():
+    global AUTO_TALLY_ID
+    AUTO_TALLY_DERIV_ID = 10000
+
+
+class TallyDerivative(EqualityMixin):
+    """A material perturbation derivative to apply to a tally.
+
+    Parameters
+    ----------
+    derivative_id : Integral, optional
+        Unique identifier for the tally derivative. If none is specified, an
+        identifier will automatically be assigned
+    variable : str, optional
+        Accepted values are 'density', 'nuclide_density', and 'temperature'
+    material : Integral, optional
+        The perturubed material ID
+    nuclide : str, optional
+        The perturbed nuclide. Only needed for 'nuclide_density' derivatives.
+        Ex: 'Xe135'
+
+    Attributes
+    ----------
+    id : Integral
+        Unique identifier for the tally derivative
+    variable : str
+        Accepted values are 'density', 'nuclide_density', and 'temperature'
+    material : Integral
+        The perturubed material ID
+    nuclide : str
+        The perturbed nuclide. Only needed for 'nuclide_density' derivatives.
+        Ex: 'Xe135'
+
+    """
+
+    def __init__(self, derivative_id=None, variable=None, material=None,
+                 nuclide=None):
+        # Initialize Tally class attributes
+        self.id = derivative_id
+        self.variable = variable
+        self.material = material
+        self.nuclide = nuclide
+
+    def __hash__(self):
+        return hash(repr(self))
+
+    def __repr__(self):
+        string = 'Tally Derivative\n'
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        string += '{: <16}=\t{}\n'.format('\tVariable', self.variable)
+
+        if self.variable == 'density':
+            string += '{: <16}=\t{}\n'.format('\tMaterial', self.material)
+        elif self.variable == 'nuclide_density':
+            string += '{: <16}=\t{}\n'.format('\tMaterial', self.material)
+            string += '{: <16}=\t{}\n'.format('\tNuclide', self.nuclide)
+        elif self.variable == 'temperature':
+            string += '{: <16}=\t{}\n'.format('\tMaterial', self.material)
+
+        return string
+
+    @property
+    def id(self):
+        return self._id
+
+    @property
+    def variable(self):
+        return self._variable
+
+    @property
+    def material(self):
+        return self._material
+
+    @property
+    def nuclide(self):
+        return self._nuclide
+
+    @id.setter
+    def id(self, deriv_id):
+        if deriv_id is None:
+            global AUTO_TALLY_DERIV_ID
+            self._id = AUTO_TALLY_DERIV_ID
+            AUTO_TALLY_DERIV_ID += 1
+        else:
+            cv.check_type('tally derivative ID', deriv_id, Integral)
+            cv.check_greater_than('tally derivative ID', deriv_id, 0,
+                                  equality=True)
+            self._id = deriv_id
+
+    @variable.setter
+    def variable(self, var):
+        if var is not None:
+            cv.check_type('derivative variable', var, string_types)
+            cv.check_value('derivative variable', var,
+                           ('density', 'nuclide_density', 'temperature'))
+        self._variable = var
+
+    @material.setter
+    def material(self, mat):
+        if mat is not None:
+            cv.check_type('derivative material', mat, Integral)
+        self._material = mat
+
+    @nuclide.setter
+    def nuclide(self, nuc):
+        if nuc is not None:
+            cv.check_type('derivative nuclide', nuc, string_types)
+        self._nuclide = nuc
+
+    def to_xml_element(self):
+        """Return XML representation of the tally derivative
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing derivative data
+
+        """
+
+        element = ET.Element("derivative")
+        element.set("id", str(self.id))
+        element.set("variable", self.variable)
+        element.set("material", str(self.material))
+        if self.variable == 'nuclide_density':
+            element.set("nuclide", self.nuclide)
+        return element

src/constants.F90

@@ -391,6 +391,12 @@ module constants
        K_TRACKLENGTH = 3, &
        LEAKAGE       = 4
 
+  ! Differential tally independent variables
+  integer, parameter :: &
+       DIFF_DENSITY = 1, &
+       DIFF_NUCLIDE_DENSITY = 2, &
+       DIFF_TEMPERATURE = 3
+
   ! ============================================================================
   ! RANDOM NUMBER STREAM CONSTANTS
 

src/cross_section.F90

@@ -7,7 +7,7 @@ module cross_section
   use global
   use list_header,      only: ListElemInt
   use material_header,  only: Material
-  use math,             only: faddeeva, broaden_wmp_polynomials
+  use math,             only: faddeeva, w_derivative, broaden_wmp_polynomials
   use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, &
                               MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,&
                               FIT_F, MultipoleArray
@@ -706,6 +706,94 @@ contains
     end if
   end subroutine multipole_eval
 
+!===============================================================================
+! MULTIPOLE_DERIV_EVAL evaluates the windowed multipole equations for the
+! derivative of cross sections in the resolved resonance regions with respect to
+! temperature.
+!===============================================================================
+
+  subroutine multipole_deriv_eval(multipole, E, sqrtkT, sigT, sigA, sigF)
+    type(MultipoleArray), intent(in) :: multipole ! The windowed multipole
+                                                  !  object to process.
+    real(8), intent(in)              :: E         ! The energy at which to
+                                                  !  evaluate the cross section
+    real(8), intent(in)              :: sqrtkT    ! The temperature in the form
+                                                  !  sqrt(kT), at which to
+                                                  !  evaluate the XS.
+    real(8), intent(out)             :: sigT      ! Total cross section
+    real(8), intent(out)             :: sigA      ! Absorption cross section
+    real(8), intent(out)             :: sigF      ! Fission cross section
+    complex(8) :: w_val    ! The faddeeva function evaluated at Z
+    complex(8) :: Z        ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole)
+    complex(8) :: sigT_factor(multipole % num_l)
+    real(8) :: sqrtE       ! sqrt(E), eV
+    real(8) :: invE        ! 1/E, eV
+    real(8) :: dopp        ! sqrt(atomic weight ratio / kT)
+    integer :: i_pole      ! index of pole
+    integer :: i_window    ! index of window
+    integer :: startw      ! window start pointer (for poles)
+    integer :: endw        ! window end pointer
+    real(8) :: T
+
+    ! ==========================================================================
+    ! Bookkeeping
+
+    ! Define some frequently used variables.
+    sqrtE = sqrt(E)
+    invE = ONE / E
+    dopp = multipole % sqrtAWR / sqrtkT
+    T = sqrtkT**2 / K_BOLTZMANN
+
+    if (sqrtkT == ZERO) call fatal_error("Windowed multipole temperature &
+         &derivatives are not implemented for 0 Kelvin cross sections.")
+
+    ! Locate us
+    i_window = floor((sqrtE - sqrt(multipole % start_E)) / multipole % spacing &
+         + ONE)
+    startw = multipole % w_start(i_window)
+    endw = multipole % w_end(i_window)
+
+    ! Fill in factors.
+    if (startw <= endw) then
+      call compute_sigT_factor(multipole, sqrtE, sigT_factor)
+    end if
+
+    ! Initialize the ouptut cross sections.
+    sigT = ZERO
+    sigA = ZERO
+    sigF = ZERO
+
+    ! TODO Polynomials: Some of the curvefit polynomials Doppler broaden so
+    ! rigorously we should be computing the derivative of those.  But in
+    ! practice, those derivatives are only large at very low energy and they
+    ! have no effect on reactor calculations.
+
+    ! ==========================================================================
+    ! Add the contribution from the poles in this window.
+
+    if (endw >= startw) then
+      do i_pole = startw, endw
+        Z = (sqrtE - multipole % data(MP_EA, i_pole)) * dopp
+        w_val = -invE * SQRT_PI * HALF * w_derivative(Z, 2)
+        if (multipole % formalism == FORM_MLBW) then
+          sigT = sigT + real((multipole % data(MLBW_RT, i_pole) * &
+                        sigT_factor(multipole%l_value(i_pole)) + &
+                        multipole % data(MLBW_RX, i_pole)) * w_val)
+          sigA = sigA + real(multipole % data(MLBW_RA, i_pole) * w_val)
+          sigF = sigF + real(multipole % data(MLBW_RF, i_pole) * w_val)
+        else if (multipole % formalism == FORM_RM) then
+          sigT = sigT + real(multipole % data(RM_RT, i_pole) * w_val * &
+                             sigT_factor(multipole % l_value(i_pole)))
+          sigA = sigA + real(multipole % data(RM_RA, i_pole) * w_val)
+          sigF = sigF + real(multipole % data(RM_RF, i_pole) * w_val)
+        end if
+      end do
+      sigT = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigT
+      sigA = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigA
+      sigF = -HALF*multipole % sqrtAWR / sqrt(K_BOLTZMANN) * T**(-1.5) * sigF
+    end if
+  end subroutine multipole_deriv_eval
+
 !===============================================================================
 ! COMPUTE_SIGT_FACTOR calculates the sigT_factor, a factor inside of the sigT
 ! equation not present in the sigA and sigF equations.

src/global.F90

@@ -20,7 +20,7 @@ module global
   use set_header,       only: SetInt
   use surface_header,   only: SurfaceContainer
   use source_header,    only: SourceDistribution
-  use tally_header,     only: TallyObject
+  use tally_header,     only: TallyObject, TallyDerivative
   use trigger_header,   only: KTrigger
   use timer_header,     only: Timer
   use volume_header,    only: VolumeCalculation
@@ -185,6 +185,10 @@ module global
   integer :: n_tallies      = 0 ! # of tallies
   integer :: n_user_tallies = 0 ! # of user tallies
 
+  ! Tally derivatives
+  type(TallyDerivative), allocatable :: tally_derivs(:)
+!$omp threadprivate(tally_derivs)
+
   ! Normalization for statistics
   integer :: n_realizations = 0 ! # of independent realizations
   real(8) :: total_weight       ! total starting particle weight in realization

src/input_xml.F90

@@ -2628,11 +2628,13 @@ contains
     type(Node), pointer :: node_mesh => null()
     type(Node), pointer :: node_tal => null()
     type(Node), pointer :: node_filt => null()
-    type(Node), pointer :: node_trigger=>null()
+    type(Node), pointer :: node_trigger => null()
+    type(Node), pointer :: node_deriv => null()
     type(NodeList), pointer :: node_mesh_list => null()
     type(NodeList), pointer :: node_tal_list => null()
     type(NodeList), pointer :: node_filt_list => null()
     type(NodeList), pointer :: node_trigger_list => null()
+    type(NodeList), pointer :: node_deriv_list => null()
     type(ElemKeyValueCI), pointer :: scores
     type(ElemKeyValueCI), pointer :: next
 
@@ -2640,6 +2642,12 @@ contains
     filename = trim(path_input) // "tallies.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
+      ! We need to allocate tally_derivs to avoid segfaults.  Also needs to be
+      ! done in parallel because tally derivs are threadprivate.
+!$omp parallel
+      allocate(tally_derivs(0))
+!$omp end parallel
+
       ! Since a tallies.xml file is optional, no error is issued here
       return
     end if
@@ -2818,6 +2826,94 @@ contains
     ! We only need the mesh info for plotting
     if (run_mode == MODE_PLOTTING) return
 
+    ! ==========================================================================
+    ! READ DATA FOR DERIVATIVES
+
+    ! Get pointer list to XML <derivative> nodes and allocate global array.
+    ! The array is threadprivate so it must be allocated in parallel.
+    call get_node_list(doc, "derivative", node_deriv_list)
+!$omp parallel
+    allocate(tally_derivs(get_list_size(node_deriv_list)))
+!$omp end parallel
+
+    ! Make sure this is not an MG run.
+    if (.not. run_CE .and. get_list_size(node_deriv_list) > 0) then
+      call fatal_error("Differential tallies not supported in multi-group mode")
+    end if
+
+    ! Read derivative attributes.
+    do i = 1, get_list_size(node_deriv_list)
+      associate(deriv => tally_derivs(i))
+        ! Get pointer to derivative node.
+        call get_list_item(node_deriv_list, i, node_deriv)
+
+        ! Copy the derivative id.
+        if (check_for_node(node_deriv, "id")) then
+          call get_node_value(node_deriv, "id", deriv % id)
+        else
+          call fatal_error("Must specify an ID for <derivative> elements in the&
+               & tally XML file")
+        end if
+
+        ! Make sure the id is > 0.
+        if (deriv % id <= 0) then
+          call fatal_error("<derivative> IDs must be an integer greater than &
+               &zero")
+        end if
+
+        ! Make sure this id has not already been used.
+        do j = 1, i-1
+          if (tally_derivs(j) % id == deriv % id) then
+            call fatal_error("Two or more <derivative>'s use the same unique &
+                 &ID: " // trim(to_str(deriv % id)))
+          end if
+        end do
+
+        ! Read the independent variable name.
+        temp_str = ""
+        call get_node_value(node_deriv, "variable", temp_str)
+        temp_str = to_lower(temp_str)
+
+        select case(temp_str)
+
+        case("density")
+          deriv % variable = DIFF_DENSITY
+          call get_node_value(node_deriv, "material", deriv % diff_material)
+
+        case("nuclide_density")
+          deriv % variable = DIFF_NUCLIDE_DENSITY
+          call get_node_value(node_deriv, "material", deriv % diff_material)
+
+          call get_node_value(node_deriv, "nuclide", word)
+          word = trim(to_lower(word))
+          pair_list => nuclide_dict % keys()
+          do while (associated(pair_list))
+            if (starts_with(pair_list % key, word)) then
+              word = pair_list % key(1:150)
+              exit
+            end if
+
+            ! Advance to next
+            pair_list => pair_list % next
+          end do
+
+          ! Check if no nuclide was found
+          if (.not. associated(pair_list)) then
+            call fatal_error("Could not find the nuclide " &
+                 // trim(word) // " specified in derivative " &
+                 // trim(to_str(deriv % id)) // " in any material.")
+          end if
+          deallocate(pair_list)
+
+          deriv % diff_nuclide = nuclide_dict % get_key(word)
+
+        case("temperature")
+          deriv % variable = DIFF_TEMPERATURE
+          call get_node_value(node_deriv, "material", deriv % diff_material)
+        end select
+      end associate
+    end do
+
     ! ==========================================================================
     ! READ TALLY DATA
 
@@ -2839,7 +2935,7 @@ contains
 
       t % estimator = ESTIMATOR_TRACKLENGTH
 
-      ! Copy material id
+      ! Copy tally id
       if (check_for_node(node_tal, "id")) then
         call get_node_value(node_tal, "id", t % id)
       else
@@ -3429,6 +3525,7 @@ contains
               call fatal_error("Cannot tally flux with an outgoing energy &
                    &filter.")
             end if
+
           case ('flux-yn')
             ! Prohibit user from tallying flux for an individual nuclide
             if (.not. (t % n_nuclide_bins == 1 .and. &
@@ -3766,6 +3863,48 @@ contains
              // trim(to_str(t % id)) // ".")
       end if
 
+      ! Check for a tally derivative.
+      if (check_for_node(node_tal, "derivative")) then
+        ! Temporarily store the derivative id.
+        call get_node_value(node_tal, "derivative", t % deriv)
+
+        ! Find the derivative with the given id, and store it's index.
+        do j = 1, size(tally_derivs)
+          if (tally_derivs(j) % id == t % deriv) then
+            t % deriv = j
+            ! Only analog or collision estimators are supported for differential
+            ! tallies.
+            if (t % estimator == ESTIMATOR_TRACKLENGTH) then
+              t % estimator = ESTIMATOR_COLLISION
+            end if
+            ! We found the derivative we were looking for; exit the do loop.
+            exit
+          end if
+          if (j == size(tally_derivs)) then
+            call fatal_error("Could not find derivative " &
+                 // trim(to_str(t % deriv)) // " specified on tally " &
+                 // trim(to_str(t % id)))
+          end if
+        end do
+
+        if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY &
+             .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then
+          if (any(t % nuclide_bins == -1)) then
+            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
+              call fatal_error("Error on tally " // trim(to_str(t % id)) &
+                   // ": Cannot use a 'nuclide_density' or 'temperature' &
+                   &derivative on a tally with an outgoing energy filter and &
+                   &'total' nuclide rate. Instead, tally each nuclide in the &
+                   &material individually.")
+              ! Note that diff tallies with these characteristics would work
+              ! correctly if no tally events occur in the perturbed material
+              ! (e.g. pertrubing moderator but only tallying fuel), but this
+              ! case would be hard to check for by only reading inputs.
+            end if
+          end if
+        end if
+      end if
+
       ! If settings.xml trigger is turned on, create tally triggers
       if (trigger_on) then
 
@@ -4787,8 +4926,8 @@ contains
         sum_percent = sum(mat % atom_density)
         mat % atom_density = mat % atom_density / sum_percent
 
-        ! Change density in g/cm^3 to atom/b-cm. Since all values are now in atom
-        ! percent, the sum needs to be re-evaluated as 1/sum(x*awr)
+        ! Change density in g/cm^3 to atom/b-cm. Since all values are now in
+        ! atom percent, the sum needs to be re-evaluated as 1/sum(x*awr)
         if (.not. density_in_atom) then
           sum_percent = ZERO
           do j = 1, mat % n_nuclides
@@ -4807,6 +4946,18 @@ contains
 
         ! Calculate nuclide atom densities
         mat % atom_density = mat % density * mat % atom_density
+
+        ! Calculate density in g/cm^3.
+        mat % density_gpcc = ZERO
+        do j = 1, mat % n_nuclides
+          if (run_CE) then
+            awr = nuclides(mat % nuclide(j)) % awr
+          else
+            awr = ONE
+          end if
+          mat % density_gpcc = mat % density_gpcc &
+               + mat % atom_density(j) * awr * MASS_NEUTRON / N_AVOGADRO
+        end do
       end associate
     end do
 

src/material_header.F90

@@ -13,6 +13,7 @@ module material_header
     integer, allocatable :: nuclide(:)      ! index in nuclides array
     real(8)              :: density         ! total atom density in atom/b-cm
     real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
+    real(8)              :: density_gpcc    ! total density in g/cm^3
 
     ! Energy grid information
     integer              :: n_grid    ! # of union material grid points

src/math.F90

@@ -753,6 +753,22 @@ contains
 
   end function faddeeva
 
+  recursive function w_derivative(z, order) result(wv)
+    complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
+    integer,                   intent(in) :: order
+    complex(8)     :: wv     ! The resulting w(z) value
+
+    select case(order)
+    case (0)
+      wv = faddeeva(z)
+    case (1)
+      wv = -TWO * z * faddeeva(z) + TWO * ONEI / SQRT_PI
+    case default
+      wv = -TWO * z * w_derivative(z, order-1) &
+           - TWO * (order-1) * w_derivative(z, order-2)
+    end select
+  end function w_derivative
+
 !===============================================================================
 ! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit.  The
 ! curvefit is a polynomial of the form

src/output.F90

@@ -777,6 +777,28 @@ contains
              // trim(t % name), unit=unit_tally, level=3)
       endif
 
+      ! Write derivative information.
+      if (t % deriv /= NONE) then
+        associate(deriv => tally_derivs(t % deriv))
+          select case (deriv % variable)
+          case (DIFF_DENSITY)
+            write(unit=unit_tally, fmt="(' Density derivative  Material ',A)") &
+                 to_str(deriv % diff_material)
+          case (DIFF_NUCLIDE_DENSITY)
+            write(unit=unit_tally, fmt="(' Nuclide density derivative  &
+                 &Material ',A,'  Nuclide ',A)") &
+                 trim(to_str(deriv % diff_material)), &
+                 trim(nuclides(deriv % diff_nuclide) % name)
+          case (DIFF_TEMPERATURE)
+            write(unit=unit_tally, fmt="(' Temperature derivative  Material ',&
+                 &A)") to_str(deriv % diff_material)
+          case default
+            call fatal_error("Differential tally dependent variable for tally "&
+                 // trim(to_str(t % id)) // " not defined in output.F90.")
+          end select
+        end associate
+      end if
+
       ! Handle surface current tallies separately
       if (t % type == TALLY_SURFACE_CURRENT) then
         call write_surface_current(t, unit_tally)

src/relaxng/tallies.rnc

@@ -15,6 +15,23 @@ element tallies {
     )
   }* &
 
+  element derivative {
+    (element id { xsd:int } | attribute id { xsd:int }) &
+    (element material { xsd:int } | attribute material { xsd:int }) &
+    ( (element variable { ( "density") }
+        | attribute variable { ( "density" ) } ) |
+      (
+        (element variable { ( "nuclide_density" ) }
+          | attribute variable { ( "nuclide_density" ) } )
+        &
+        (element nuclide { xsd:string { maxLength = "12" } }
+          | attribute nuclide { xsd:string { maxLength = "12" } } ) |
+      )
+      (element variable { ( "temperature") }
+        | attribute variable { ( "temperature" ) } )
+    )
+  }* &
+
   element tally {
     (element id { xsd:int } | attribute id { xsd:int }) &
     (element name { xsd:string { maxLength="52" } } |
@@ -41,7 +58,8 @@ element tallies {
       (element type { xsd:string } | attribute type { xsd:string }) &
       (element threshold { xsd:double} | attribute threshold { xsd:double }) &
       (element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )?
-    }*
+    }* &
+    (element derivative { xsd:int } | attribute derivative { xsd:int } )?
   }* &
 
   element assume_separate { xsd:boolean }?

src/relaxng/tallies.rng

@@ -89,6 +89,68 @@
         </interleave>
       </element>
     </zeroOrMore>
+    <zeroOrMore>
+      <element name="derivative">
+        <interleave>
+          <choice>
+            <element name="id">
+              <data type="int"/>
+            </element>
+            <attribute name="id">
+              <data type="int"/>
+            </attribute>
+          </choice>
+          <choice>
+            <element name="material">
+              <data type="int"/>
+            </element>
+            <attribute name="material">
+              <data type="int"/>
+            </attribute>
+          </choice>
+          <choice>
+            <choice>
+              <element name="variable">
+                <value>density</value>
+              </element>
+              <attribute name="variable">
+                <value>density</value>
+              </attribute>
+            </choice>
+            <interleave>
+              <choice>
+                <element name="variable">
+                  <value>nuclide_density</value>
+                </element>
+                <attribute name="variable">
+                  <value>nuclide_density</value>
+                </attribute>
+              </choice>
+              <choice>
+                <element name="nuclide">
+                  <data type="string">
+                    <param name="maxLength">12</param>
+                  </data>
+                </element>
+                <attribute name="nuclide">
+                  <data type="string">
+                    <param name="maxLength">12</param>
+                  </data>
+                </attribute>
+              </choice>
+            </interleave>
+            <choice>
+              <element name="variable">
+                <value>temperature</value>
+              </element>
+              <attribute name="variable">
+                <value>temperature</value>
+              </attribute>
+            </choice>
+          </choice>
+        </interleave>
+      </element>
+    </zeroOrMore>
     <zeroOrMore>
       <element name="tally">
         <interleave>
@@ -254,6 +316,16 @@
               </interleave>
             </element>
           </zeroOrMore>
+          <optional>
+            <choice>
+              <element name="derivative">
+                <data type="int"/>
+              </element>
+              <attribute name="derivative">
+                <data type="int"/>
+              </attribute>
+            </choice>
+          </optional>
         </interleave>
       </element>
     </zeroOrMore>

src/simulation.F90

@@ -83,7 +83,7 @@ contains
 
         ! ====================================================================
         ! LOOP OVER PARTICLES
-!$omp parallel do schedule(static) firstprivate(p)
+!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs)
         PARTICLE_LOOP: do i_work = 1, work
           current_work = i_work
 

src/state_point.F90

@@ -50,7 +50,8 @@ contains
     integer, allocatable :: key_array(:)
     integer(HID_T) :: file_id
     integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
-                      mesh_group, filter_group, runtime_group
+                      mesh_group, filter_group, derivs_group, deriv_group, &
+                      runtime_group
     character(MAX_WORD_LEN), allocatable :: str_array(:)
     character(MAX_FILE_LEN)    :: filename
     type(RegularMesh), pointer :: meshp
@@ -198,6 +199,38 @@ contains
 
       call close_group(meshes_group)
 
+      ! Write information for derivatives.
+      if (size(tally_derivs) > 0) then
+        derivs_group = create_group(tallies_group, "derivatives")
+        do i = 1, size(tally_derivs)
+          associate(deriv => tally_derivs(i))
+            deriv_group = create_group(derivs_group, "derivative " &
+                 // trim(to_str(deriv % id)))
+            select case (deriv % variable)
+            case (DIFF_DENSITY)
+              call write_dataset(deriv_group, "independent variable", "density")
+              call write_dataset(deriv_group, "material", deriv % diff_material)
+            case (DIFF_NUCLIDE_DENSITY)
+              call write_dataset(deriv_group, "independent variable", &
+                   "nuclide_density")
+              call write_dataset(deriv_group, "material", deriv % diff_material)
+              call write_dataset(deriv_group, "nuclide", &
+                   nuclides(deriv % diff_nuclide) % name)
+            case (DIFF_TEMPERATURE)
+              call write_dataset(deriv_group, "independent variable", &
+                   "temperature")
+              call write_dataset(deriv_group, "material", deriv % diff_material)
+            case default
+              call fatal_error("Independent variable for derivative " &
+                   // trim(to_str(deriv % id)) // " not defined in &
+                   &state_point.F90.")
+            end select
+            call close_group(deriv_group)
+          end associate
+        end do
+        call close_group(derivs_group)
+      end if
+
       ! Write number of tallies
       call write_dataset(tallies_group, "n_tallies", n_tallies)
 
@@ -273,6 +306,13 @@ contains
           call write_dataset(tally_group, "nuclides", str_array)
           deallocate(str_array)
 
+          ! Write derivative information.
+          if (tally % deriv /= NONE) then
+            call write_dataset(tally_group, "derivative", &
+                 tally_derivs(tally % deriv) % id)
+          end if
+
+          ! Write scores.
           call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
           allocate(str_array(size(tally % score_bins)))
           do j = 1, size(tally % score_bins)

src/tally.F90

@@ -8,6 +8,7 @@ module tally
 
   use algorithm,        only: binary_search
   use constants
+  use cross_section,    only: multipole_deriv_eval
   use error,            only: fatal_error
   use geometry_header
   use global
@@ -1087,6 +1088,14 @@ contains
 
       end select
 
+      !#########################################################################
+      ! Add derivative information on score for differential tallies.
+
+      if (t % deriv /= NONE) then
+        call apply_derivative_to_score(p, t, i_nuclide, atom_density, &
+             score_bin, score)
+      end if
+
       !#########################################################################
       ! Expand score if necessary and add to tally results.
       call expand_and_score(p, t, score_index, filter_index, score_bin, &
@@ -2420,6 +2429,13 @@ contains
         ! determine score based on bank site weight and keff
         score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
 
+        ! Add derivative information for differential tallies.  Note that the
+        ! i_nuclide and atom_density arguments do not matter since this is an
+        ! analog estimator.
+        if (t % deriv /= NONE) then
+          call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score)
+        end if
+
         if (.not. run_CE .and. eo_filt % matches_transport_groups) then
 
           ! determine outgoing energy from fission bank
@@ -3107,6 +3123,778 @@ contains
 
   end subroutine score_surface_current
 
+!===============================================================================
+! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative
+!===============================================================================
+
+  subroutine apply_derivative_to_score(p, t, i_nuclide, atom_density, &
+                                       score_bin, score)
+    type(Particle),    intent(in)    :: p
+    type(TallyObject), intent(in)    :: t
+    integer,           intent(in)    :: i_nuclide
+    real(8),           intent(in)    :: atom_density   ! atom/b-cm
+    integer,           intent(in)    :: score_bin
+    real(8),           intent(inout) :: score
+
+    integer :: l
+    logical :: scoring_diff_nuclide
+    real(8) :: flux_deriv
+    real(8) :: dsigT, dsigA, dsigF, cum_dsig
+
+    if (score == ZERO) return
+
+    ! If our score was previously c then the new score is
+    ! c * (1/f * d_f/d_p + 1/c * d_c/d_p)
+    ! where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the
+    ! perturbated variable.
+
+    associate(deriv => tally_derivs(t % deriv))
+      flux_deriv = deriv % flux_deriv
+
+      select case (tally_derivs(t % deriv) % variable)
+
+      !=========================================================================
+      ! Density derivative:
+      ! c = Sigma_MT
+      ! c = sigma_MT * N
+      ! c = sigma_MT * rho * const
+      ! d_c / d_rho = sigma_MT * const
+      ! (1 / c) * (d_c / d_rho) = 1 / rho
+
+      case (DIFF_DENSITY)
+        select case (t % estimator)
+
+        case (ESTIMATOR_ANALOG)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
+                SCORE_NU_FISSION)
+            if (materials(p % material) % id == deriv % diff_material) then
+              score = score * (flux_deriv + ONE &
+                   / materials(p % material) % density_gpcc)
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case (ESTIMATOR_COLLISION)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
+                SCORE_NU_FISSION)
+            if (materials(p % material) % id == deriv % diff_material) then
+              score = score * (flux_deriv + ONE &
+                   / materials(p % material) % density_gpcc)
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case default
+            call fatal_error("Differential tallies are only implemented for &
+                 &analog and collision estimators.")
+        end select
+
+      !=========================================================================
+      ! Nuclide density derivative:
+      ! If we are scoring a reaction rate for a single nuclide then
+      ! c = Sigma_MT_i
+      ! c = sigma_MT_i * N_i
+      ! d_c / d_N_i = sigma_MT_i
+      ! (1 / c) * (d_c / d_N_i) = 1 / N_i
+      ! If the score is for the total material (i_nuclide = -1)
+      ! c = Sum_i(Sigma_MT_i)
+      ! d_c / d_N_i = sigma_MT_i
+      ! (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT
+      ! where i is the perturbed nuclide.
+
+      case (DIFF_NUCLIDE_DENSITY)
+        select case (t % estimator)
+
+        case (ESTIMATOR_ANALOG)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, &
+                SCORE_NU_FISSION)
+            if (materials(p % material) % id == deriv % diff_material &
+                 .and. p % event_nuclide == deriv % diff_nuclide) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == deriv % diff_nuclide) exit
+                end do
+
+                score = score * (flux_deriv &
+                     + ONE / mat % atom_density(l))
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case (ESTIMATOR_COLLISION)
+          scoring_diff_nuclide = &
+               (materials(p % material) % id == deriv % diff_material) &
+               .and. (i_nuclide == deriv % diff_nuclide)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % total /= ZERO) then
+              score = score * (flux_deriv &
+                   + micro_xs(deriv % diff_nuclide) % total &
+                   / material_xs % total)
+            else if (scoring_diff_nuclide .and. &
+                 micro_xs(deriv % diff_nuclide) % total /= ZERO) then
+              score = score * (flux_deriv + ONE / atom_density)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_SCATTER)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % total - material_xs % absorption /= ZERO) then
+              score = score * (flux_deriv &
+                   + (micro_xs(deriv % diff_nuclide) % total &
+                   - micro_xs(deriv % diff_nuclide) % absorption) &
+                   / (material_xs % total - material_xs % absorption))
+            else if (scoring_diff_nuclide .and. &
+                 (micro_xs(deriv % diff_nuclide) % total &
+                 - micro_xs(deriv % diff_nuclide) % absorption) /= ZERO) then
+              score = score * (flux_deriv + ONE / atom_density)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_ABSORPTION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % absorption /= ZERO) then
+              score = score * (flux_deriv &
+                   + micro_xs(deriv % diff_nuclide) % absorption &
+                   / material_xs % absorption )
+            else if (scoring_diff_nuclide .and. &
+                 micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then
+              score = score * (flux_deriv + ONE / atom_density)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_FISSION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % fission /= ZERO) then
+              score = score * (flux_deriv &
+                   + micro_xs(deriv % diff_nuclide) % fission &
+                   / material_xs % fission)
+            else if (scoring_diff_nuclide .and. &
+                 micro_xs(deriv % diff_nuclide) % fission /= ZERO) then
+              score = score * (flux_deriv + ONE / atom_density)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_NU_FISSION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % nu_fission /= ZERO) then
+              score = score * (flux_deriv &
+                   + micro_xs(deriv % diff_nuclide) % nu_fission &
+                   / material_xs % nu_fission)
+            else if (scoring_diff_nuclide .and. &
+                 micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then
+              score = score * (flux_deriv + ONE / atom_density)
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case default
+            call fatal_error("Differential tallies are only implemented for &
+                 &analog and collision estimators.")
+        end select
+
+      !=========================================================================
+      ! Temperature derivative:
+      ! If we are scoring a reaction rate for a single nuclide then
+      ! c = Sigma_MT_i
+      ! c = sigma_MT_i * N_i
+      ! d_c / d_T = (d_sigma_Mt_i / d_T) * N_i
+      ! (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i
+      ! If the score is for the total material (i_nuclide = -1)
+      ! (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i
+      ! where i is the perturbed nuclide.  The d_sigma_MT_i / d_T term is
+      ! computed by multipole_deriv_eval.  It only works for the resolved
+      ! resonance range and requires multipole data.
+
+      case (DIFF_TEMPERATURE)
+        select case (t % estimator)
+
+        case (ESTIMATOR_ANALOG)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL)
+            if (materials(p % material) % id == deriv % diff_material .and. &
+                 micro_xs(p % event_nuclide) % total > ZERO) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == p % event_nuclide) exit
+                end do
+
+                dsigT = ZERO
+                associate (nuc => nuclides(p % event_nuclide))
+                  if (nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                  end if
+                end associate
+                score = score * (flux_deriv &
+                     + dsigT * mat % atom_density(l) / material_xs % total)
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_SCATTER)
+            if (materials(p % material) % id == deriv % diff_material .and. &
+                 (micro_xs(p % event_nuclide) % total &
+                 - micro_xs(p % event_nuclide) % absorption) > ZERO) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == p % event_nuclide) exit
+                end do
+
+                dsigT = ZERO
+                dsigA = ZERO
+                associate (nuc => nuclides(p % event_nuclide))
+                  if (nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                  end if
+                end associate
+                score = score * (flux_deriv + (dsigT - dsigA) &
+                     * mat % atom_density(l) / &
+                     (material_xs % total - material_xs % absorption))
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_ABSORPTION)
+            if (materials(p % material) % id == deriv % diff_material .and. &
+                 micro_xs(p % event_nuclide) % absorption > ZERO) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == p % event_nuclide) exit
+                end do
+
+                dsigA = ZERO
+                associate (nuc => nuclides(p % event_nuclide))
+                  if (nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                  end if
+                end associate
+                score = score * (flux_deriv &
+                     + dsigA * mat % atom_density(l) / material_xs % absorption)
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_FISSION)
+            if (materials(p % material) % id == deriv % diff_material .and. &
+                 micro_xs(p % event_nuclide) % fission > ZERO) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == p % event_nuclide) exit
+                end do
+
+                dsigF = ZERO
+                associate (nuc => nuclides(p % event_nuclide))
+                  if (nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                  end if
+                end associate
+                score = score * (flux_deriv &
+                     + dsigF * mat % atom_density(l) / material_xs % fission)
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_NU_FISSION)
+            if (materials(p % material) % id == deriv % diff_material .and. &
+                 micro_xs(p % event_nuclide) % nu_fission > ZERO) then
+              associate(mat => materials(p % material))
+                ! Search for the index of the perturbed nuclide.
+                do l = 1, mat % n_nuclides
+                  if (mat % nuclide(l) == p % event_nuclide) exit
+                end do
+
+                dsigF = ZERO
+                associate (nuc => nuclides(p % event_nuclide))
+                  if (nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                  end if
+                end associate
+                score = score * (flux_deriv &
+                     + dsigF * mat % atom_density(l) / material_xs % nu_fission&
+                     * micro_xs(p % event_nuclide) % nu_fission &
+                     / micro_xs(p % event_nuclide) % fission)
+              end associate
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case (ESTIMATOR_COLLISION)
+
+          select case (score_bin)
+
+          case (SCORE_FLUX)
+            score = score * flux_deriv
+
+          case (SCORE_TOTAL)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % total > ZERO) then
+              cum_dsig = ZERO
+              associate(mat => materials(p % material))
+                do l = 1, mat % n_nuclides
+                  associate (nuc => nuclides(mat % nuclide(l)))
+                    if (nuc % mp_present .and. &
+                         p % last_E >= nuc % multipole % start_E .and. &
+                         p % last_E <= nuc % multipole % end_E .and. &
+                         micro_xs(mat % nuclide(l)) % total > ZERO) then
+                      call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                           p % sqrtkT, dsigT, dsigA, dsigF)
+                      cum_dsig = cum_dsig + dsigT * mat % atom_density(l)
+                    end if
+                  end associate
+                end do
+              end associate
+              score = score * (flux_deriv &
+                   + cum_dsig / material_xs % total)
+            else if (materials(p % material) % id == deriv % diff_material &
+                 .and. material_xs % total > ZERO) then
+              dsigT = ZERO
+              associate (nuc => nuclides(i_nuclide))
+                if (nuc % mp_present .and. &
+                     p % last_E >= nuc % multipole % start_E .and. &
+                     p % last_E <= nuc % multipole % end_E) then
+                  call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                       p % sqrtkT, dsigT, dsigA, dsigF)
+                end if
+              end associate
+              score = score * (flux_deriv &
+                   + dsigT / micro_xs(i_nuclide) % total)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_SCATTER)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 (material_xs % total - material_xs % absorption) > ZERO) then
+              cum_dsig = ZERO
+              associate(mat => materials(p % material))
+                do l = 1, mat % n_nuclides
+                  associate (nuc => nuclides(mat % nuclide(l)))
+                    if (nuc % mp_present .and. &
+                         p % last_E >= nuc % multipole % start_E .and. &
+                         p % last_E <= nuc % multipole % end_E .and. &
+                         (micro_xs(mat % nuclide(l)) % total &
+                         - micro_xs(mat % nuclide(l)) % absorption) > ZERO) then
+                      call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                           p % sqrtkT, dsigT, dsigA, dsigF)
+                      cum_dsig = cum_dsig &
+                           + (dsigT - dsigA) * mat % atom_density(l)
+                    end if
+                  end associate
+                end do
+              end associate
+              score = score * (flux_deriv + cum_dsig &
+                   / (material_xs % total - material_xs % absorption))
+            else if ( materials(p % material) % id == deriv % diff_material &
+                 .and. (material_xs % total - material_xs % absorption) > ZERO)&
+                 then
+              dsigT = ZERO
+              dsigA = ZERO
+              associate (nuc => nuclides(i_nuclide))
+                if (nuc % mp_present .and. &
+                     p % last_E >= nuc % multipole % start_E .and. &
+                     p % last_E <= nuc % multipole % end_E) then
+                  call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                       p % sqrtkT, dsigT, dsigA, dsigF)
+                end if
+              end associate
+              score = score * (flux_deriv + (dsigT - dsigA) &
+                   / (micro_xs(i_nuclide) % total &
+                   - micro_xs(i_nuclide) % absorption))
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_ABSORPTION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % absorption > ZERO) then
+              cum_dsig = ZERO
+              associate(mat => materials(p % material))
+                do l = 1, mat % n_nuclides
+                  associate (nuc => nuclides(mat % nuclide(l)))
+                    if (nuc % mp_present .and. &
+                         p % last_E >= nuc % multipole % start_E .and. &
+                         p % last_E <= nuc % multipole % end_E .and. &
+                         micro_xs(mat % nuclide(l)) % absorption > ZERO) then
+                      call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                           p % sqrtkT, dsigT, dsigA, dsigF)
+                      cum_dsig = cum_dsig + dsigA * mat % atom_density(l)
+                    end if
+                  end associate
+                end do
+              end associate
+              score = score * (flux_deriv &
+                   + cum_dsig / material_xs % absorption)
+            else if (materials(p % material) % id == deriv % diff_material &
+                 .and. material_xs % absorption > ZERO) then
+              dsigA = ZERO
+              associate (nuc => nuclides(i_nuclide))
+                if (nuc % mp_present .and. &
+                     p % last_E >= nuc % multipole % start_E .and. &
+                     p % last_E <= nuc % multipole % end_E) then
+                  call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                       p % sqrtkT, dsigT, dsigA, dsigF)
+                end if
+              end associate
+              score = score * (flux_deriv &
+                   + dsigA / micro_xs(i_nuclide) % absorption)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_FISSION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % fission > ZERO) then
+              cum_dsig = ZERO
+              associate(mat => materials(p % material))
+                do l = 1, mat % n_nuclides
+                  associate (nuc => nuclides(mat % nuclide(l)))
+                    if (nuc % mp_present .and. &
+                         p % last_E >= nuc % multipole % start_E .and. &
+                         p % last_E <= nuc % multipole % end_E .and. &
+                         micro_xs(mat % nuclide(l)) % fission > ZERO) then
+                      call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                           p % sqrtkT, dsigT, dsigA, dsigF)
+                      cum_dsig = cum_dsig + dsigF * mat % atom_density(l)
+                    end if
+                  end associate
+                end do
+              end associate
+              score = score * (flux_deriv &
+                   + cum_dsig / material_xs % fission)
+            else if (materials(p % material) % id == deriv % diff_material &
+                 .and. material_xs % fission > ZERO) then
+              dsigF = ZERO
+              associate (nuc => nuclides(i_nuclide))
+                if (nuc % mp_present .and. &
+                     p % last_E >= nuc % multipole % start_E .and. &
+                     p % last_E <= nuc % multipole % end_E) then
+                  call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                       p % sqrtkT, dsigT, dsigA, dsigF)
+                end if
+              end associate
+              score = score * (flux_deriv &
+                   + dsigF / micro_xs(i_nuclide) % fission)
+            else
+              score = score * flux_deriv
+            end if
+
+          case (SCORE_NU_FISSION)
+            if (i_nuclide == -1 .and. &
+                 materials(p % material) % id == deriv % diff_material .and. &
+                 material_xs % nu_fission > ZERO) then
+              cum_dsig = ZERO
+              associate(mat => materials(p % material))
+                do l = 1, mat % n_nuclides
+                  associate (nuc => nuclides(mat % nuclide(l)))
+                    if (nuc % mp_present .and. &
+                         p % last_E >= nuc % multipole % start_E .and. &
+                         p % last_E <= nuc % multipole % end_E .and. &
+                         micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then
+                      call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                           p % sqrtkT, dsigT, dsigA, dsigF)
+                      cum_dsig = cum_dsig + dsigF * mat % atom_density(l) &
+                           * micro_xs(mat % nuclide(l)) % nu_fission &
+                           / micro_xs(mat % nuclide(l)) % fission
+                    end if
+                  end associate
+                end do
+              end associate
+              score = score * (flux_deriv &
+                   + cum_dsig / material_xs % nu_fission)
+            else if (materials(p % material) % id == deriv % diff_material &
+                 .and. material_xs % nu_fission > ZERO) then
+              dsigF = ZERO
+              associate (nuc => nuclides(i_nuclide))
+                if (nuc % mp_present .and. &
+                     p % last_E >= nuc % multipole % start_E .and. &
+                     p % last_E <= nuc % multipole % end_E) then
+                  call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                       p % sqrtkT, dsigT, dsigA, dsigF)
+                end if
+              end associate
+              score = score * (flux_deriv &
+                   + dsigF / micro_xs(i_nuclide) % fission)
+            else
+              score = score * flux_deriv
+            end if
+
+          case default
+            call fatal_error('Tally derivative not defined for a score on &
+                 &tally ' // trim(to_str(t % id)))
+          end select
+
+        case default
+            call fatal_error("Differential tallies are only implemented for &
+                 &analog and collision estimators.")
+        end select
+      end select
+    end associate
+  end subroutine apply_derivative_to_score
+
+!===============================================================================
+! SCORE_TRACK_DERIVATIVE Adjust flux derivatives on differential tallies to
+! account for a neutron travelling through a perturbed material.
+!===============================================================================
+
+  subroutine score_track_derivative(p, distance)
+    type(Particle), intent(in) :: p
+    real(8),        intent(in) :: distance ! Neutron flight distance
+
+    integer :: i, l
+    real(8) :: dsigT, dsigA, dsigF
+
+    ! A void material cannot be perturbed so it will not affect flux derivatives
+    if (p % material == MATERIAL_VOID) return
+
+    do i = 1, size(tally_derivs)
+      associate(deriv => tally_derivs(i))
+        select case (deriv % variable)
+
+        case (DIFF_DENSITY)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material) then
+              ! phi is proportional to e^(-Sigma_tot * dist)
+              ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
+              ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
+              deriv % flux_deriv = deriv % flux_deriv &
+                   - distance * material_xs % total / mat % density_gpcc
+            end if
+          end associate
+
+        case (DIFF_NUCLIDE_DENSITY)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material) then
+              ! phi is proportional to e^(-Sigma_tot * dist)
+              ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
+              ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
+              deriv % flux_deriv = deriv % flux_deriv &
+                   - distance * micro_xs(deriv % diff_nuclide) % total
+            end if
+          end associate
+
+        case (DIFF_TEMPERATURE)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material) then
+              do l=1, mat % n_nuclides
+                associate (nuc => nuclides(mat % nuclide(l)))
+                  if (nuc % mp_present .and. &
+                       p % E >= nuc % multipole % start_E .and. &
+                       p % E <= nuc % multipole % end_E) then
+                    ! phi is proportional to e^(-Sigma_tot * dist)
+                    ! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
+                    ! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
+                    call multipole_deriv_eval(nuc % multipole, p % E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                    deriv % flux_deriv = deriv % flux_deriv &
+                         - distance * dsigT * mat % atom_density(l)
+                  end if
+                end associate
+              end do
+            end if
+          end associate
+        end select
+      end associate
+    end do
+  end subroutine score_track_derivative
+
+!===============================================================================
+! SCORE_COLLISION_DERIVATIVE Adjust flux derivatives on differential tallies to
+! account for a neutron scattering in the perturbed material.  Note that this
+! subroutine will be called after absorption events in addition to scattering
+! events, but any flux derivatives scored after an absorption will never be
+! tallied.  This is because the order of operations is
+! 1. Particle is moved.
+! 2. score_track_derivative is called.
+! 3. Collision physics are computed, and the particle is labeled absorbed.
+! 4. Analog- and collision-estimated tallies are scored.
+! 5. This subroutine is called.
+! 6. Particle is killed and no more tallies are scored.
+! Hence, it is safe to assume that only derivative of the scattering cross
+! section need to be computed here.
+!===============================================================================
+
+  subroutine score_collision_derivative(p)
+    type(Particle), intent(in) :: p
+
+    integer :: i, j, l
+    real(8) :: dsigT, dsigA, dsigF
+
+    ! A void material cannot be perturbed so it will not affect flux derivatives
+    if (p % material == MATERIAL_VOID) return
+
+    do i = 1, size(tally_derivs)
+      associate(deriv => tally_derivs(i))
+        select case (deriv % variable)
+
+        case (DIFF_DENSITY)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material) then
+              ! phi is proportional to Sigma_s
+              ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
+              ! (1 / phi) * (d_phi / d_rho) = 1 / rho
+              deriv % flux_deriv = deriv % flux_deriv &
+                   + ONE / mat % density_gpcc
+            end if
+          end associate
+
+        case (DIFF_NUCLIDE_DENSITY)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material &
+                 .and. p % event_nuclide == deriv % diff_nuclide) then
+              ! Find the index in this material for the diff_nuclide.
+              do j = 1, mat % n_nuclides
+                if (mat % nuclide(j) == deriv % diff_nuclide) exit
+              end do
+              ! Make sure we found the nuclide.
+              if (mat % nuclide(j) /= deriv % diff_nuclide) then
+                call fatal_error("Couldn't find the right nuclide.")
+              end if
+              ! phi is proportional to Sigma_s
+              ! (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s
+              ! (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s
+              ! (1 / phi) * (d_phi / d_N) = 1 / N
+              deriv % flux_deriv = deriv % flux_deriv &
+                   + ONE / mat % atom_density(j)
+            end if
+          end associate
+
+        case (DIFF_TEMPERATURE)
+          associate (mat => materials(p % material))
+            if (mat % id == deriv % diff_material) then
+              do l=1, mat % n_nuclides
+                associate (nuc => nuclides(mat % nuclide(l)))
+                  if (mat % nuclide(l) == p % event_nuclide .and. &
+                       nuc % mp_present .and. &
+                       p % last_E >= nuc % multipole % start_E .and. &
+                       p % last_E <= nuc % multipole % end_E) then
+                    ! phi is proportional to Sigma_s
+                    ! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
+                    ! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
+                    call multipole_deriv_eval(nuc % multipole, p % last_E, &
+                         p % sqrtkT, dsigT, dsigA, dsigF)
+                    deriv % flux_deriv = deriv % flux_deriv + (dsigT - dsigA)&
+                         / (micro_xs(mat % nuclide(l)) % total &
+                         - micro_xs(mat % nuclide(l)) % absorption)
+                    ! Note that this is an approximation!  The real scattering
+                    ! cross section is Sigma_s(E'->E, uvw'->uvw) =
+                    ! Sigma_s(E') * P(E'->E, uvw'->uvw).  We are assuming that
+                    ! d_P(E'->E, uvw'->uvw) / d_T = 0 and only computing
+                    ! d_S(E') / d_T.  Using this approximation in the vicinity
+                    ! of low-energy resonances causes errors (~2-5% for PWR
+                    ! pincell eigenvalue derivatives).
+                  end if
+                end associate
+              end do
+            end if
+          end associate
+        end select
+      end associate
+    end do
+  end subroutine score_collision_derivative
+
+!===============================================================================
+! ZERO_FLUX_DERIVS Set the flux derivatives on differential tallies to zero.
+!===============================================================================
+
+  subroutine zero_flux_derivs()
+    integer :: i
+    do i = 1, size(tally_derivs)
+      tally_derivs(i) % flux_deriv = ZERO
+    end do
+  end subroutine zero_flux_derivs
+
 !===============================================================================
 ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
 ! within the batch to a new random variable

src/tally_header.F90

@@ -10,6 +10,19 @@ module tally_header
 
   implicit none
 
+!===============================================================================
+! TALLYDERIVATIVE describes a first-order derivative that can be applied to
+! tallies.
+!===============================================================================
+
+  type TallyDerivative
+    integer :: id
+    integer :: variable
+    integer :: diff_material
+    integer :: diff_nuclide
+    real(8) :: flux_deriv
+  end type TallyDerivative
+
 !===============================================================================
 ! TALLYOBJECT describes a user-specified tally. The region of phase space to
 ! tally in is given by the TallyFilters and the results are stored in a
@@ -72,6 +85,9 @@ module tally_header
     integer                           :: n_triggers = 0  ! # of triggers
     type(TriggerObject),  allocatable :: triggers(:)     ! Array of triggers
 
+    ! Index for the TallyDerivative for differential tallies.
+    integer :: deriv = NONE
+
   contains
     procedure :: write_results_hdf5
     procedure :: read_results_hdf5

src/tracking.F90

@@ -14,7 +14,9 @@ module tracking
   use random_lcg,         only: prn
   use string,             only: to_str
   use tally,              only: score_analog_tally, score_tracklength_tally, &
-                                score_collision_tally, score_surface_current
+                                score_collision_tally, score_surface_current, &
+                                score_track_derivative, &
+                                score_collision_derivative, zero_flux_derivs
   use track_output,       only: initialize_particle_track, write_particle_track, &
                                 add_particle_track, finalize_particle_track
 
@@ -63,6 +65,9 @@ contains
       call initialize_particle_track()
     endif
 
+    ! Every particle starts with no accumulated flux derivative.
+    if (active_tallies % size() > 0) call zero_flux_derivs()
+
     EVENT_LOOP: do
       ! If the cell hasn't been determined based on the particle's location,
       ! initiate a search for the current cell. This generally happens at the
@@ -137,6 +142,9 @@ contains
              distance * material_xs % nu_fission
       end if
 
+      ! Score flux derivative accumulators for differential tallies.
+      if (active_tallies % size() > 0) call score_track_derivative(p, distance)
+
       if (d_collision > d_boundary) then
         ! ====================================================================
         ! PARTICLE CROSSES SURFACE
@@ -213,6 +221,9 @@ contains
             p % coord(j + 1) % uvw = p % coord(j) % uvw
           end if
         end do
+
+        ! Score flux derivative accumulators for differential tallies.
+        if (active_tallies % size() > 0) call score_collision_derivative(p)
       end if
 
       ! Save coordinates for tallying purposes

tests/test_diff_tally/inputs_true.dat

@@ -0,0 +1 @@
+df5a0ee7d353fe25344496058cea42e3244a7ab32c13c34c03b4b2f6adf88579f08a71ec0d22d8d264e50e3663717068ff54308e2d16aed8e085c222a8c2fdbd

tests/test_diff_tally/results_true.dat

@@ -0,0 +1,177 @@
+d_material,d_nuclide,d_variable,score,mean,std. dev.
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,3.9902949e-02,9.1258428e-03
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,9.8213745e-04,9.8213745e-04
+1,,density,nu-fission,2.7945389e-02,2.0999368e-02
+1,,density,scatter,4.0626155e-01,1.7017674e-02
+1,,density,nu-fission,2.4595279e-02,2.1224443e-02
+1,,density,scatter,2.6505962e-03,2.3698970e-03
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,-1.2721728e-01,1.8810986e-01
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,0.0000000e+00,0.0000000e+00
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,4.1950592e-02,7.3680059e-02
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,scatter,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,scatter,-1.9798576e-02,1.9798576e-02
+1,O16,nuclide_density,nu-fission,9.3632921e-01,2.0322286e+00
+1,O16,nuclide_density,scatter,-2.2042881e-01,1.6781767e-01
+1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,scatter,2.3761276e+00,2.3761276e+00
+1,U235,nuclide_density,nu-fission,7.7441690e+02,8.6460933e+01
+1,U235,nuclide_density,scatter,9.6917103e+01,1.0750582e+01
+1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
+1,,temperature,scatter,-2.0607998e-06,2.0607998e-06
+1,,temperature,nu-fission,3.8845881e-06,3.5405083e-05
+1,,temperature,scatter,-6.6772431e-07,2.0750826e-07
+1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
+1,,temperature,scatter,0.0000000e+00,0.0000000e+00
+1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
+1,,temperature,scatter,0.0000000e+00,0.0000000e+00
+3,,density,flux,-7.6179774e+00,5.0326447e+00
+3,,density,flux,-1.6242019e+01,6.6231774e+00
+1,,density,flux,-2.2394746e-01,5.4534297e-02
+1,,density,flux,-5.7132854e-02,1.6898134e-01
+1,O16,nuclide_density,flux,-1.3403658e+01,1.3017127e+01
+1,O16,nuclide_density,flux,-1.0499715e+01,2.1684953e+01
+1,U235,nuclide_density,flux,-2.4741168e+03,5.7986812e+01
+1,U235,nuclide_density,flux,-1.9955533e+03,4.3254820e+02
+1,,temperature,flux,1.0601815e-04,3.3076550e-04
+1,,temperature,flux,1.6664450e-04,2.8761112e-04
+3,,density,total,-3.3161585e+00,2.5615932e+00
+3,,density,absorption,-4.3491248e-01,5.1559021e-01
+3,,density,scatter,-2.8812460e+00,2.0540305e+00
+3,,density,fission,-2.3060366e-01,3.1191109e-01
+3,,density,nu-fission,-5.6817321e-01,7.6143389e-01
+3,,density,total,-2.8908007e-01,3.9383184e-01
+3,,density,absorption,-2.5237485e-01,3.6565053e-01
+3,,density,scatter,-3.6705221e-02,2.9275371e-02
+3,,density,fission,-2.1271120e-01,3.0873698e-01
+3,,density,nu-fission,-5.1866339e-01,7.5235533e-01
+3,,density,total,2.7320569e+00,8.7709465e+00
+3,,density,absorption,8.4322466e-02,1.2807898e-01
+3,,density,scatter,2.6477344e+00,8.6432567e+00
+3,,density,fission,0.0000000e+00,0.0000000e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,total,0.0000000e+00,0.0000000e+00
+3,,density,absorption,0.0000000e+00,0.0000000e+00
+3,,density,scatter,0.0000000e+00,0.0000000e+00
+3,,density,fission,0.0000000e+00,0.0000000e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,total,4.7523439e-01,2.5652004e-02
+1,,density,absorption,3.1301234e-02,7.6793142e-03
+1,,density,scatter,4.4393315e-01,1.8023252e-02
+1,,density,fission,1.5559355e-02,3.2310294e-03
+1,,density,nu-fission,3.8658109e-02,7.8012233e-03
+1,,density,total,2.2062939e-02,4.3689061e-03
+1,,density,absorption,1.6845626e-02,4.0750601e-03
+1,,density,scatter,5.2173132e-03,2.9663452e-04
+1,,density,fission,1.3503649e-02,3.2642875e-03
+1,,density,nu-fission,3.2945474e-02,7.9507685e-03
+1,,density,total,-9.4735691e-02,2.5584665e-01
+1,,density,absorption,-4.0246397e-03,5.2111770e-03
+1,,density,scatter,-9.0711052e-02,2.5073762e-01
+1,,density,fission,0.0000000e+00,0.0000000e+00
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,density,total,0.0000000e+00,0.0000000e+00
+1,,density,absorption,0.0000000e+00,0.0000000e+00
+1,,density,scatter,0.0000000e+00,0.0000000e+00
+1,,density,fission,0.0000000e+00,0.0000000e+00
+1,,density,nu-fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,total,4.4488215e+01,5.3981603e+00
+1,O16,nuclide_density,absorption,-2.9372611e-01,9.3282813e-01
+1,O16,nuclide_density,scatter,4.4781941e+01,4.5257555e+00
+1,O16,nuclide_density,fission,9.9827872e-02,4.0698927e-01
+1,O16,nuclide_density,nu-fission,2.3342880e-01,9.8749212e-01
+1,O16,nuclide_density,total,4.0914051e-02,6.2301812e-01
+1,O16,nuclide_density,absorption,8.8289354e-02,5.5570922e-01
+1,O16,nuclide_density,scatter,-4.7375303e-02,6.7714071e-02
+1,O16,nuclide_density,fission,1.3678733e-01,4.4006413e-01
+1,O16,nuclide_density,nu-fission,3.3257520e-01,1.0719313e+00
+1,O16,nuclide_density,total,-1.6055918e+01,2.3439335e+01
+1,O16,nuclide_density,absorption,-3.9601231e-01,3.0131622e-01
+1,O16,nuclide_density,scatter,-1.5659906e+01,2.3140407e+01
+1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,total,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,absorption,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,fission,0.0000000e+00,0.0000000e+00
+1,O16,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,total,-6.1444902e+02,6.0111692e+01
+1,U235,nuclide_density,absorption,2.3977134e+02,4.5868405e+01
+1,U235,nuclide_density,scatter,-8.5422036e+02,2.4877920e+01
+1,U235,nuclide_density,fission,3.1319132e+02,3.1465126e+01
+1,U235,nuclide_density,nu-fission,7.6415465e+02,7.6501604e+01
+1,U235,nuclide_density,total,5.0916179e+02,3.5499319e+01
+1,U235,nuclide_density,absorption,3.9477360e+02,3.4326165e+01
+1,U235,nuclide_density,scatter,1.1438818e+02,2.7378217e+00
+1,U235,nuclide_density,fission,3.1360739e+02,3.2243605e+01
+1,U235,nuclide_density,nu-fission,7.6527228e+02,7.8683697e+01
+1,U235,nuclide_density,total,-4.1815711e+03,8.7796611e+02
+1,U235,nuclide_density,absorption,-1.0563226e+02,2.5308627e+01
+1,U235,nuclide_density,scatter,-4.0759388e+03,8.5332185e+02
+1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,total,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,absorption,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,scatter,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,fission,0.0000000e+00,0.0000000e+00
+1,U235,nuclide_density,nu-fission,0.0000000e+00,0.0000000e+00
+1,,temperature,total,2.1367130e-04,1.8632816e-04
+1,,temperature,absorption,6.9456249e-05,3.2199390e-05
+1,,temperature,scatter,1.4421505e-04,1.5754180e-04
+1,,temperature,fission,3.0674981e-06,1.8048377e-05
+1,,temperature,nu-fission,7.4768877e-06,4.3976766e-05
+1,,temperature,total,5.8223116e-06,2.3786474e-05
+1,,temperature,absorption,5.2142698e-06,2.1901911e-05
+1,,temperature,scatter,6.0804185e-07,1.9834890e-06
+1,,temperature,fission,3.0674200e-06,1.8048542e-05
+1,,temperature,nu-fission,7.4767028e-06,4.3977153e-05
+1,,temperature,total,2.1510156e-04,4.6388707e-04
+1,,temperature,absorption,2.2409527e-06,8.2901060e-06
+1,,temperature,scatter,2.1286061e-04,4.5560034e-04
+1,,temperature,fission,0.0000000e+00,0.0000000e+00
+1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
+1,,temperature,total,0.0000000e+00,0.0000000e+00
+1,,temperature,absorption,0.0000000e+00,0.0000000e+00
+1,,temperature,scatter,0.0000000e+00,0.0000000e+00
+1,,temperature,fission,0.0000000e+00,0.0000000e+00
+1,,temperature,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,absorption,-1.6517201e-01,3.0984738e-01
+3,,density,absorption,8.0402344e-03,1.4689335e-01
+1,,density,absorption,2.9069882e-02,2.2139609e-03
+1,,density,absorption,-9.4690065e-03,8.6605392e-03
+1,O16,nuclide_density,absorption,7.6911962e-01,4.1945687e-01
+1,O16,nuclide_density,absorption,-6.3724795e-01,6.6341488e-01
+1,U235,nuclide_density,absorption,1.4109543e+02,1.8518890e+01
+1,U235,nuclide_density,absorption,-1.2052822e+02,3.2084002e+01
+1,,temperature,absorption,3.9995404e-05,2.1703384e-05
+1,,temperature,absorption,2.7274361e-06,8.9339257e-06
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,-8.7934551e-01,1.0210652e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,-3.5236516e-03,3.5236516e-03
+3,,density,nu-fission,7.2953012e-02,2.9725863e-01
+3,,density,scatter,-2.2716410e+00,1.2546347e+00
+3,,density,nu-fission,8.9171481e-02,3.0634856e-01
+3,,density,scatter,-1.0151747e-03,9.6296775e-03
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,2.3163923e+00,4.8766690e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,0.0000000e+00,0.0000000e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,4.0762434e-01,3.8504141e+00
+3,,density,nu-fission,0.0000000e+00,0.0000000e+00
+3,,density,scatter,0.0000000e+00,0.0000000e+00

tests/test_diff_tally/test_diff_tally.py

@@ -0,0 +1,140 @@
+#!/usr/bin/env python
+
+import glob
+import os
+import sys
+
+import pandas as pd
+
+sys.path.insert(0, os.pardir)
+from testing_harness import PyAPITestHarness
+import openmc
+
+class DiffTallyTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        # Build default materials/geometry
+        self._input_set.build_default_materials_and_geometry()
+
+        # Set settings explicitly
+        self._input_set.settings.batches = 3
+        self._input_set.settings.inactive = 0
+        self._input_set.settings.particles = 100
+        self._input_set.settings.source = openmc.Source(space=openmc.stats.Box(
+            [-160, -160, -183], [160, 160, 183]))
+        self._input_set.settings.temperature['multipole'] = True
+
+        self._input_set.tallies = openmc.Tallies()
+
+        filt_mats = openmc.MaterialFilter((1, 3))
+        filt_eout = openmc.EnergyoutFilter((0.0, 0.625, 20.0e6))
+
+        # We want density derivatives for both water and fuel to get coverage
+        # for both fissile and non-fissile materials.
+        d1 = openmc.TallyDerivative(derivative_id=1)
+        d1.variable = 'density'
+        d1.material = 3
+        d2 = openmc.TallyDerivative(derivative_id=2)
+        d2.variable = 'density'
+        d2.material = 1
+
+        # O-16 is a good nuclide to test against because it is present in both
+        # water and fuel.  Some routines need to recognize that they have the
+        # perturbed nuclide but not the perturbed material.
+        d3 = openmc.TallyDerivative(derivative_id=3)
+        d3.variable = 'nuclide_density'
+        d3.material = 1
+        d3.nuclide = 'O16'
+
+        # A fissile nuclide, just for good measure.
+        d4 = openmc.TallyDerivative(derivative_id=4)
+        d4.variable = 'nuclide_density'
+        d4.material = 1
+        d4.nuclide = 'U235'
+
+        # Temperature derivatives.
+        d5 = openmc.TallyDerivative(derivative_id=5)
+        d5.variable = 'temperature'
+        d5.material = 1
+
+        derivs = [d1, d2, d3, d4, d5]
+
+        # Cover the flux score.
+        for i in range(5):
+            t = openmc.Tally()
+            t.add_score('flux')
+            t.add_filter(filt_mats)
+            t.derivative = derivs[i]
+            self._input_set.tallies.append(t)
+
+        # Cover supported scores with a collision estimator.
+        for i in range(5):
+            t = openmc.Tally()
+            t.add_score('total')
+            t.add_score('absorption')
+            t.add_score('scatter')
+            t.add_score('fission')
+            t.add_score('nu-fission')
+            t.add_filter(filt_mats)
+            t.add_nuclide('total')
+            t.add_nuclide('U235')
+            t.derivative = derivs[i]
+            self._input_set.tallies.append(t)
+
+        # Cover an analog estimator.
+        for i in range(5):
+            t = openmc.Tally()
+            t.add_score('absorption')
+            t.add_filter(filt_mats)
+            t.estimator = 'analog'
+            t.derivative = derivs[i]
+            self._input_set.tallies.append(t)
+
+        # Energyout filter and total nuclide for the density derivatives.
+        for i in range(2):
+            t = openmc.Tally()
+            t.add_score('nu-fission')
+            t.add_score('scatter')
+            t.add_filter(filt_mats)
+            t.add_filter(filt_eout)
+            t.add_nuclide('total')
+            t.add_nuclide('U235')
+            t.derivative = derivs[i]
+            self._input_set.tallies.append(t)
+
+        # Energyout filter without total nuclide for other derivatives.
+        for i in range(2, 5):
+            t = openmc.Tally()
+            t.add_score('nu-fission')
+            t.add_score('scatter')
+            t.add_filter(filt_mats)
+            t.add_filter(filt_eout)
+            t.add_nuclide('U235')
+            t.derivative = derivs[i]
+            self._input_set.tallies.append(t)
+
+        self._input_set.export()
+
+    def _get_results(self):
+        # Read the statepoint and summary files.
+        statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
+        sp = openmc.StatePoint(statepoint)
+
+        # Extract the tally data as a Pandas DataFrame.
+        df = pd.DataFrame()
+        for t in sp.tallies.values():
+            df = df.append(t.get_pandas_dataframe(), ignore_index=True)
+
+        # Extract the relevant data as a CSV string.
+        cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean',
+                'std. dev.')
+        return df.to_csv(None, columns=cols, index=False, float_format='%.7e')
+
+    def _cleanup(self):
+        super(DiffTallyTestHarness, self)._cleanup()
+        f = os.path.join(os.getcwd(), 'tallies.xml')
+        if os.path.exists(f): os.remove(f)
+
+
+if __name__ == '__main__':
+    harness = DiffTallyTestHarness('statepoint.3.h5', True)
+    harness.main()