diff --git a/src/RHEO/README b/src/RHEO/README
index 00fff2d694..7090fc828c 100644
--- a/src/RHEO/README
+++ b/src/RHEO/README
@@ -2,6 +2,6 @@ RHEO or Reproducing Hydrodynamics and Elastic Objects is a package to model mult
 systems. The authors include Joel Clemmer (Sandia), Thomas O'Connor (Carnegie Mellon), and
 Eric Palermo (Carnegie Mellon).
 
-This package requires the GNU scientific library (GSL) version 2.7 or later. To build this
-package, one must first separately install GSL in a location that can be found by your
-environment.
+This package requires the GNU scientific library (GSL). We recommend version 2.7 or later. To
+build this package, one must first separately install GSL in a location that can be found by
+your environment.
diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp
new file mode 100644
index 0000000000..a23ba5a639
--- /dev/null
+++ b/src/RHEO/compute_rheo_property_atom.cpp
@@ -0,0 +1,162 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+----------------------------------------------------------------------- */
+
+#include "compute_rheo_property_atom.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "error.h"
+#include "memory.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  index(nullptr), colindex(nullptr), pack_choice(nullptr)
+{
+  if (narg < 4)  utils::missing_cmd_args(FLERR, "compute property/atom", error);
+
+  peratom_flag = 1;
+  nvalues = narg - 3;
+  if (nvalues == 1) size_peratom_cols = 0;
+  else size_peratom_cols = nvalues;
+
+  // parse input values
+  // customize a new keyword by adding to if statement
+
+  pack_choice = new FnPtrPack[nvalues];
+
+  int i;
+  for (int iarg = 3; iarg < narg; iarg++) {
+    i = iarg-3;
+
+    if (strcmp(arg[iarg],"id") == 0) {
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_id;
+    } else if (strcmp(arg[iarg],"mol") == 0) {
+      if (!atom->molecule_flag)
+        error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_molecule;
+    } else if (strcmp(arg[iarg],"proc") == 0) {
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_proc;
+    } else if (strcmp(arg[iarg],"type") == 0) {
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_type;
+    } else if (strcmp(arg[iarg],"mass") == 0) {
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_mass;
+
+
+    } else {
+      error->all(FLERR,"Invalid keyword {} for compute rheo/property/atom command ", arg[iarg]);
+    }
+  }
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
+{
+  delete[] pack_choice;
+  memory->destroy(vector_atom);
+  memory->destroy(array_atom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOPropertyAtom::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow vector or array if necessary
+
+  if (atom->nmax > nmax) {
+    nmax = atom->nmax;
+    if (nvalues == 1) {
+      memory->destroy(vector_atom);
+      memory->create(vector_atom,nmax,"rheo/property/atom:vector");
+    } else {
+      memory->destroy(array_atom);
+      memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array");
+    }
+  }
+
+  // fill vector or array with per-atom values
+
+  if (nvalues == 1) {
+    buf = vector_atom;
+    (this->*pack_choice[0])(0);
+  } else {
+    if (nmax) buf = &array_atom[0][0];
+    else buf = nullptr;
+    for (int n = 0; n < nvalues; n++)
+      (this->*pack_choice[n])(n);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeRHEOPropertyAtom::memory_usage()
+{
+  double bytes = (double)nmax * nvalues * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   one method for every keyword compute rheo/property/atom can output
+   the atom property is packed into buf starting at n with stride nvalues
+   customize a new keyword by adding a method
+------------------------------------------------------------------------- */
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOPropertyAtom::pack_rho(int n)
+{
+  double *rho = atom->rho;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = rho[i];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOPropertyAtom::pack_status(int n)
+{
+  int *status = atom->status;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = status[i];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h
new file mode 100644
index 0000000000..26ca004da0
--- /dev/null
+++ b/src/RHEO/compute_rheo_property_atom.h
@@ -0,0 +1,68 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(rheo/property/atom,ComputeRHEOPropertyAtom);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
+#define LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeRHEOPropertyAtom : public Compute {
+ public:
+  ComputeRHEOPropertyAtom(class LAMMPS *, int, char **);
+  ~ComputeRHEOPropertyAtom() override;
+  void compute_peratom() override;
+  double memory_usage() override;
+
+ private:
+  int nvalues;
+  int nmax;
+  double *buf;
+
+  typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int);
+  FnPtrPack *pack_choice;    // ptrs to pack functions
+
+  void pack_rho(int);
+  void pack_drho(int);
+  void pack_temperature(int);
+  void pack_heatflow(int);
+  void pack_status(int);
+  void pack_phase(int);
+  void pack_surface(int);
+  void pack_r_surface(int);
+  void pack_divr_surface(int);
+  void pack_nx_surface(int);
+  void pack_ny_surface(int);
+  void pack_nz_surface(int);
+  void pack_coordination(int);
+  void pack_viscosity(int);
+  void pack_pressure(int);
+  void pack_conductivity(int);
+  void pack_cv(int);
+  void pack_vx_shift(int);
+  void pack_vy_shift(int);
+  void pack_vz_shift(int);
+
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h
index f936206811..a74696e68c 100644
--- a/src/RHEO/fix_rheo.h
+++ b/src/RHEO/fix_rheo.h
@@ -75,8 +75,8 @@ namespace RHEO_NS {
   enum Status{
     // Phase status
     STATUS_FLUID = 1 << 0,
-    STATUS_REACTIVE = 1 << 1,
-    STATUS_SOLID = 1 << 2,
+    STATUS_SOLID = 1 << 1,
+    STATUS_REACTIVE = 1 << 2,
     STATUS_FREEZING = 1 << 3,
     // Surface status
     STATUS_BULK = 1 << 4,
diff --git a/src/set.cpp b/src/set.cpp
index 3a175cbfe2..92033b772e 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -49,7 +49,7 @@ enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,
      DIPOLE,DIPOLE_RANDOM,SPIN_ATOM,SPIN_RANDOM,SPIN_ELECTRON,RADIUS_ELECTRON,
      QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,TEMPERATURE,
      DIAMETER,RADIUS_ATOM,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
-     SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
+     RHEO_STATUS,SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
      SMD_CONTACT_RADIUS,DPDTHETA,EPSILON,IVEC,DVEC,IARRAY,DARRAY};
 
 #define BIG INT_MAX
@@ -515,6 +515,24 @@ void Set::command(int narg, char **arg)
       topology(IMPROPER);
       iarg += 2;
 
+    } else if (strcmp(arg[iarg],"rheo/rho") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/rho", error);
+      if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
+      else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+      if (!atom->rho_flag)
+        error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
+      set(SPH_RHO);
+      iarg += 2;
+
+    } else if (strcmp(arg[iarg],"rheo/status") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error);
+      if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
+      else ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+      if (!atom->status_flag)
+        error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
+      set(RHEO_STATUS);
+      iarg += 2;
+
     } else if (strcmp(arg[iarg],"sph/e") == 0) {
       if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/e", error);
       if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
@@ -881,6 +899,13 @@ void Set::set(int keyword)
       if (dvalue <= 0.0) error->one(FLERR,"Invalid volume in set command");
       atom->vfrac[i] = dvalue;
     }
+
+    else if (keyword == RHEO_STATUS) {
+      if (ivalue != 0 && ivalue !=2)
+        error->one(FLERR,"Invalid value {} in set command for rheo/status", ivalue);
+      atom->status[i] = ivalue;
+    }
+
     else if (keyword == SPH_E) atom->esph[i] = dvalue;
     else if (keyword == SPH_CV) atom->cv[i] = dvalue;
     else if (keyword == SPH_RHO) atom->rho[i] = dvalue;