diff --git a/doc/src/Howto_granular.rst b/doc/src/Howto_granular.rst
index c22cab66bc..b0c801be11 100644
--- a/doc/src/Howto_granular.rst
+++ b/doc/src/Howto_granular.rst
@@ -45,10 +45,15 @@ atoms, and should be used for granular system instead of the fix style
 
 To model heat conduction, one must add the temperature and heatflow
 atom variables with:
+
 * :doc:`fix property/atom <fix_property_atom>`
+
 a temperature integration fix
+
 * :doc:`fix heat/flow <fix_heat_flow>`
+
 and a heat conduction option defined in both
+
 * :doc:`pair_style granular <pair_granular>`
 * :doc:`fix wall/gran <fix_wall_gran>`
 
diff --git a/doc/src/compute_property_atom.rst b/doc/src/compute_property_atom.rst
index b03d6eb74e..4484c3b697 100644
--- a/doc/src/compute_property_atom.rst
+++ b/doc/src/compute_property_atom.rst
@@ -23,8 +23,9 @@ Syntax
                              spx, spy, spz, sp, fmx, fmy, fmz,
                              nbonds,
                              radius, diameter, omegax, omegay, omegaz,
+                             temperature, heatflow,
                              angmomx, angmomy, angmomz,
-                             shapex,shapey, shapez,
+                             shapex, shapey, shapez,
                              quatw, quati, quatj, quatk, tqx, tqy, tqz,
                              end1x, end1y, end1z, end2x, end2y, end2z,
                              corner1x, corner1y, corner1z,
@@ -56,6 +57,8 @@ Syntax
            *nbonds* = number of bonds assigned to an atom
            *radius,diameter* = radius,diameter of spherical particle
            *omegax,omegay,omegaz* = angular velocity of spherical particle
+           *temperature* = internal temperature of spherical particle
+           *heatflow* = internal heat flow of spherical particle
            *angmomx,angmomy,angmomz* = angular momentum of aspherical particle
            *shapex,shapey,shapez* = 3 diameters of aspherical particle
            *quatw,quati,quatj,quatk* = quaternion components for aspherical or body particles
diff --git a/doc/src/dump.rst b/doc/src/dump.rst
index 2d1598e493..82faba6c81 100644
--- a/doc/src/dump.rst
+++ b/doc/src/dump.rst
@@ -104,7 +104,6 @@ Syntax
                                q, mux, muy, muz, mu,
                                radius, diameter, omegax, omegay, omegaz,
                                angmomx, angmomy, angmomz, tqx, tqy, tqz,
-                               heatflow, temperature,
                                c_ID, c_ID[I], f_ID, f_ID[I], v_name,
                                i_name, d_name, i2_name[I], d2_name[I]
 
@@ -131,8 +130,6 @@ Syntax
            omegax,omegay,omegaz = angular velocity of spherical particle
            angmomx,angmomy,angmomz = angular momentum of aspherical particle
            tqx,tqy,tqz = torque on finite-size particles
-           heatflow = rate of heat flow into particle
-           temperature = temperature of particle
            c_ID = per-atom vector calculated by a compute with ID
            c_ID[I] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below)
            f_ID = per-atom vector calculated by a fix with ID
diff --git a/doc/src/fix_property_atom.rst b/doc/src/fix_property_atom.rst
index b177fe9a96..d20358b7a7 100644
--- a/doc/src/fix_property_atom.rst
+++ b/doc/src/fix_property_atom.rst
@@ -22,6 +22,8 @@ Syntax
        *mol* = molecule IDs
        *q* = charge
        *rmass* = per-atom mass
+       *temperature* = internal temperature of atom
+       *heatflow* = internal heat flow of atom
        i_name = new integer vector referenced by name
        d_name = new floating-point vector referenced by name
        i2_name = new integer array referenced by name
@@ -59,14 +61,18 @@ these properties for each atom in the system when a data file is read.
 This fix augments the set of per-atom properties with new custom
 ones. This can be useful in several scenarios.
 
-If the atom style does not define molecule IDs, per-atom charge, or
-per-atom mass, they can be added using the *mol*\ , *q* or *rmass*
+If the atom style does not define molecule IDs, per-atom charge,
+per-atom mass, internal temperature, or internal heat flow, they can
+be added using the *mol*\ , *q*, *rmass*, *temperature*, or *heatflow*
 keywords.  This could be useful to define "molecules" to use as rigid
 bodies with the :doc:`fix rigid <fix_rigid>` command, or to carry
 around an extra flag with atoms (stored as a molecule ID) that can be
 used by various commands like :doc:`compute chunk/atom
 <compute_chunk_atom>` to group atoms without having to use the group
 command (which is limited to a total of 32 groups including *all*\ ).
+For finite-size particles, an internal temperature and heat flow can
+be used to model heat conduction as in the
+:doc:`GRANULAR package <Howto_granular>`.
 
 Another application is to use the *rmass* flag in order to have
 per-atom masses instead of per-type masses.  This could be used to
@@ -85,9 +91,10 @@ properties that are not needed such as bond lists, which incurs some
 overhead when there are no bonds.
 
 In the future, we may add additional existing per-atom properties to
-fix property/atom, similar to *mol*\ , *q* or *rmass*\ , which
-"turn-on" specific properties defined by some atom styles, so they can
-be easily used by atom styles that do not define them.
+fix property/atom, similar to *mol*\ , *q*, *rmass*\ , *temperature*\ ,
+or *heatflow* which "turn-on" specific properties defined by some atom
+styles, so they can be easily used by atom styles that do not define
+them.
 
 More generally, the *i_name* and *d_name* options allow one or more
 new custom per-atom vectors to be defined.  Likewise the *i2_name* and
diff --git a/examples/granular/in.pour.heat b/examples/granular/in.pour.heat
index 907e56dc39..cc6b03f7d0 100644
--- a/examples/granular/in.pour.heat
+++ b/examples/granular/in.pour.heat
@@ -73,7 +73,8 @@ thermo          100
 
 timestep        0.001
 
-#dump           1 all custom 1000 ${name}.dump id type radius mass x y z temperature heatflow
+compute        1 all property/atom temperature heatflow
+#dump           1 all custom 1000 ${name}.dump id type radius mass x y z c_1[*]
 
 run             100000
 
diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp
index b484df7fab..f1482d4203 100644
--- a/src/BPM/bond_bpm.cpp
+++ b/src/BPM/bond_bpm.cpp
@@ -357,7 +357,6 @@ void BondBPM::process_broken(int i, int j)
   if (i < nlocal) {
     for (m = 0; m < num_bond[i]; m++) {
       if (bond_atom[i][m] == tag[j]) {
-        bond_type[i][m] = 0;
         n = num_bond[i];
         bond_type[i][m] = bond_type[i][n - 1];
         bond_atom[i][m] = bond_atom[i][n - 1];
@@ -372,7 +371,6 @@ void BondBPM::process_broken(int i, int j)
   if (j < nlocal) {
     for (m = 0; m < num_bond[j]; m++) {
       if (bond_atom[j][m] == tag[i]) {
-        bond_type[j][m] = 0;
         n = num_bond[j];
         bond_type[j][m] = bond_type[j][n - 1];
         bond_atom[j][m] = bond_atom[j][n - 1];
diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp
index c3c101b995..b95b7267dc 100644
--- a/src/compute_property_atom.cpp
+++ b/src/compute_property_atom.cpp
@@ -205,6 +205,14 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
       if (!atom->omega_flag)
         error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
       pack_choice[i] = &ComputePropertyAtom::pack_omegaz;
+    } else if (strcmp(arg[iarg],"temperature") == 0) {
+      if (!atom->temperature_flag)
+        error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
+      pack_choice[i] = &ComputePropertyAtom::pack_temperature;
+    } else if (strcmp(arg[iarg],"heatflow") == 0) {
+      if (!atom->heatflow_flag)
+        error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
+      pack_choice[i] = &ComputePropertyAtom::pack_heatflow;
     } else if (strcmp(arg[iarg],"angmomx") == 0) {
       if (!atom->angmom_flag)
         error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
@@ -1213,6 +1221,36 @@ void ComputePropertyAtom::pack_omegaz(int n)
 
 /* ---------------------------------------------------------------------- */
 
+void ComputePropertyAtom::pack_temperature(int n)
+{
+  double *temperature = atom->temperature;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = temperature[i];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_heatflow(int n)
+{
+  double *heatflow = atom->heatflow;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) buf[n] = heatflow[i];
+    else buf[n] = 0.0;
+    n += nvalues;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputePropertyAtom::pack_angmomx(int n)
 {
   double **angmom = atom->angmom;
diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h
index 034b2901c2..d4f4db564f 100644
--- a/src/compute_property_atom.h
+++ b/src/compute_property_atom.h
@@ -95,6 +95,8 @@ class ComputePropertyAtom : public Compute {
   void pack_omegax(int);
   void pack_omegay(int);
   void pack_omegaz(int);
+  void pack_temperature(int);
+  void pack_heatflow(int);
   void pack_angmomx(int);
   void pack_angmomy(int);
   void pack_angmomz(int);
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index ed70c7413d..e9935b774f 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -41,7 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
      XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI,
      IX,IY,IZ,
      VX,VY,VZ,FX,FY,FZ,
-     Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,HEATFLOW,TEMPERATURE,
+     Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
      OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
      TQX,TQY,TQZ,
      COMPUTE,FIX,VARIABLE,IVEC,DVEC,IARRAY,DARRAY};
@@ -929,18 +929,6 @@ int DumpCustom::count()
         for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
         ptr = dchoose;
         nstride = 1;
-      } else if (thresh_array[ithresh] == HEATFLOW) {
-        if (!atom->heatflow_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
-        ptr = atom->heatflow;
-        nstride = 1;
-      } else if (thresh_array[ithresh] == TEMPERATURE) {
-        if (!atom->temperature_flag)
-          error->all(FLERR,
-                     "Threshold for an atom property that isn't allocated");
-        ptr = atom->temperature;
-        nstride = 1;
       } else if (thresh_array[ithresh] == OMEGAX) {
         if (!atom->omega_flag)
           error->all(FLERR,
@@ -1395,16 +1383,6 @@ int DumpCustom::parse_fields(int narg, char **arg)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
       pack_choice[iarg] = &DumpCustom::pack_diameter;
       vtype[iarg] = Dump::DOUBLE;
-    } else if (strcmp(arg[iarg],"heatflow") == 0) {
-      if (!atom->heatflow_flag)
-        error->all(FLERR,"Dumping an atom property that isn't allocated");
-      pack_choice[iarg] = &DumpCustom::pack_heatflow;
-      vtype[iarg] = Dump::DOUBLE;
-    } else if (strcmp(arg[iarg],"temperature") == 0) {
-      if (!atom->temperature_flag)
-        error->all(FLERR,"Dumping an atom property that isn't allocated");
-      pack_choice[iarg] = &DumpCustom::pack_temperature;
-      vtype[iarg] = Dump::DOUBLE;
     } else if (strcmp(arg[iarg],"omegax") == 0) {
       if (!atom->omega_flag)
         error->all(FLERR,"Dumping an atom property that isn't allocated");
@@ -1875,8 +1853,6 @@ int DumpCustom::modify_param(int narg, char **arg)
 
     else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
     else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER;
-    else if (strcmp(arg[1],"heatflow") == 0) thresh_array[nthresh] = HEATFLOW;
-    else if (strcmp(arg[1],"temperature") == 0) thresh_array[nthresh] = TEMPERATURE;
     else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX;
     else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY;
     else if (strcmp(arg[1],"omegaz") == 0) thresh_array[nthresh] = OMEGAZ;
@@ -2791,30 +2767,6 @@ void DumpCustom::pack_diameter(int n)
 
 /* ---------------------------------------------------------------------- */
 
-void DumpCustom::pack_heatflow(int n)
-{
-  double *heatflow = atom->heatflow;
-
-  for (int i = 0; i < nchoose; i++) {
-    buf[n] = heatflow[clist[i]];
-    n += size_one;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_temperature(int n)
-{
-  double *temperature = atom->temperature;
-
-  for (int i = 0; i < nchoose; i++) {
-    buf[n] = temperature[clist[i]];
-    n += size_one;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
 void DumpCustom::pack_omegax(int n)
 {
   double **omega = atom->omega;
diff --git a/src/dump_custom.h b/src/dump_custom.h
index 2b04944ec3..b600bd60b8 100644
--- a/src/dump_custom.h
+++ b/src/dump_custom.h
@@ -188,8 +188,6 @@ class DumpCustom : public Dump {
   void pack_mu(int);
   void pack_radius(int);
   void pack_diameter(int);
-  void pack_heatflow(int);
-  void pack_temperature(int);
 
   void pack_omegax(int);
   void pack_omegay(int);