diff --git a/doc/content/source/auxkernels/MDNParticleAux.md b/doc/content/source/auxkernels/MDNParticleAux.md
new file mode 100644
index 00000000..fdba1adc
--- /dev/null
+++ b/doc/content/source/auxkernels/MDNParticleAux.md
@@ -0,0 +1,14 @@
+# MDNParticleAux
+
+!syntax description /AuxKernels/MDNParticleAux
+
+Injects the number of MD particles per element computed by
+MD runner object into an auxiliary variable.
+
+!syntax parameters /AuxKernels/MDNParticleAux
+
+!syntax inputs /AuxKernels/MDNParticleAux
+
+!syntax children /AuxKernels/MDNParticleAux
+
+!bibtex bibliography
diff --git a/doc/content/source/userobjects/LAMMPSFileRunner.md b/doc/content/source/userobjects/LAMMPSFileRunner.md
new file mode 100644
index 00000000..3a758b52
--- /dev/null
+++ b/doc/content/source/userobjects/LAMMPSFileRunner.md
@@ -0,0 +1,21 @@
+# LAMMPSFileRunner
+
+`LAMMPSFileRunner` reads in LAMMPS dumps file and maps them on an FEM mesh.
+The `time_sequence` parameter determines if a sequence of LAMMPS files is read
+or a single file is loaded and used throughout. If a sequence of files is used,
+the `lammps_file` parameter should contain the file base, e.g. `path/to/file/filebase`
+if the files are called `path/to/file/filebase.<time>.xyz`.
+
+Time conversion from FEM time to LAMMPS time is facilitated by the `time_conversion`
+parameter that accepts a MOOSE function object. Denoting FEM time by `t` and LAMMPS time
+by `T`, `time_conversion` is a function `T = F(t)`.
+
+!syntax description /UserObjects/LAMMPSFileRunner
+
+!syntax parameters /UserObjects/LAMMPSFileRunner
+
+!syntax inputs /UserObjects/LAMMPSFileRunner
+
+!syntax children /UserObjects/LAMMPSFileRunner
+
+!bibtex bibliography
diff --git a/include/auxkernels/MDNParticleAux.h b/include/auxkernels/MDNParticleAux.h
new file mode 100644
index 00000000..8e7f6b23
--- /dev/null
+++ b/include/auxkernels/MDNParticleAux.h
@@ -0,0 +1,33 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#ifndef MDNPARTICLEAUX_H
+#define MDNPARTICLEAUX_H
+
+#include "AuxKernel.h"
+
+// forward declarations
+class MDRunBase;
+class MDNParticleAux;
+
+template <>
+InputParameters validParams<MDNParticleAux>();
+
+class MDNParticleAux : public AuxKernel
+{
+public:
+  MDNParticleAux(const InputParameters & params);
+  virtual ~MDNParticleAux() {}
+
+  virtual Real computeValue();
+
+protected:
+  const MDRunBase & _md_uo;
+};
+
+#endif // MDNPARTICLEAUX_H
diff --git a/include/userobjects/LAMMPSFileRunner.h b/include/userobjects/LAMMPSFileRunner.h
new file mode 100644
index 00000000..3756fa46
--- /dev/null
+++ b/include/userobjects/LAMMPSFileRunner.h
@@ -0,0 +1,67 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#ifndef LAMMPSFILERUNNER_H
+#define LAMMPSFILERUNNER_H
+
+#include "MDRunBase.h"
+
+class LAMMPSFileRunner;
+
+template <>
+InputParameters validParams<LAMMPSFileRunner>();
+
+/**
+ * Reads lammps dump files to emulate MD simulation
+ */
+class LAMMPSFileRunner : public MDRunBase
+{
+public:
+  LAMMPSFileRunner(const InputParameters & parameters);
+
+  virtual void initialize() override {}
+  virtual void execute() override {}
+  virtual void finalize() override {}
+  virtual void initialSetup() override;
+
+  virtual void updateParticleList() override;
+
+protected:
+  /// reads a LAMMPS file
+  void readLAMMPSFile(FileName filename);
+
+  /// reads two LAMMPS files and interpolates times
+  void readLAMMPSFileHistory(std::pair<FileName, FileName> filenames,
+                             std::pair<Real, Real> timestamps,
+                             Real current_time);
+
+  /// helper function that finds two files, one right before and one right after md_time
+  void findBracketingLAMMPSFiles(Real md_time,
+                                 std::pair<std::string, std::string> & filenames,
+                                 std::pair<Real, Real> & timestamps);
+
+  /// whether a sequence of input files or a single input file is read
+  bool _time_sequence;
+
+  /// name of LAMMPS file or file base if _time_sequence == true
+  FileName _lammps_file;
+
+  /// column of , y, z coordinate in LAMMPS files
+  std::vector<unsigned int> _pos_columns;
+
+  /// get velocity from file or compute it from positions
+  bool _velocity_from_file;
+
+  /// column of , y, z coordinate in LAMMPS files
+  std::vector<unsigned int> _vel_columns;
+
+  /// Conversion from FEM time to MD time_stamp
+  Function * _time_conversion;
+};
+
+#endif // LAMMPSFILERUNNER_H
diff --git a/include/userobjects/MDRunBase.h b/include/userobjects/MDRunBase.h
new file mode 100644
index 00000000..95acc97b
--- /dev/null
+++ b/include/userobjects/MDRunBase.h
@@ -0,0 +1,93 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#ifndef MDRUNBASE_H
+#define MDRUNBASE_H
+
+#include "GeneralUserObject.h"
+#include "KDTree.h"
+#include "libmesh/bounding_box.h"
+
+class MDRunBase;
+class MooseMesh;
+
+template <>
+InputParameters validParams<MDRunBase>();
+
+/**
+ * Base class for molecular dynamics runs in Magpie
+ */
+class MDRunBase : public GeneralUserObject
+{
+public:
+  MDRunBase(const InputParameters & parameters);
+
+  void initialSetup() override;
+  void timestepSetup() override;
+
+  class MDParticles
+  {
+  public:
+    /// particle's position and history
+    std::vector<Point> pos;
+
+    /// particle's velocity
+    std::vector<Point> vel;
+
+    /// the id of particle in the MD calculation
+    std::vector<unsigned int> id;
+
+    /// the mesh element the particles are in
+    std::vector<dof_id_type> elem_id;
+  };
+
+  /// access to the MDParticles
+  const MDParticles & particles() const { return _md_particles; }
+
+  /// access to the element to particle map
+  const std::vector<unsigned int> elemParticles(dof_id_type elem_id) const;
+
+protected:
+  /// call back function to update the particle list
+  virtual void updateParticleList() = 0;
+
+  /// updates the KDTree object
+  void updateKDTree();
+
+  /// map MDParticles to elements
+  void mapMDParticles();
+
+  /// The Mesh we're using
+  MooseMesh & _mesh;
+
+  /// dimension of the mesh
+  const unsigned int _dim;
+
+  /// dimension of the mesh
+  const unsigned int _nproc;
+
+  /// stores bounding boxes of all other processors
+  std::vector<BoundingBox> _bbox;
+
+  /// total number of particles
+  unsigned int _n_particles;
+
+  /// number of local particles
+  unsigned int _n_local_particles;
+
+  /// stores the location of
+  MDParticles _md_particles;
+
+  /// a map from elem pointer to particles in this element
+  std::map<dof_id_type, std::vector<unsigned int>> _elem_particles;
+
+  /// A KDTree object to handle md_particles
+  std::unique_ptr<KDTree> _kd_tree;
+};
+
+#endif // MDRUNBASE_H
diff --git a/src/auxkernels/MDNParticleAux.C b/src/auxkernels/MDNParticleAux.C
new file mode 100644
index 00000000..fc965d6c
--- /dev/null
+++ b/src/auxkernels/MDNParticleAux.C
@@ -0,0 +1,34 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "MDNParticleAux.h"
+#include "MDRunBase.h"
+
+registerMooseObject("MagpieApp", MDNParticleAux);
+
+template <>
+InputParameters
+validParams<MDNParticleAux>()
+{
+  InputParameters params = validParams<AuxKernel>();
+  params.addRequiredParam<UserObjectName>("user_object", "Name of MD runner UserObject");
+  params.addClassDescription("Injects the number of MD particles from MDRunBase object user_object "
+                             "into auxiliary variable.");
+  return params;
+}
+
+MDNParticleAux::MDNParticleAux(const InputParameters & parameters)
+  : AuxKernel(parameters), _md_uo(getUserObject<MDRunBase>("user_object"))
+{
+}
+
+Real
+MDNParticleAux::computeValue()
+{
+  return _md_uo.elemParticles(_current_elem->id()).size();
+}
diff --git a/src/userobjects/LAMMPSFileRunner.C b/src/userobjects/LAMMPSFileRunner.C
new file mode 100644
index 00000000..86c30060
--- /dev/null
+++ b/src/userobjects/LAMMPSFileRunner.C
@@ -0,0 +1,332 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "LAMMPSFileRunner.h"
+#include "MooseUtils.h"
+#include <sstream>
+
+registerMooseObject("MagpieApp", LAMMPSFileRunner);
+
+template <>
+InputParameters
+validParams<LAMMPSFileRunner>()
+{
+  InputParameters params = validParams<MDRunBase>();
+  params.addClassDescription("Reads a LAMMPS dump file or sequence of dump files and "
+                             "maps LAMMPS particles to MOOSE FEM mesh.");
+  params.addParam<bool>("time_sequence",
+                        false,
+                        "If true, a sequence of dump files is read, "
+                        "else a single dump file is read.");
+  params.addRequiredParam<FileName>("lammps_file",
+                                    "Name of LAMMPS dump file or file base of LAMMPS file for "
+                                    "a sequence.");
+  params.addRequiredParam<std::vector<unsigned int>>(
+      "xyz_columns", "Column ids of the x, y, and z coordinates of particles.");
+
+  params.addParam<std::vector<unsigned int>>(
+      "vel_columns", "Column ids of the vx, vy, and vz velocity components of particles.");
+
+  params.addParam<FunctionName>("time_conversion",
+                                "A conversion from FEM simulation time to MD time stamps.");
+  params.addClassDescription("Allows coupling FEM calculations to LAMMPS dump files.");
+  return params;
+}
+
+LAMMPSFileRunner::LAMMPSFileRunner(const InputParameters & parameters)
+  : MDRunBase(parameters),
+    _time_sequence(getParam<bool>("time_sequence")),
+    _lammps_file(getParam<FileName>("lammps_file")),
+    _pos_columns(getParam<std::vector<unsigned int>>("xyz_columns")),
+    _velocity_from_file(isParamValid("vel_columns")),
+    _time_conversion(nullptr)
+{
+  ;
+  if (_pos_columns.size() != _dim)
+    mooseError(
+        "Dimensionality is ", _dim, " but xyz_columns has ", _pos_columns.size(), " entries.");
+
+  if (_velocity_from_file)
+  {
+    _vel_columns = getParam<std::vector<unsigned int>>("vel_columns");
+    if (_vel_columns.size() != _dim)
+      mooseError(
+          "Dimensionality is ", _dim, " but vel_columns has ", _vel_columns.size(), " entries.");
+  }
+
+  if (isParamValid("time_conversion"))
+    _time_conversion = &getFunction("time_conversion");
+}
+
+void
+LAMMPSFileRunner::initialSetup()
+{
+  // call MDRunBase initialSetup
+  MDRunBase::initialSetup();
+
+  if (_time_sequence)
+    return;
+
+  readLAMMPSFile(_lammps_file);
+  updateKDTree();
+}
+
+void
+LAMMPSFileRunner::updateParticleList()
+{
+  // read LAMMPS file if a time sequence is read
+  if (_time_sequence)
+  {
+    // map to MD time
+    Real md_time = _t;
+    if (_time_conversion)
+      md_time = _time_conversion->value(_t, Point());
+
+    // find files w/ timestamps "bracketing" md_time
+    std::pair<std::string, std::string> files;
+    std::pair<Real, Real> timestamps;
+    findBracketingLAMMPSFiles(md_time, files, timestamps);
+    readLAMMPSFileHistory(files, timestamps, md_time);
+    updateKDTree();
+  }
+
+  // update the mapping to mesh if mesh has changed or ...
+  mapMDParticles();
+}
+
+void
+LAMMPSFileRunner::findBracketingLAMMPSFiles(Real md_time,
+                                            std::pair<std::string, std::string> & filenames,
+                                            std::pair<Real, Real> & timestamps)
+{
+  /// the return values
+  std::string before;
+  std::string after;
+
+  /// extract path and file_base from provided lammps file entry
+  std::pair<std::string, std::string> p = MooseUtils::splitFileName(_lammps_file);
+
+  /// get all files in the directory
+  std::list<std::string> files = MooseUtils::getFilesInDirs({p.first});
+
+  /// loop over all files, retain the ones that begin with base and extract their timestamps
+  std::vector<Real> times;
+  times.reserve(files.size());
+
+  unsigned int lammps_files_found = 0;
+  Real closest_after_time = std::numeric_limits<Real>::max();
+  Real closest_before_time = -std::numeric_limits<Real>::max();
+  for (auto & file : files)
+  {
+    std::pair<std::string, std::string> f = MooseUtils::splitFileName(file);
+    if (f.second.compare(0, p.second.size(), p.second))
+      continue;
+
+    // split the file at dots and check the second item ==> timestamp
+    Real timestamp;
+    std::vector<std::string> elements;
+    MooseUtils::tokenize(f.second, elements, 1, ".");
+    if (elements.size() < 2)
+      mooseError("Tokenizing filename ",
+                 f.second,
+                 " failed. LAMMPS filename must be base.<t>.xyz where <t> is the timestamp.");
+    std::stringstream sstr(elements[1]);
+    sstr >> timestamp;
+
+    // increase the counter & check if this file is a candidate for before/after
+    ++lammps_files_found;
+    if (timestamp > md_time && timestamp < closest_after_time)
+    {
+      closest_after_time = timestamp;
+      after = file;
+    }
+    else if (timestamp <= md_time && timestamp > closest_before_time)
+    {
+      closest_before_time = timestamp;
+      before = file;
+    }
+  }
+
+  if (lammps_files_found == 0)
+    mooseError("No LAMMPS file with name base ", p.second, " in folder ", p.first);
+
+  // guard against the case that only before or after were found => set the two equal
+  if (before.size() == 0)
+    before = after;
+  else if (after.size() == 0)
+    after = before;
+
+  filenames = std::pair<std::string, std::string>(before, after);
+  timestamps = std::pair<Real, Real>(closest_before_time, closest_after_time);
+}
+
+void
+LAMMPSFileRunner::readLAMMPSFile(FileName filename)
+{
+  // check if this file is readable
+  MooseUtils::checkFileReadable(filename);
+
+  // read LAMMPS file
+  std::ifstream file(filename);
+
+  // skip first three lines
+  for (unsigned int j = 0; j < 3; ++j)
+    file.ignore(10000, '\n');
+
+  // get number of particles
+  file >> _n_particles;
+
+  // skip another five lines
+  for (unsigned int j = 0; j < 6; ++j)
+    file.ignore(10000, '\n');
+
+  // attach the right _md_particles.pos
+  std::vector<Point> & position = _md_particles.pos;
+
+  // zero out the position & id vector
+  position.clear();
+  _md_particles.id.clear();
+
+  // size the md particle vector considering not all processors have all particles
+  position.reserve(std::ceil(2.0 * _n_particles / _nproc));
+  _md_particles.id.reserve(std::ceil(2.0 * _n_particles / _nproc));
+
+  // read particle entries
+  _n_local_particles = 0;
+  for (unsigned int j = 0; j < _n_particles; ++j)
+  {
+    // check if this particle is in the local bounding box
+    std::string line;
+    std::getline(file, line);
+    std::vector<std::string> elements;
+    MooseUtils::tokenize<std::string>(line, elements, 1, " ");
+
+    Point pos;
+    for (unsigned int i = 0; i < _dim; ++i)
+    {
+      std::stringstream sstr(elements[_pos_columns[i]]);
+      sstr >> pos(i);
+    }
+
+    // check if this particle is in this processor BBox
+    if (!_bbox[processor_id()].contains_point(pos))
+      continue;
+
+    ++_n_local_particles;
+    position.push_back(pos);
+    _md_particles.id.push_back(j);
+  }
+  file.close();
+
+  // size back
+  position.shrink_to_fit();
+  _md_particles.id.shrink_to_fit();
+}
+
+void
+LAMMPSFileRunner::readLAMMPSFileHistory(std::pair<FileName, FileName> filenames,
+                                        std::pair<Real, Real> timestamps,
+                                        Real current_time)
+{
+  // check if this file is readable
+  MooseUtils::checkFileReadable(filenames.first);
+  MooseUtils::checkFileReadable(filenames.second);
+
+  // open LAMMPS file
+  std::ifstream file_before(filenames.first);
+  std::ifstream file_after(filenames.second);
+
+  // compute weights for particle
+  Real weight = (current_time - timestamps.first) / (timestamps.second - timestamps.first);
+
+  // skip first three lines
+  for (unsigned int j = 0; j < 3; ++j)
+  {
+    file_before.ignore(10000, '\n');
+    file_after.ignore(10000, '\n');
+  }
+
+  // get number of particles
+  unsigned int np;
+  file_before >> _n_particles;
+  file_after >> np;
+  if (np != _n_particles)
+    mooseError("Files ",
+               filenames.first,
+               " and ",
+               filenames.second,
+               " have different number of particles.");
+
+  // skip another five lines
+  for (unsigned int j = 0; j < 6; ++j)
+  {
+    file_before.ignore(10000, '\n');
+    file_after.ignore(10000, '\n');
+  }
+
+  // attach the right _md_particles.pos
+  std::vector<Point> & position = _md_particles.pos;
+
+  // zero out the position & id vector
+  position.clear();
+  _md_particles.id.clear();
+
+  // size the md particle vector considering not all processors have all particles
+  position.reserve(std::ceil(2.0 * _n_particles / _nproc));
+  _md_particles.id.reserve(std::ceil(2.0 * _n_particles / _nproc));
+
+  // read particle entries
+  _n_local_particles = 0;
+  for (unsigned int j = 0; j < _n_particles; ++j)
+  {
+    // check if this particle is in the local bounding box
+    // either in before or after state
+
+    // get the MD particle from the before file
+    std::vector<std::string> elements;
+    std::string line_before;
+    std::getline(file_before, line_before);
+    MooseUtils::tokenize<std::string>(line_before, elements, 1, " ");
+
+    Point pos_before;
+    for (unsigned int i = 0; i < _dim; ++i)
+    {
+      std::stringstream sstr(elements[_pos_columns[i]]);
+      sstr >> pos_before(i);
+    }
+
+    // get the MD particle from the after file
+    elements.clear();
+    std::string line_after;
+    std::getline(file_after, line_after);
+
+    MooseUtils::tokenize<std::string>(line_after, elements, 1, " ");
+
+    Point pos_after;
+    for (unsigned int i = 0; i < _dim; ++i)
+    {
+      std::stringstream sstr(elements[_pos_columns[i]]);
+      sstr >> pos_after(i);
+    }
+
+    // check if this particle is in this processor BBox
+    if (!_bbox[processor_id()].contains_point(pos_after) &&
+        !_bbox[processor_id()].contains_point(pos_before))
+      continue;
+
+    ++_n_local_particles;
+    position.push_back((1 - weight) * pos_before + weight * pos_after);
+    _md_particles.id.push_back(j);
+  }
+  file_before.close();
+  file_after.close();
+
+  // size back
+  position.shrink_to_fit();
+  _md_particles.id.shrink_to_fit();
+}
diff --git a/src/userobjects/MDRunBase.C b/src/userobjects/MDRunBase.C
new file mode 100644
index 00000000..f2489af4
--- /dev/null
+++ b/src/userobjects/MDRunBase.C
@@ -0,0 +1,131 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "MDRunBase.h"
+#include "MooseMesh.h"
+
+#include "libmesh/mesh_tools.h"
+
+// custom data load and data store methods for MDParticles class
+template <>
+inline void
+dataStore(std::ostream & stream, MDRunBase::MDParticles & pl, void * context)
+{
+  dataStore(stream, pl.pos, context);
+  dataStore(stream, pl.vel, context);
+  dataStore(stream, pl.id, context);
+  dataStore(stream, pl.elem_id, context);
+}
+
+template <>
+inline void
+dataLoad(std::istream & stream, MDRunBase::MDParticles & pl, void * context)
+{
+  dataLoad(stream, pl.pos, context);
+  dataLoad(stream, pl.vel, context);
+  dataLoad(stream, pl.id, context);
+  dataLoad(stream, pl.elem_id, context);
+}
+
+template <>
+InputParameters
+validParams<MDRunBase>()
+{
+  InputParameters params = validParams<GeneralUserObject>();
+  params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
+  params.suppressParameter<ExecFlagEnum>("execute_on");
+
+  params.addClassDescription(
+      "Base class for execution of coupled molecular dynamics MOOSE calculations.");
+  return params;
+}
+
+MDRunBase::MDRunBase(const InputParameters & parameters)
+  : GeneralUserObject(parameters),
+    _mesh(_subproblem.mesh()),
+    _dim(_mesh.dimension()),
+    _nproc(_app.n_processors()),
+    _bbox(_nproc)
+{
+}
+
+void
+MDRunBase::initialSetup()
+{
+  for (unsigned int j = 0; j < _nproc; ++j)
+    _bbox[j] = MeshTools::create_processor_bounding_box(_mesh, j);
+}
+
+void
+MDRunBase::timestepSetup()
+{
+  // update/init subdomain bounding boxes
+  for (unsigned int j = 0; j < _nproc; ++j)
+    _bbox[j] = MeshTools::create_processor_bounding_box(_mesh, j);
+
+  // callback for updating md particle list
+  updateParticleList();
+}
+
+void
+MDRunBase::updateKDTree()
+{
+  _kd_tree = libmesh_make_unique<KDTree>(_md_particles.pos, 50);
+}
+
+const std::vector<unsigned int>
+MDRunBase::elemParticles(dof_id_type elem_id) const
+{
+  if (_elem_particles.find(elem_id) != _elem_particles.end())
+    return _elem_particles.find(elem_id)->second;
+  return std::vector<unsigned int>(0);
+}
+
+void
+MDRunBase::mapMDParticles()
+{
+  // clear data structures
+  _elem_particles.clear();
+  _md_particles.elem_id.assign(_md_particles.pos.size(), libMesh::invalid_uint);
+
+  // loop over elements to check against bounding box
+  ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
+  ConstElemRange::const_iterator el = elem_range.begin();
+  for (el = elem_range.begin(); el != elem_range.end(); ++el)
+  {
+    const Elem * elem = *el;
+
+    // find all points within an inflated bounding box
+    std::vector<std::pair<std::size_t, Real>> indices_dist;
+    BoundingBox bbox = elem->loose_bounding_box();
+    Point center = 0.5 * (bbox.min() + bbox.max());
+    Real radius = (bbox.max() - center).norm();
+    _kd_tree->radiusSearch(center, radius, indices_dist);
+
+    for (auto & p : indices_dist)
+    {
+      Point candidate = _md_particles.pos[p.first];
+
+      // sort out whatever is not in the bbox & avoid double counting of elements
+      // on element boundaries (first come first serve for elems)
+      if (_md_particles.elem_id[p.first] != libMesh::invalid_uint ||
+          !bbox.contains_point(candidate))
+        continue;
+
+      // serious candidates need to be checked using expensive contains_point
+      if (elem->contains_point(candidate))
+      {
+        _md_particles.elem_id[p.first] = elem->id();
+        if (_elem_particles.find(elem->id()) != _elem_particles.end())
+          _elem_particles[elem->id()].push_back(elem->id());
+        else
+          _elem_particles[elem->id()] = {elem->id()};
+      }
+    }
+  }
+}
diff --git a/tests/userobjects/lammps_file/gold/lammps_file_runner_out.e b/tests/userobjects/lammps_file/gold/lammps_file_runner_out.e
new file mode 100644
index 00000000..532cf0b0
Binary files /dev/null and b/tests/userobjects/lammps_file/gold/lammps_file_runner_out.e differ
diff --git a/tests/userobjects/lammps_file/gold/lammps_file_runner_sequence_out.e b/tests/userobjects/lammps_file/gold/lammps_file_runner_sequence_out.e
new file mode 100644
index 00000000..4928abc8
Binary files /dev/null and b/tests/userobjects/lammps_file/gold/lammps_file_runner_sequence_out.e differ
diff --git a/tests/userobjects/lammps_file/lammps_file_runner.i b/tests/userobjects/lammps_file/lammps_file_runner.i
new file mode 100644
index 00000000..0e03c4db
--- /dev/null
+++ b/tests/userobjects/lammps_file/lammps_file_runner.i
@@ -0,0 +1,54 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  xmin = -1
+  xmax = 2
+  nx = 3
+  ymin = -2
+  ymax = 1
+  ny = 3
+  zmin = 0
+  zmax = 2
+  nz = 2
+[]
+
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [./d]
+  [../]
+[]
+
+[AuxVariables]
+  [./nparticles]
+    family = MONOMIAL
+    order = CONSTANT
+  [../]
+[]
+
+[AuxKernels]
+  [./np]
+    type = MDNParticleAux
+    variable = nparticles
+    user_object = lammps_runner
+  [../]
+[]
+
+[UserObjects]
+  [./lammps_runner]
+    type = LAMMPSFileRunner
+    lammps_file = 'simple.0.xyz'
+    xyz_columns = '0 1 2'
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 1
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/userobjects/lammps_file/lammps_file_runner_sequence.i b/tests/userobjects/lammps_file/lammps_file_runner_sequence.i
new file mode 100644
index 00000000..f6697928
--- /dev/null
+++ b/tests/userobjects/lammps_file/lammps_file_runner_sequence.i
@@ -0,0 +1,64 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  xmin = 0
+  xmax = 2
+  nx = 2
+  ymin = 0
+  ymax = 2
+  ny = 2
+  zmin = 0
+  zmax = 2
+  nz = 2
+[]
+
+[Problem]
+  kernel_coverage_check = false
+[]
+
+[Variables]
+  [./d]
+  [../]
+[]
+
+[AuxVariables]
+  [./nparticles]
+    family = MONOMIAL
+    order = CONSTANT
+  [../]
+[]
+
+[AuxKernels]
+  [./np]
+    type = MDNParticleAux
+    variable = nparticles
+    user_object = lammps_runner
+  [../]
+[]
+
+[UserObjects]
+  [./lammps_runner]
+    type = LAMMPSFileRunner
+    lammps_file = 'sequence/simple'
+    time_sequence = true
+    xyz_columns = '0 1 2'
+    time_conversion = tc
+  [../]
+[]
+
+[Functions]
+  [./tc]
+    type = ParsedFunction
+    value = '2 * t + 5'
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  end_time = 3
+  dt = 0.1
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/tests/userobjects/lammps_file/sequence/simple.11.xyz b/tests/userobjects/lammps_file/sequence/simple.11.xyz
new file mode 100644
index 00000000..e9a54d35
--- /dev/null
+++ b/tests/userobjects/lammps_file/sequence/simple.11.xyz
@@ -0,0 +1,10 @@
+ITEM: TIMESTEP
+11
+ITEM: NUMBER OF ATOMS
+1
+ITEM: BOX BOUNDS pp pp pp
+ 0.0 2.0
+ 0.0 2.0
+ 0.0 2.0
+ITEM: ATOMS x y z
+ 1.5 1.5 1.5
diff --git a/tests/userobjects/lammps_file/sequence/simple.5.xyz b/tests/userobjects/lammps_file/sequence/simple.5.xyz
new file mode 100644
index 00000000..3a8f893a
--- /dev/null
+++ b/tests/userobjects/lammps_file/sequence/simple.5.xyz
@@ -0,0 +1,10 @@
+ITEM: TIMESTEP
+5
+ITEM: NUMBER OF ATOMS
+1
+ITEM: BOX BOUNDS pp pp pp
+ 0.0 2.0
+ 0.0 2.0
+ 0.0 2.0
+ITEM: ATOMS x y z
+ 0.5 0.5 0.5
diff --git a/tests/userobjects/lammps_file/sequence/simple.7.xyz b/tests/userobjects/lammps_file/sequence/simple.7.xyz
new file mode 100644
index 00000000..6bf541fa
--- /dev/null
+++ b/tests/userobjects/lammps_file/sequence/simple.7.xyz
@@ -0,0 +1,10 @@
+ITEM: TIMESTEP
+7
+ITEM: NUMBER OF ATOMS
+1
+ITEM: BOX BOUNDS pp pp pp
+ 0.0 2.0
+ 0.0 2.0
+ 0.0 2.0
+ITEM: ATOMS x y z
+ 1.5 0.5 0.5
diff --git a/tests/userobjects/lammps_file/sequence/simple.9.xyz b/tests/userobjects/lammps_file/sequence/simple.9.xyz
new file mode 100644
index 00000000..582540b1
--- /dev/null
+++ b/tests/userobjects/lammps_file/sequence/simple.9.xyz
@@ -0,0 +1,10 @@
+ITEM: TIMESTEP
+9
+ITEM: NUMBER OF ATOMS
+1
+ITEM: BOX BOUNDS pp pp pp
+ 0.0 2.0
+ 0.0 2.0
+ 0.0 2.0
+ITEM: ATOMS x y z
+ 1.5 1.5 0.5
diff --git a/tests/userobjects/lammps_file/simple.0.xyz b/tests/userobjects/lammps_file/simple.0.xyz
new file mode 100644
index 00000000..5d6d7925
--- /dev/null
+++ b/tests/userobjects/lammps_file/simple.0.xyz
@@ -0,0 +1,14 @@
+ITEM: TIMESTEP
+0
+ITEM: NUMBER OF ATOMS
+5
+ITEM: BOX BOUNDS pp pp pp
+-1.0000000000000000e+00 2.0000000000000000e+00
+-2.0000000000000000e+00 1.0000000000000000e+00
+ 0.0000000000000000e+00 2.0000000000000000e+00
+ITEM: ATOMS x y z
+-0.5 -1.5 0.1
+ 0.1  0.9 0.9
+ 0.99 0.8 0.1
+ 0.99 0.8 1.1
+ 2.1  0.0 0.0
diff --git a/tests/userobjects/lammps_file/tests b/tests/userobjects/lammps_file/tests
new file mode 100644
index 00000000..ded569df
--- /dev/null
+++ b/tests/userobjects/lammps_file/tests
@@ -0,0 +1,14 @@
+[Tests]
+  [./lammps_file_runner]
+    type = 'Exodiff'
+    input = 'lammps_file_runner.i'
+    exodiff = 'lammps_file_runner_out.e'
+  [../]
+
+  [./lammps_file_runner_sequence]
+    type = 'Exodiff'
+    input = 'lammps_file_runner_sequence.i'
+    exodiff = 'lammps_file_runner_sequence_out.e'
+    max_parallel = 1
+  [../]
+[]