diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst
index 745d898d01..0a44ccb197 100644
--- a/doc/src/pair_granular.rst
+++ b/doc/src/pair_granular.rst
@@ -282,7 +282,9 @@ stress components, total surface area involved in
 contact, and individual contact areas. In the input script, these quantities are
 initialized by calling *run 0* and can then be accessed using subsequent *compute*
 commands. The last *compute* command uses *pair/local p13* to calculate the pairwise
-contact areas for each active contact in the *group-ID*.
+contact areas for each active contact in the *group-ID*. Due to the use of an apparent
+radius in the *mdr* model, the keyword/arg pair *cutoff radius* must be specified for
+*pair/local* to properly detect existing contacts.
 
 .. code-block:: LAMMPS
 
@@ -292,7 +294,7 @@ contact areas for each active contact in the *group-ID*.
    compute ID group-ID property/atom d_sigmayy
    compute ID group-ID property/atom d_sigmazz
    compute ID group-ID property/atom d_Acon1
-   compute ID group-ID pair/local p13
+   compute ID group-ID pair/local p13 cutoff radius
 
 .. note::
 
diff --git a/examples/granular/in.triaxial.compaction.12 b/examples/granular/in.triaxial.compaction.12
index 6befe55241..eee9a2bfd8 100644
--- a/examples/granular/in.triaxial.compaction.12
+++ b/examples/granular/in.triaxial.compaction.12
@@ -98,7 +98,7 @@ dump dumpParticles all custom ${output_rate} triaxial_compaction_12.dump id type
 run 0
 
 # print out contact area evolution for particles 1 and 12
-compute Ac_1_12 particles_1_12 pair/local p13
+compute Ac_1_12 particles_1_12 pair/local p13 cutoff radius
 compute Ac_1_12_sum particles_1_12 reduce sum c_Ac_1_12 inputs local
 variable Ac_1_12 equal c_Ac_1_12_sum
 fix logArea all print 100 "${plane_disp} ${Ac_1_12}" file pair_1_12_contact_area_triaxial12.csv screen no
diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp
index 79b6c2406e..bae3371ebc 100644
--- a/src/GRANULAR/pair_granular.cpp
+++ b/src/GRANULAR/pair_granular.cpp
@@ -728,6 +728,10 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   model->xj = x[j];
   model->radi = radius[i];
   model->radj = radius[j];
+  model->i = i;
+  model->j = j;
+  model->itype = itype;
+  model->jtype = jtype;
   model->history_update = 0; // Don't update history
 
   // If history is needed