diff --git a/include/userobjects/DPAUserObject.h b/include/userobjects/DPAUserObject.h
new file mode 100644
index 00000000..1f69acd2
--- /dev/null
+++ b/include/userobjects/DPAUserObject.h
@@ -0,0 +1,84 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+#ifdef GSL_ENABLED
+
+#pragma once
+
+#include "ParkinCoulterBase.h"
+
+class DPAUserObject;
+
+template <>
+InputParameters validParams<DPAUserObject>();
+
+class DPAUserObject : public ParkinCoulterBase
+{
+public:
+  DPAUserObject(const InputParameters & parameters);
+  void finalize() override;
+  void execute() override;
+  void initialSetup() override;
+
+protected:
+  virtual void initDamageFunctions() override;
+  virtual std::vector<unsigned int> atomicNumbers() const override;
+  virtual std::vector<Real> massNumbers() const override;
+  virtual std::vector<Real> numberFractions() const override;
+  virtual Real maxEnergy() const override;
+
+  /// determines if an update of the displacement function is required
+  bool recomputeDisplacementFunction() const;
+
+  /// a helper function that sets _ns and checks consistency of A, Z, N
+  void initAZNHelper();
+
+  /// a helper function that computes the neutron damage efficiency
+  Real
+  neutronDamageEfficiency(unsigned int i, unsigned int j, unsigned int g, unsigned int x) const;
+
+  /// tolerance for recomputing the displacement function
+  Real _tol = 1e-10;
+
+  /// the computed dose rates
+  Real _dpa = 0;
+
+  /// is damage accumulated during a transient or computed for steady state
+  bool _is_transient_irradiation;
+
+  /// irradiation_time used when dpa is estimated from steady-state calculations
+  Real _irradiation_time;
+
+  /// the neutron reaction types considered for computing dpa
+  MultiMooseEnum _neutron_reaction_types;
+
+  /// convenience copy of this parameter
+  std::vector<std::vector<std::string>> _xs_names;
+
+  ///@{ data used for computing dpa value
+  const VectorPostprocessorValue & _atomic_numbers;
+  const VectorPostprocessorValue & _mass_numbers;
+  const VectorPostprocessorValue & _number_densities;
+  const VectorPostprocessorValue & _energy_group_boundaries;
+  const VectorPostprocessorValue & _scalar_flux;
+  std::vector<std::vector<const VectorPostprocessorValue *>> _cross_sections;
+  ///@}
+
+  ///@{ the "old" versions of the data; used for determining if disp function update is required
+  std::vector<Real> _atomic_numbers_old;
+  std::vector<Real> _mass_numbers_old;
+  std::vector<Real> _number_densities_old;
+  ///@}
+
+  /// number of neutron energy groups
+  unsigned int _ng;
+
+  /// number of reaction types creating radiation damage
+  unsigned int _nr;
+};
+
+#endif // GSL_ENABLED
diff --git a/src/userobjects/DPAUserObject.C b/src/userobjects/DPAUserObject.C
new file mode 100644
index 00000000..64b76715
--- /dev/null
+++ b/src/userobjects/DPAUserObject.C
@@ -0,0 +1,318 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+#ifdef GSL_ENABLED
+
+#include "DPAUserObject.h"
+#include "PolyatomicRecoil.h"
+#include "PolyatomicDisplacementFunction.h"
+#include "PolyatomicDamageEnergyFunction.h"
+#include "PolyatomicDisplacementDerivativeFunction.h"
+#include "MooseMesh.h"
+
+// various other includes
+#include <mytrim/element.h>
+#include <limits>
+
+registerMooseObject("MagpieApp", DPAUserObject);
+
+template <>
+InputParameters
+validParams<DPAUserObject>()
+{
+  InputParameters params = validParams<ParkinCoulterBase>();
+  params.addParam<Real>("irradiation_time", "Irradiation time used ");
+  MultiMooseEnum damage_reaction_types("elastic");
+  params.addParam<MultiMooseEnum>("damage_reaction_types",
+                                  damage_reaction_types,
+                                  "The neutron reaction causing radiation damage");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "Z", "The name of the VectorPostprocessor containing the atomic numbers");
+  params.addRequiredParam<std::string>("Z_column",
+                                       "VectorPostprocessor column containing the atomic numbers");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "A", "The name of the VectorPostprocessor containing the mass numbers");
+  params.addRequiredParam<std::string>("A_column",
+                                       "VectorPostprocessor column containing the mass numbers");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "number_densities", "The name of the VectorPostprocessor containing the number densities");
+  params.addRequiredParam<std::string>(
+      "number_densities_column", "VectorPostprocessor column containing the number densities");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "energy_group_boundaries",
+      "The name of the VectorPostprocessor containing the neutron flux energy group boundaries in "
+      "units of eV starting with the highest energy group");
+  params.addRequiredParam<std::string>(
+      "energy_group_boundaries_column",
+      "VectorPostprocessor column containing the neutron flux energy group boundaries");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "scalar_flux", "The name of the VectorPostprocessor containing the neutron scalar flux");
+  params.addRequiredParam<std::string>(
+      "scalar_flux_column", "VectorPostprocessor column containing the neutron scalar flux");
+  params.addRequiredParam<VectorPostprocessorName>(
+      "cross_section", "The name of the VectorPostprocessor containing the cross sections");
+  params.addRequiredParam<std::vector<std::vector<std::string>>>(
+      "cross_section_column",
+      "VectorPostprocessor columns containing cross sections. Parameter is a vector of vector of "
+      "strings. The first index runs over all reaction types, the second index runs over all "
+      "isotopes. Each vector must be number of groups long");
+  params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_END;
+  params.suppressParameter<ExecFlagEnum>("execute_on");
+  params.addClassDescription(
+      "Computes the dose in dpa from composition, cross section, damage type, and neutron flux.");
+  return params;
+}
+
+DPAUserObject::DPAUserObject(const InputParameters & parameters)
+  : ParkinCoulterBase(parameters),
+    _is_transient_irradiation(!isParamValid("irradiation_time")),
+    _irradiation_time(_is_transient_irradiation ? 0 : getParam<Real>("irradiation_time")),
+    _neutron_reaction_types(getParam<MultiMooseEnum>("damage_reaction_types")),
+    _xs_names(getParam<std::vector<std::vector<std::string>>>("cross_section_column")),
+    _atomic_numbers(getVectorPostprocessorValue("Z", getParam<std::string>("Z_column"))),
+    _mass_numbers(getVectorPostprocessorValue("A", getParam<std::string>("A_column"))),
+    _number_densities(getVectorPostprocessorValue(
+        "number_densities", getParam<std::string>("number_densities_column"))),
+    _energy_group_boundaries(getVectorPostprocessorValue(
+        "energy_group_boundaries", getParam<std::string>("energy_group_boundaries_column"))),
+    _scalar_flux(
+        getVectorPostprocessorValue("scalar_flux", getParam<std::string>("scalar_flux_column")))
+{
+  _cross_sections.resize(_xs_names.size());
+  for (unsigned int i = 0; i < _xs_names.size(); ++i)
+  {
+    _cross_sections[i].resize(_xs_names[i].size());
+    for (unsigned int j = 0; j < _xs_names[i].size(); ++j)
+      _cross_sections[i][j] = &getVectorPostprocessorValue("cross_section", _xs_names[i][j]);
+  }
+}
+
+void
+DPAUserObject::initialSetup()
+{
+  // checks on scalar flux & energy groups
+  _ng = _scalar_flux.size();
+  if (_ng == 0)
+    paramError("scalar_flux_column",
+               "Size of VPP is zero. scalar_flux_column",
+               getParam<std::string>("scalar_flux_column"),
+               "may not exist in VPP.");
+  if (_energy_group_boundaries.size() != _ng + 1)
+    paramError("energy_group_boundaries",
+               "Parameter requires number of neutron energy groups + 1 entries. G = ",
+               _ng);
+  for (unsigned int j = 1; j < _energy_group_boundaries.size(); ++j)
+    if (_energy_group_boundaries[j - 1] <= _energy_group_boundaries[j])
+      paramError("energy_group_boundaries", "Entries must be strictly decreasing");
+
+  // check on A, Z, number densities
+  initAZNHelper();
+
+  // checks on cross sections
+  _nr = _neutron_reaction_types.size();
+  if (_cross_sections.size() != _nr)
+    paramError("cross_section_column",
+               "Leading dimension is ",
+               _cross_sections.size(),
+               " but should be equal to number of damage reaction types ",
+               _nr);
+  for (unsigned int i = 0; i < _nr; ++i)
+  {
+    if (_cross_sections[i].size() != _atomic_numbers.size())
+      paramError("cross_section_column",
+                 "Second dimension for index ",
+                 i,
+                 " does not match number of isotope ",
+                 _atomic_numbers.size());
+    for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
+      if (_cross_sections[i][j]->size() != _ng)
+        paramError("cross_section",
+                   "Column named ",
+                   _xs_names[i][j],
+                   " has ",
+                   _cross_sections[i][j]->size(),
+                   " entries when number of groups is ",
+                   _ng);
+  }
+
+  if (_nr == 0)
+    paramError("damage_reaction_types", "At least one damage mechanism must be provided");
+}
+
+void
+DPAUserObject::initAZNHelper()
+{
+  unsigned int ns = atomicNumbers().size();
+  if (ns == 0)
+    paramError("Z_column",
+               "Size of VPP is zero. Z_column",
+               getParam<std::string>("Z_column"),
+               "may not exist in VPP.");
+
+  if (_mass_numbers.size() != ns)
+    paramError("A", "Size of VPP is not equal to size of atomic number array.");
+
+  if (_number_densities.size() != ns)
+    paramError("number_densities", "Size of VPP is not equal to size of atomic number array.");
+}
+
+std::vector<unsigned int>
+DPAUserObject::atomicNumbers() const
+{
+  // in this kind we need to work a bit since
+  // VPPs are Reals and we need to return unsigned int
+  std::vector<unsigned int> Zs(_atomic_numbers.size());
+  for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
+  {
+    unsigned int i = static_cast<unsigned int>(_atomic_numbers[j]);
+    if (std::abs(_atomic_numbers[j] - i) > 1e-12)
+      paramError("Z", "Entry ", j, ":", _atomic_numbers[j], " is not a non-negative integer.");
+    Zs[j] = i;
+  }
+  return Zs;
+}
+
+std::vector<Real>
+DPAUserObject::massNumbers() const
+{
+  return _mass_numbers;
+}
+
+std::vector<Real>
+DPAUserObject::numberFractions() const
+{
+  // need to normalize here
+  Real sum = 0;
+  for (unsigned int j = 0; j < _number_densities.size(); ++j)
+    sum += _number_densities[j];
+  std::vector<Real> nf(_number_densities.size());
+  for (unsigned int j = 0; j < _number_densities.size(); ++j)
+    nf[j] = _number_densities[j] / sum;
+  return nf;
+}
+
+Real
+DPAUserObject::maxEnergy() const
+{
+  Real max_neutron_energy = _energy_group_boundaries[0];
+
+  // find maximum gamma
+  Real min_mass = std::numeric_limits<Real>::max();
+  for (unsigned int j = 0; j < _mass_numbers.size(); ++j)
+    if (min_mass > _mass_numbers[j])
+      min_mass = _mass_numbers[j];
+  Real gamma = 4 * min_mass / (min_mass + 1) / (min_mass + 1);
+  return gamma * max_neutron_energy;
+}
+
+void
+DPAUserObject::initDamageFunctions()
+{
+  _padf = libmesh_make_unique<PolyatomicDisplacementFunction>(polyMat(), NET, _Ecap);
+}
+
+void
+DPAUserObject::execute()
+{
+  if (recomputeDisplacementFunction())
+  {
+    initAZNHelper();
+    computeDamageFunctions();
+    _padf->computeDisplacementFunctionIntegral();
+
+    // copy over A, Z, N to A, Z, N_old
+    _atomic_numbers_old.resize(_atomic_numbers.size());
+    _mass_numbers_old.resize(_atomic_numbers.size());
+    _number_densities_old.resize(_atomic_numbers.size());
+    for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
+    {
+      _atomic_numbers_old[j] = _atomic_numbers[j];
+      _mass_numbers_old[j] = _mass_numbers[j];
+      _number_densities_old[j] = _number_densities[j];
+    }
+  }
+
+  // damage function have been computed
+  // so now we can start computing dpa
+  if (!_is_transient_irradiation)
+    _dpa = 0;
+
+  Real del_t = _is_transient_irradiation ? _dt : _irradiation_time;
+
+  for (unsigned int x = 0; x < _nr; ++x)
+  {
+    for (unsigned int g = 0; g < _ng; ++g)
+    {
+      Real del_E = _energy_group_boundaries[g] - _energy_group_boundaries[g + 1];
+
+      for (unsigned int i = 0; i < _atomic_numbers.size(); ++i)
+        for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
+          _dpa += del_t / del_E * _scalar_flux[g] * _number_densities[i] *
+                  (*_cross_sections[x][i])[g] * neutronDamageEfficiency(i, j, g, x);
+    }
+  }
+}
+
+bool
+DPAUserObject::recomputeDisplacementFunction() const
+{
+  // the size could have changed
+  if (_atomic_numbers.size() != _atomic_numbers_old.size() ||
+      _mass_numbers.size() != _mass_numbers_old.size() ||
+      _number_densities.size() != _number_densities_old.size())
+    return true;
+
+  // the size is still the same
+  for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
+    if (!MooseUtils::absoluteFuzzyEqual(_atomic_numbers[j], _atomic_numbers_old[j], _tol) ||
+        !MooseUtils::absoluteFuzzyEqual(_mass_numbers[j], _mass_numbers_old[j], _tol) ||
+        !MooseUtils::absoluteFuzzyEqual(_number_densities[j], _number_densities_old[j], _tol))
+      return true;
+
+  return false;
+}
+
+Real
+DPAUserObject::neutronDamageEfficiency(unsigned int i,
+                                       unsigned int j,
+                                       unsigned int g,
+                                       unsigned int x) const
+{
+  std::vector<Real> points;
+  std::vector<Real> weights;
+  PolyatomicDisplacementFunction::gslQuadRule(100, points, weights);
+
+  if (_neutron_reaction_types[x] == "elastic")
+  {
+    Real gamma = 4 * _mass_numbers[i] / (_mass_numbers[i] + 1) / (_mass_numbers[i] + 1);
+    Real from_E = _energy_group_boundaries[g + 1];
+    Real to_E = _energy_group_boundaries[g];
+    Real delta_E = to_E - from_E;
+
+    // integral over E
+    Real integral = 0;
+    for (unsigned int qp = 0; qp < points.size(); ++qp)
+    {
+      Real energy = 0.5 * (points[qp] + 1) * delta_E + from_E;
+      Real w = 0.5 * weights[qp] * delta_E;
+      integral +=
+          w * _padf->linearInterpolationIntegralDamageFunction(gamma * energy, i, j) / energy;
+    }
+
+    return integral / gamma;
+  }
+
+  mooseError("Neutron reaction type not recognized. Should never get here.");
+  return 0;
+}
+
+void
+DPAUserObject::finalize()
+{
+}
+
+#endif