diff --git a/src/constants.F90 b/src/constants.F90
index 375c517e76..8e8e2db76c 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -259,7 +259,7 @@ module constants
        EVENT_ABSORB  =  2
 
   ! Tally score type
-  integer, parameter :: N_SCORE_TYPES = 22
+  integer, parameter :: N_SCORE_TYPES = 23
   integer, parameter :: &
        SCORE_FLUX               = -1,  & ! flux
        SCORE_TOTAL              = -2,  & ! total reaction rate
@@ -282,7 +282,8 @@ module constants
        SCORE_NU_SCATTER_YN      = -19, & ! angular flux-weighted nu-scattering moment (0:N)
        SCORE_EVENTS             = -20, & ! number of events
        SCORE_DELAYED_NU_FISSION = -21, & ! delayed neutron production rate
-       SCORE_INVERSE_VELOCITY   = -22    ! flux-weighted inverse velocity
+       SCORE_INVERSE_VELOCITY   = -22, & ! flux-weighted inverse velocity
+       SCORE_KEFF               = -23    ! multiplication factor
 
   ! Maximum scattering order supported
   integer, parameter :: MAX_ANG_ORDER = 10
@@ -348,6 +349,11 @@ module constants
        K_TRACKLENGTH = 3, &
        LEAKAGE       = 4
 
+  ! Differential tally dependent variables
+  integer, parameter :: &
+       DIFF_DENSITY = 1, &
+       DIFF_NUCLIDE_DENSITY = 2
+
   ! ============================================================================
   ! RANDOM NUMBER STREAM CONSTANTS
 
diff --git a/src/initialize.F90 b/src/initialize.F90
index 9d63eb4e93..6130be0e25 100644
--- a/src/initialize.F90
+++ b/src/initialize.F90
@@ -831,6 +831,7 @@ contains
       ! 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
+        mat % density_gpcc = -mat % density
         sum_percent = ZERO
         do j = 1, mat%n_nuclides
           index_list = xs_listing_dict%get_key(mat%names(j))
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 58e4013355..052748d139 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2197,6 +2197,7 @@ contains
     type(Node), pointer :: node_tal => null()
     type(Node), pointer :: node_filt => 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()
@@ -2407,7 +2408,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
@@ -2843,6 +2844,52 @@ contains
         end do
       end if
 
+      ! =======================================================================
+      ! READ DATA FOR DERIVATIVES
+
+      if (check_for_node(node_tal, "derivative")) then
+        call get_node_ptr(node_tal, "derivative", node_deriv)
+        allocate(t % deriv)
+
+        temp_str = ""
+        call get_node_value(node_deriv, "variable", temp_str)
+        temp_str = to_lower(temp_str)
+
+        select case(temp_str)
+
+        case("density")
+          t % deriv % dep_var = DIFF_DENSITY
+          call get_node_value(node_deriv, "material", t % deriv % diff_material)
+
+        case("nuclide_density")
+          t % deriv % dep_var = DIFF_NUCLIDE_DENSITY
+          call get_node_value(node_deriv, "material", t % 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 tally " &
+                 &// trim(to_str(t % id)) // " in any material.")
+          end if
+          deallocate(pair_list)
+
+          t % deriv % diff_nuclide = nuclide_dict % get_key(word)
+        end select
+      end if
+
       ! =======================================================================
       ! READ DATA FOR SCORES
 
@@ -2991,6 +3038,11 @@ contains
               call fatal_error("Cannot tally flux with an outgoing energy &
                    &filter.")
             end if
+
+            if (allocated(t % deriv)) then
+              t % estimator = ESTIMATOR_COLLISION
+            end if
+
           case ('flux-yn')
             ! Prohibit user from tallying flux for an individual nuclide
             if (.not. (t % n_nuclide_bins == 1 .and. &
@@ -3188,6 +3240,12 @@ contains
           case ('events')
             t % score_bins(j) = SCORE_EVENTS
 
+          case ('keff')
+            t % score_bins(j) = SCORE_KEFF
+            if (allocated(t % deriv)) then
+              t % estimator = ESTIMATOR_COLLISION
+            end if
+
           case ('elastic', '(n,elastic)')
             t % score_bins(j) = ELASTIC
           case ('(n,2nd)')
diff --git a/src/material_header.F90 b/src/material_header.F90
index c2f4b51bc0..24e459cbb5 100644
--- a/src/material_header.F90
+++ b/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
diff --git a/src/output.F90 b/src/output.F90
index 55cd5b2a5b..0f667c121b 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -986,6 +986,7 @@ contains
     score_names(abs(SCORE_NU_SCATTER_YN))      = "Scattering Prod. Rate Moment"
     score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
     score_names(abs(SCORE_INVERSE_VELOCITY))   = "Flux-Weighted Inverse Velocity"
+    score_names(abs(SCORE_KEFF))               = "k_eff, multiplaction factor"
 
     ! Create filename for tally output
     filename = trim(path_output) // "tallies.out"
diff --git a/src/tally.F90 b/src/tally.F90
index fcd6dfdd9b..7fce90cce5 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -654,6 +654,13 @@ contains
           end if
         end if
 
+      case (SCORE_KEFF)
+        if (t % estimator == ESTIMATOR_COLLISION) then
+          score = p % last_wgt * material_xs % nu_fission / material_xs % total
+        else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
+          score = flux * material_xs % nu_fission
+        end if
+
       case default
         if (t % estimator == ESTIMATOR_ANALOG) then
           ! Any other score is assumed to be a MT number. Thus, we just need
@@ -732,6 +739,35 @@ contains
 
       end select
 
+      !#########################################################################
+      ! Add derivative information on score for differential tallies.
+
+      if (allocated(t % deriv) .and. t % estimator == ESTIMATOR_COLLISION) then
+        select case (score_bin)
+
+        case (SCORE_KEFF)
+          select case (t % deriv % dep_var)
+
+          case (DIFF_DENSITY)
+            score = score * t % deriv % accumulator
+
+          case (DIFF_NUCLIDE_DENSITY)
+            mat => materials(p % material)
+            if (mat % id == t % deriv % diff_material .and. &
+                 micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO) then
+              do l = 1, mat % n_nuclides
+                if (mat % nuclide(l) == t % deriv % diff_nuclide) exit
+              end do
+              score = score * (t % deriv % accumulator + ONE &
+                   / mat % atom_density(l))
+            else
+              score = score * t % deriv % accumulator
+            end if
+
+          end select
+        end select
+      end if
+
       !#########################################################################
       ! Expand score if necessary and add to tally results.
 
@@ -2210,6 +2246,80 @@ contains
 
   end subroutine score_surface_current
 
+!===============================================================================
+! SCORE_DIFF_ACCUMULATORS
+!===============================================================================
+
+  subroutine score_diff_accumulators(p, distance, collision)
+    type(particle), intent(in) :: p
+    real(8),        intent(in) :: distance
+    logical,        intent(in) :: collision
+
+    integer :: i, j
+    type(TallyObject), pointer :: t
+    type(Material), pointer :: mat
+
+    do i = 1, n_user_tallies
+      t => user_tallies(i)
+
+      if (allocated(t % deriv)) then
+        select case (t % deriv % dep_var)
+
+        case (DIFF_DENSITY)
+          if (p % material == MATERIAL_VOID) cycle
+          mat => materials(p % material)
+          if (mat % id == t % deriv % diff_material) then
+            if (collision) then
+              t % deriv % accumulator = t % deriv % accumulator &
+                   + ONE / mat % density_gpcc &
+                   - distance * material_xs % total / mat % density_gpcc
+            else
+              t % deriv % accumulator = t % deriv % accumulator &
+                   - distance * material_xs % total / mat % density_gpcc
+            end if
+          end if
+
+        case (DIFF_NUCLIDE_DENSITY)
+          if (p % material == MATERIAL_VOID) cycle
+          mat => materials(p % material)
+          if (mat % id == t % deriv % diff_material) then
+            do j = 1, mat % n_nuclides
+              if (mat % nuclide(j) == t % deriv % diff_nuclide) exit
+            end do
+            if (mat % nuclide(j) /= t % deriv % diff_nuclide) then
+              call fatal_error("Couldn't find the right nuclide.")
+            end if
+            t % deriv % accumulator = t % deriv % accumulator &
+                 - distance * micro_xs(t % deriv % diff_nuclide) % total
+            if (collision) then
+              if (p % event == EVENT_SCATTER) then
+                  t % deriv % accumulator = t % deriv % accumulator &
+                       + micro_xs(t % deriv % diff_nuclide) % elastic &
+                       / material_xs % elastic
+              else if (p % event == EVENT_ABSORB) then
+                  t % deriv % accumulator = t % deriv % accumulator &
+                       + micro_xs(t % deriv % diff_nuclide) % absorption &
+                       / material_xs % absorption
+              end if
+            end if
+          end if
+        end select
+      end if
+    end do
+  end subroutine
+
+  subroutine clear_diff_accumulators()
+    integer :: i
+    type(TallyObject), pointer :: t
+
+    do i = 1, n_user_tallies
+      t => user_tallies(i)
+      if (allocated(t % deriv)) then
+        t % deriv % accumulator = ZERO
+      end if
+    end do
+  end subroutine
+
 !===============================================================================
 ! GET_NEXT_BIN determines the next scoring bin for a particular filter variable
 !===============================================================================
diff --git a/src/tally_header.F90 b/src/tally_header.F90
index 01dcd9bdb5..cfbe0d5376 100644
--- a/src/tally_header.F90
+++ b/src/tally_header.F90
@@ -64,6 +64,18 @@ module tally_header
       procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter
   end type TallyFilter
 
+
+!===============================================================================
+! TALLYDERIVATIVE
+!===============================================================================
+
+  type TallyDerivative
+    real(8) :: accumulator
+    integer :: dep_var
+    integer :: diff_material
+    integer :: diff_nuclide
+  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
@@ -130,6 +142,9 @@ module tally_header
     integer                           :: n_triggers = 0  ! # of triggers
     type(TriggerObject),  allocatable :: triggers(:)     ! Array of triggers
 
+    ! Derivative for differentially tallies
+    type(TallyDerivative), allocatable :: deriv
+
     ! Type-Bound procedures
     contains
       procedure :: clear => tallyobject_clear ! Deallocates TallyObject
@@ -204,6 +219,9 @@ module tally_header
 
       this % n_triggers = 0
 
+      if (allocated(this % deriv)) &
+           deallocate (this % deriv)
+
     end subroutine tallyobject_clear
 
 end module tally_header
diff --git a/src/tracking.F90 b/src/tracking.F90
index 2e4503e139..9502746da7 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -13,7 +13,8 @@ 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_diff_accumulators, clear_diff_accumulators
   use track_output,    only: initialize_particle_track, write_particle_track, &
                              add_particle_track, finalize_particle_track
 
@@ -60,6 +61,9 @@ contains
       call initialize_particle_track()
     endif
 
+    ! Every particle starts with no accumulated flux derivative.
+    call clear_diff_accumulators()
+
     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
@@ -114,6 +118,9 @@ contains
              distance * material_xs % nu_fission
       end if
 
+      ! Score flux derivative accumulators for differential tallies.
+      call score_diff_accumulators(p, distance, d_boundary > d_collision)
+
       if (d_collision > d_boundary) then
         ! ====================================================================
         ! PARTICLE CROSSES SURFACE
diff --git a/tests/test_diff_density/build_xml.py b/tests/test_diff_density/build_xml.py
new file mode 100644
index 0000000000..d5e50c7e82
--- /dev/null
+++ b/tests/test_diff_density/build_xml.py
@@ -0,0 +1,150 @@
+import openmc
+
+
+def make_mats(mod_density=0.7420582):
+    water = openmc.Material(name='water')
+    water.set_density('g/cm3', mod_density)
+    water.add_nuclide(openmc.Nuclide('H-1'), 2.0)
+    water.add_nuclide(openmc.Nuclide('O-16'), 1.0)
+
+    fuel = openmc.Material(name='fuel 2.4%')
+    fuel.set_density('g/cm3', 10.29769)
+    fuel.add_nuclide(openmc.Nuclide('U-234'), 5.7987e-06)
+    fuel.add_nuclide(openmc.Nuclide('U-235'), 7.2175e-04)
+    fuel.add_nuclide(openmc.Nuclide('U-238'), 2.2253e-02)
+    fuel.add_nuclide(openmc.Nuclide('O-16'), 4.5750e-02)
+    fuel.add_nuclide(openmc.Nuclide('O-17'), 9.4222e-05)
+
+    zirc4 = openmc.Material(name='zircaloy 4')
+    zirc4.set_density('g/cm3', 6.55)
+    zirc4.add_nuclide(openmc.Nuclide('O-16'), 3.0743e-04)
+    zirc4.add_nuclide(openmc.Nuclide('O-17'), 1.1711e-07)
+    # O-18 omitted
+    zirc4.add_nuclide(openmc.Nuclide('Cr-50'), 3.2962e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Cr-52'), 6.3564e-05)
+    zirc4.add_nuclide(openmc.Nuclide('Cr-53'), 7.2076e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Cr-54'), 1.7941e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Fe-54'), 8.6699e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Fe-56'), 1.3610e-04)
+    zirc4.add_nuclide(openmc.Nuclide('Fe-57'), 3.1431e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Fe-58'), 4.1829e-07)
+    zirc4.add_nuclide(openmc.Nuclide('Zr-90'), 2.1827e-02)
+    zirc4.add_nuclide(openmc.Nuclide('Zr-91'), 4.7600e-03)
+    zirc4.add_nuclide(openmc.Nuclide('Zr-92'), 7.2758e-03)
+    zirc4.add_nuclide(openmc.Nuclide('Zr-94'), 7.3734e-03)
+    zirc4.add_nuclide(openmc.Nuclide('Zr-96'), 1.1879e-03)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-112'), 4.6735e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-114'), 3.1799e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-115'), 1.6381e-06)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-116'), 7.0055e-05)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-117'), 3.7003e-05)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-118'), 1.1669e-04)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-119'), 4.1387e-05)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-120'), 1.5697e-04)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-122'), 2.2308e-05)
+    zirc4.add_nuclide(openmc.Nuclide('Sn-124'), 2.7897e-05)
+
+    materials_file = openmc.MaterialsFile()
+    materials_file.default_xs = '71c'
+    materials_file.add_materials([water, fuel, zirc4])
+
+    return (materials_file,
+            {'mod':water,
+            'fuel':fuel,
+            'clad':zirc4})
+
+
+def make_geom(mats):
+    # Instantiate surfaces.
+    rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer')
+    rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner')
+    rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer')
+
+    x0 = openmc.XPlane(x0=-0.62992)
+    x1 = openmc.XPlane(x0=0.62992)
+    y0 = openmc.YPlane(y0=-0.62992)
+    y1 = openmc.YPlane(y0=0.62992)
+    z0 = openmc.ZPlane(z0=-1.0)
+    z1 = openmc.ZPlane(z0=1.0)
+    x0.boundary_type = 'reflective'
+    x1.boundary_type = 'reflective'
+    y0.boundary_type = 'reflective'
+    y1.boundary_type = 'reflective'
+    z0.boundary_type = 'reflective'
+    z1.boundary_type = 'reflective'
+
+    # Instantiate cells.
+    fuel_c = openmc.Cell(name='fuel')
+    fuel_c.add_surface(rfo, halfspace=-1)
+    fuel_c.add_surface(z0, halfspace=+1)
+    fuel_c.add_surface(z1, halfspace=-1)
+
+    gap = openmc.Cell(name='gap')
+    gap.add_surface(rci, halfspace=-1)
+    gap.add_surface(rfo, halfspace=+1)
+    gap.add_surface(z0, halfspace=+1)
+    gap.add_surface(z1, halfspace=-1)
+
+    clad_c = openmc.Cell(name='clad')
+    clad_c.add_surface(rco, halfspace=-1)
+    clad_c.add_surface(rci, halfspace=+1)
+    clad_c.add_surface(z0, halfspace=+1)
+    clad_c.add_surface(z1, halfspace=-1)
+
+    mod_c = openmc.Cell(name='moderator')
+    mod_c.add_surface(rco, halfspace=+1)
+    mod_c.add_surface(x0, halfspace=+1)
+    mod_c.add_surface(x1, halfspace=-1)
+    mod_c.add_surface(y0, halfspace=+1)
+    mod_c.add_surface(y1, halfspace=-1)
+    mod_c.add_surface(z0, halfspace=+1)
+    mod_c.add_surface(z1, halfspace=-1)
+
+    # Add materials to cells.
+    fuel_c.fill = mats['fuel']
+    gap.fill = 'void'
+    clad_c.fill = mats['clad']
+    mod_c.fill = mats['mod']
+
+    # Instantiate universes.
+    u0 = openmc.Universe(universe_id=0)
+    u0.add_cells([fuel_c, gap, clad_c, mod_c])
+
+    # Write the XML file.
+    geometry = openmc.Geometry()
+    geometry.root_universe = u0
+    geometry_file = openmc.GeometryFile()
+    geometry_file.geometry = geometry
+
+    return geometry_file
+
+
+def make_settings(**kwargs):
+    if 'batches' in kwargs:
+        batches = kwargs['batches']
+
+    settings_file = openmc.SettingsFile()
+    settings_file.batches = kwargs.setdefault('batches', 100)
+    settings_file.inactive = kwargs.setdefault('inactive', 10)
+    settings_file.particles = kwargs.setdefault('particles', 1000)
+    settings_file.set_source_space('box', (-0.6, -0.6, -0.9,
+                                           0.6, 0.6, 0.9))
+    settings_file.entropy_lower_left = (-0.7, -0.7, -1.1)
+    settings_file.entropy_upper_right = (0.7, 0.7, 1.1)
+    settings_file.entropy_dimension = (10, 10, 10)
+    return settings_file
+
+
+def make_inputs(mod_density=0.7420582, **kwargs):
+    (mat_file, mats) = make_mats(mod_density)
+    mat_file.export_to_xml()
+
+    geo_file = make_geom(mats)
+    geo_file.export_to_xml()
+
+    sets_file = make_settings(**kwargs)
+    sets_file.export_to_xml()
+
+
+if __name__ == '__main__':
+    make_inputs(inactive=5, batches=10, particles=100)
diff --git a/tests/test_diff_density/geometry.xml b/tests/test_diff_density/geometry.xml
new file mode 100644
index 0000000000..e7e5b097ba
--- /dev/null
+++ b/tests/test_diff_density/geometry.xml
@@ -0,0 +1,16 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+    <cell id="10000" material="10001" name="fuel" region="(-10000 10007 -10008)" universe="0" />
+    <cell id="10001" material="void" name="gap" region="(-10001 10000 10007 -10008)" universe="0" />
+    <cell id="10002" material="10002" name="clad" region="(-10002 10001 10007 -10008)" universe="0" />
+    <cell id="10003" material="10000" name="moderator" region="(10002 10003 -10004 10005 -10006 10007 -10008)" universe="0" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.39218" id="10000" name="fuel outer" type="z-cylinder" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.40005" id="10001" name="clad inner" type="z-cylinder" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.4572" id="10002" name="clad outer" type="z-cylinder" />
+    <surface boundary="reflective" coeffs="-0.62992" id="10003" type="x-plane" />
+    <surface boundary="reflective" coeffs="0.62992" id="10004" type="x-plane" />
+    <surface boundary="reflective" coeffs="-0.62992" id="10005" type="y-plane" />
+    <surface boundary="reflective" coeffs="0.62992" id="10006" type="y-plane" />
+    <surface boundary="reflective" coeffs="-1.0" id="10007" type="z-plane" />
+    <surface boundary="reflective" coeffs="1.0" id="10008" type="z-plane" />
+</geometry>
diff --git a/tests/test_diff_density/materials.xml b/tests/test_diff_density/materials.xml
new file mode 100644
index 0000000000..955ea92e7d
--- /dev/null
+++ b/tests/test_diff_density/materials.xml
@@ -0,0 +1,45 @@
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+    <default_xs>71c</default_xs>
+    <material id="10000" name="water">
+        <density units="g/cm3" value="0.7420582" />
+        <nuclide ao="2.0" name="H-1" />
+        <nuclide ao="1.0" name="O-16" />
+    </material>
+    <material id="10001" name="fuel 2.4%">
+        <density units="g/cm3" value="10.29769" />
+        <nuclide ao="5.7987e-06" name="U-234" />
+        <nuclide ao="0.00072175" name="U-235" />
+        <nuclide ao="0.022253" name="U-238" />
+        <nuclide ao="0.04575" name="O-16" />
+        <nuclide ao="9.4222e-05" name="O-17" />
+    </material>
+    <material id="10002" name="zircaloy 4">
+        <density units="g/cm3" value="6.55" />
+        <nuclide ao="0.00030743" name="O-16" />
+        <nuclide ao="1.1711e-07" name="O-17" />
+        <nuclide ao="3.2962e-06" name="Cr-50" />
+        <nuclide ao="6.3564e-05" name="Cr-52" />
+        <nuclide ao="7.2076e-06" name="Cr-53" />
+        <nuclide ao="1.7941e-06" name="Cr-54" />
+        <nuclide ao="8.6699e-06" name="Fe-54" />
+        <nuclide ao="0.0001361" name="Fe-56" />
+        <nuclide ao="3.1431e-06" name="Fe-57" />
+        <nuclide ao="4.1829e-07" name="Fe-58" />
+        <nuclide ao="0.021827" name="Zr-90" />
+        <nuclide ao="0.00476" name="Zr-91" />
+        <nuclide ao="0.0072758" name="Zr-92" />
+        <nuclide ao="0.0073734" name="Zr-94" />
+        <nuclide ao="0.0011879" name="Zr-96" />
+        <nuclide ao="4.6735e-06" name="Sn-112" />
+        <nuclide ao="3.1799e-06" name="Sn-114" />
+        <nuclide ao="1.6381e-06" name="Sn-115" />
+        <nuclide ao="7.0055e-05" name="Sn-116" />
+        <nuclide ao="3.7003e-05" name="Sn-117" />
+        <nuclide ao="0.00011669" name="Sn-118" />
+        <nuclide ao="4.1387e-05" name="Sn-119" />
+        <nuclide ao="0.00015697" name="Sn-120" />
+        <nuclide ao="2.2308e-05" name="Sn-122" />
+        <nuclide ao="2.7897e-05" name="Sn-124" />
+    </material>
+</materials>
diff --git a/tests/test_diff_density/results_true.dat b/tests/test_diff_density/results_true.dat
new file mode 100644
index 0000000000..5e515c4613
--- /dev/null
+++ b/tests/test_diff_density/results_true.dat
@@ -0,0 +1,5 @@
+k-combined:
+1.502605E+00 2.273838E-02
+tally 1:
+5.280104E+00
+1.379258E+01
diff --git a/tests/test_diff_density/settings.xml b/tests/test_diff_density/settings.xml
new file mode 100644
index 0000000000..a41ba9725b
--- /dev/null
+++ b/tests/test_diff_density/settings.xml
@@ -0,0 +1,18 @@
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+    <eigenvalue>
+        <particles>100</particles>
+        <batches>10</batches>
+        <inactive>5</inactive>
+    </eigenvalue>
+    <source>
+        <space type="box">
+            <parameters>-0.6 -0.6 -0.9 0.6 0.6 0.9</parameters>
+        </space>
+    </source>
+    <entropy>
+        <dimension>10 10 10</dimension>
+        <lower_left>-0.7 -0.7 -1.1</lower_left>
+        <upper_right>0.7 0.7 1.1</upper_right>
+    </entropy>
+</settings>
diff --git a/tests/test_diff_density/tallies.xml b/tests/test_diff_density/tallies.xml
new file mode 100644
index 0000000000..214f33ab37
--- /dev/null
+++ b/tests/test_diff_density/tallies.xml
@@ -0,0 +1,7 @@
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <tally id="1">
+    <scores>keff</scores>
+    <derivative variable="density" material="10000"/>
+  </tally>
+</tallies>
diff --git a/tests/test_diff_density/test_diff_density.py b/tests/test_diff_density/test_diff_density.py
new file mode 100755
index 0000000000..ed6addec45
--- /dev/null
+++ b/tests/test_diff_density/test_diff_density.py
@@ -0,0 +1,11 @@
+#!/usr/bin/env python
+
+import os
+import sys
+sys.path.insert(0, os.pardir)
+from testing_harness import TestHarness
+
+
+if __name__ == '__main__':
+    harness = TestHarness('statepoint.10.*', True)
+    harness.main()
diff --git a/tests/test_diff_nuclide_density/build_xml.py b/tests/test_diff_nuclide_density/build_xml.py
new file mode 100644
index 0000000000..f77a337cb1
--- /dev/null
+++ b/tests/test_diff_nuclide_density/build_xml.py
@@ -0,0 +1,164 @@
+import openmc
+
+
+def make_mats(**kwargs):
+    water_dict = {
+         'B-10': 8.0042e-06,
+         'B-11': 3.2218e-05,
+         'H-1': 4.9457e-02,
+         'H-2': 7.4196e-06,
+         'O-16': 2.4672e-02,
+         'O-17': 9.3982e-06 + 5.0701e-05}
+    water = openmc.Material(name='water')
+    water.set_density('sum')
+    for nuclide in water_dict:
+        water.add_nuclide(nuclide, water_dict[nuclide])
+
+    U235_dens = kwargs.setdefault('U235_dens', 5.5814e-04)
+    fuel_dict = {
+         'B-10': 8.0042e-06,
+         'U-234': 4.4842e-06,
+         'U-235': U235_dens,
+         'U-238': 2.2407e-02,
+         'O-16': 4.5828e-02,
+         'O-17': 1.7457e-05 + 9.4176e-05}
+    fuel = openmc.Material(name='fuel 2.4%')
+    fuel.set_density('sum')
+    for nuclide in fuel_dict:
+        fuel.add_nuclide(nuclide, fuel_dict[nuclide])
+
+    zirc4_dict = {
+         'O-16': 3.0743e-04,
+         'O-17': 1.1711e-07 + 6.3176e-07,
+         'Cr-50': 3.2962e-06,
+         'Cr-52': 6.3564e-05,
+         'Cr-53': 7.2076e-06,
+         'Cr-54': 1.7941e-06,
+         'Fe-54': 8.6699e-06,
+         'Fe-56': 1.3610e-04,
+         'Fe-57': 3.1431e-06,
+         'Fe-58': 4.1829e-07,
+         'Zr-90': 2.1827e-02,
+         'Zr-91': 4.7600e-03,
+         'Zr-92': 7.2758e-03,
+         'Zr-94': 7.3734e-03,
+         'Zr-96': 1.1879e-03,
+         'Sn-112': 4.6735e-06,
+         'Sn-114': 3.1799e-06,
+         'Sn-115': 1.6381e-06,
+         'Sn-116': 7.0055e-05,
+         'Sn-117': 3.7003e-05,
+         'Sn-118': 1.1669e-04,
+         'Sn-119': 4.1387e-05,
+         'Sn-120': 1.5697e-04,
+         'Sn-122': 2.2308e-05,
+         'Sn-124': 2.7897e-05}
+    zirc4 = openmc.Material(name='zircaloy 4')
+    zirc4.set_density('sum')
+    for nuclide in zirc4_dict:
+        zirc4.add_nuclide(nuclide, zirc4_dict[nuclide])
+
+    materials_file = openmc.MaterialsFile()
+    materials_file.default_xs = '71c'
+    materials_file.add_materials([water, fuel, zirc4])
+
+    return (materials_file,
+            {'mod':water,
+            'fuel':fuel,
+            'clad':zirc4})
+
+
+def make_geom(mats):
+    # Instantiate surfaces.
+    rfo = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218, name='fuel outer')
+    rci = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.40005, name='clad inner')
+    rco = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720, name='clad outer')
+
+    x0 = openmc.XPlane(x0=-0.62992)
+    x1 = openmc.XPlane(x0=0.62992)
+    y0 = openmc.YPlane(y0=-0.62992)
+    y1 = openmc.YPlane(y0=0.62992)
+    z0 = openmc.ZPlane(z0=-1.0)
+    z1 = openmc.ZPlane(z0=1.0)
+    x0.boundary_type = 'reflective'
+    x1.boundary_type = 'reflective'
+    y0.boundary_type = 'reflective'
+    y1.boundary_type = 'reflective'
+    z0.boundary_type = 'reflective'
+    z1.boundary_type = 'reflective'
+
+    # Instantiate cells.
+    fuel_c = openmc.Cell(name='fuel')
+    fuel_c.add_surface(rfo, halfspace=-1)
+    fuel_c.add_surface(z0, halfspace=+1)
+    fuel_c.add_surface(z1, halfspace=-1)
+
+    gap = openmc.Cell(name='gap')
+    gap.add_surface(rci, halfspace=-1)
+    gap.add_surface(rfo, halfspace=+1)
+    gap.add_surface(z0, halfspace=+1)
+    gap.add_surface(z1, halfspace=-1)
+
+    clad_c = openmc.Cell(name='clad')
+    clad_c.add_surface(rco, halfspace=-1)
+    clad_c.add_surface(rci, halfspace=+1)
+    clad_c.add_surface(z0, halfspace=+1)
+    clad_c.add_surface(z1, halfspace=-1)
+
+    mod_c = openmc.Cell(name='moderator')
+    mod_c.add_surface(rco, halfspace=+1)
+    mod_c.add_surface(x0, halfspace=+1)
+    mod_c.add_surface(x1, halfspace=-1)
+    mod_c.add_surface(y0, halfspace=+1)
+    mod_c.add_surface(y1, halfspace=-1)
+    mod_c.add_surface(z0, halfspace=+1)
+    mod_c.add_surface(z1, halfspace=-1)
+
+    # Add materials to cells.
+    fuel_c.fill = mats['fuel']
+    gap.fill = 'void'
+    clad_c.fill = mats['clad']
+    mod_c.fill = mats['mod']
+
+    # Instantiate universes.
+    u0 = openmc.Universe(universe_id=0)
+    u0.add_cells([fuel_c, gap, clad_c, mod_c])
+
+    # Write the XML file.
+    geometry = openmc.Geometry()
+    geometry.root_universe = u0
+    geometry_file = openmc.GeometryFile()
+    geometry_file.geometry = geometry
+
+    return geometry_file
+
+
+def make_settings(**kwargs):
+    if 'batches' in kwargs:
+        batches = kwargs['batches']
+
+    settings_file = openmc.SettingsFile()
+    settings_file.batches = kwargs.setdefault('batches', 1610)
+    settings_file.inactive = kwargs.setdefault('inactive', 10)
+    settings_file.particles = kwargs.setdefault('particles', 1000)
+    settings_file.set_source_space('box', (-0.6, -0.6, -0.9,
+                                           0.6, 0.6, 0.9))
+    settings_file.entropy_lower_left = (-0.7, -0.7, -1.1)
+    settings_file.entropy_upper_right = (0.7, 0.7, 1.1)
+    settings_file.entropy_dimension = (10, 10, 10)
+    return settings_file
+
+
+def make_inputs(**kwargs):
+    (mat_file, mats) = make_mats(**kwargs)
+    mat_file.export_to_xml()
+
+    geo_file = make_geom(mats)
+    geo_file.export_to_xml()
+
+    sets_file = make_settings(**kwargs)
+    sets_file.export_to_xml()
+
+
+if __name__ == '__main__':
+    make_inputs(inactive=5, batches=10, particles=100)
diff --git a/tests/test_diff_nuclide_density/geometry.xml b/tests/test_diff_nuclide_density/geometry.xml
new file mode 100644
index 0000000000..e7e5b097ba
--- /dev/null
+++ b/tests/test_diff_nuclide_density/geometry.xml
@@ -0,0 +1,16 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+    <cell id="10000" material="10001" name="fuel" region="(-10000 10007 -10008)" universe="0" />
+    <cell id="10001" material="void" name="gap" region="(-10001 10000 10007 -10008)" universe="0" />
+    <cell id="10002" material="10002" name="clad" region="(-10002 10001 10007 -10008)" universe="0" />
+    <cell id="10003" material="10000" name="moderator" region="(10002 10003 -10004 10005 -10006 10007 -10008)" universe="0" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.39218" id="10000" name="fuel outer" type="z-cylinder" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.40005" id="10001" name="clad inner" type="z-cylinder" />
+    <surface boundary="transmission" coeffs="0.0 0.0 0.4572" id="10002" name="clad outer" type="z-cylinder" />
+    <surface boundary="reflective" coeffs="-0.62992" id="10003" type="x-plane" />
+    <surface boundary="reflective" coeffs="0.62992" id="10004" type="x-plane" />
+    <surface boundary="reflective" coeffs="-0.62992" id="10005" type="y-plane" />
+    <surface boundary="reflective" coeffs="0.62992" id="10006" type="y-plane" />
+    <surface boundary="reflective" coeffs="-1.0" id="10007" type="z-plane" />
+    <surface boundary="reflective" coeffs="1.0" id="10008" type="z-plane" />
+</geometry>
diff --git a/tests/test_diff_nuclide_density/materials.xml b/tests/test_diff_nuclide_density/materials.xml
new file mode 100644
index 0000000000..02b55f4523
--- /dev/null
+++ b/tests/test_diff_nuclide_density/materials.xml
@@ -0,0 +1,50 @@
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+    <default_xs>71c</default_xs>
+    <material id="10000" name="water">
+        <density units="sum" />
+        <nuclide ao="6.00992e-05" name="O-17" />
+        <nuclide ao="0.024672" name="O-16" />
+        <nuclide ao="8.0042e-06" name="B-10" />
+        <nuclide ao="3.2218e-05" name="B-11" />
+        <nuclide ao="0.049457" name="H-1" />
+        <nuclide ao="7.4196e-06" name="H-2" />
+    </material>
+    <material id="10001" name="fuel 2.4%">
+        <density units="sum" />
+        <nuclide ao="0.022407" name="U-238" />
+        <nuclide ao="0.000111633" name="O-17" />
+        <nuclide ao="8.0042e-06" name="B-10" />
+        <nuclide ao="0.045828" name="O-16" />
+        <nuclide ao="4.4842e-06" name="U-234" />
+        <nuclide ao="0.00055814" name="U-235" />
+    </material>
+    <material id="10002" name="zircaloy 4">
+        <density units="sum" />
+        <nuclide ao="0.00011669" name="Sn-118" />
+        <nuclide ao="4.1387e-05" name="Sn-119" />
+        <nuclide ao="4.6735e-06" name="Sn-112" />
+        <nuclide ao="3.1799e-06" name="Sn-114" />
+        <nuclide ao="1.6381e-06" name="Sn-115" />
+        <nuclide ao="7.0055e-05" name="Sn-116" />
+        <nuclide ao="3.7003e-05" name="Sn-117" />
+        <nuclide ao="0.0001361" name="Fe-56" />
+        <nuclide ao="3.1431e-06" name="Fe-57" />
+        <nuclide ao="8.6699e-06" name="Fe-54" />
+        <nuclide ao="4.1829e-07" name="Fe-58" />
+        <nuclide ao="0.0072758" name="Zr-92" />
+        <nuclide ao="2.7897e-05" name="Sn-124" />
+        <nuclide ao="0.00015697" name="Sn-120" />
+        <nuclide ao="2.2308e-05" name="Sn-122" />
+        <nuclide ao="1.7941e-06" name="Cr-54" />
+        <nuclide ao="6.3564e-05" name="Cr-52" />
+        <nuclide ao="7.2076e-06" name="Cr-53" />
+        <nuclide ao="3.2962e-06" name="Cr-50" />
+        <nuclide ao="7.4887e-07" name="O-17" />
+        <nuclide ao="0.00030743" name="O-16" />
+        <nuclide ao="0.00476" name="Zr-91" />
+        <nuclide ao="0.021827" name="Zr-90" />
+        <nuclide ao="0.0011879" name="Zr-96" />
+        <nuclide ao="0.0073734" name="Zr-94" />
+    </material>
+</materials>
diff --git a/tests/test_diff_nuclide_density/results_true.dat b/tests/test_diff_nuclide_density/results_true.dat
new file mode 100644
index 0000000000..68b611ae15
--- /dev/null
+++ b/tests/test_diff_nuclide_density/results_true.dat
@@ -0,0 +1,5 @@
+k-combined:
+1.055198E+00 5.785829E-02
+tally 1:
+3.081784E+03
+1.905482E+06
diff --git a/tests/test_diff_nuclide_density/settings.xml b/tests/test_diff_nuclide_density/settings.xml
new file mode 100644
index 0000000000..a41ba9725b
--- /dev/null
+++ b/tests/test_diff_nuclide_density/settings.xml
@@ -0,0 +1,18 @@
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+    <eigenvalue>
+        <particles>100</particles>
+        <batches>10</batches>
+        <inactive>5</inactive>
+    </eigenvalue>
+    <source>
+        <space type="box">
+            <parameters>-0.6 -0.6 -0.9 0.6 0.6 0.9</parameters>
+        </space>
+    </source>
+    <entropy>
+        <dimension>10 10 10</dimension>
+        <lower_left>-0.7 -0.7 -1.1</lower_left>
+        <upper_right>0.7 0.7 1.1</upper_right>
+    </entropy>
+</settings>
diff --git a/tests/test_diff_nuclide_density/tallies.xml b/tests/test_diff_nuclide_density/tallies.xml
new file mode 100644
index 0000000000..321b2c9b0e
--- /dev/null
+++ b/tests/test_diff_nuclide_density/tallies.xml
@@ -0,0 +1,10 @@
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <tally id="1">
+    <scores>keff</scores>
+    <!--
+    <derivative variable="density" material="10000"/>
+    -->
+    <derivative variable="nuclide_density" material="10001" nuclide="U-235.71c"/>
+  </tally>
+</tallies>
diff --git a/tests/test_diff_nuclide_density/test_diff_nuclide_density.py b/tests/test_diff_nuclide_density/test_diff_nuclide_density.py
new file mode 100755
index 0000000000..ed6addec45
--- /dev/null
+++ b/tests/test_diff_nuclide_density/test_diff_nuclide_density.py
@@ -0,0 +1,11 @@
+#!/usr/bin/env python
+
+import os
+import sys
+sys.path.insert(0, os.pardir)
+from testing_harness import TestHarness
+
+
+if __name__ == '__main__':
+    harness = TestHarness('statepoint.10.*', True)
+    harness.main()