diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c
index 9d7238f..07ed0e5 100644
--- a/gromacs/force_lj.c
+++ b/gromacs/force_lj.c
@@ -61,6 +61,7 @@ double computeForceLJ_ref(Parameter *param, Atom *atom, Neighbor *neighbor, Stat
 
         for(int k = 0; k < numneighs; k++) {
             int cj = neighs[k];
+            int any = 0;
             MD_FLOAT *cjptr = cluster_pos_ptr(cj);
             for(int cii = 0; cii < CLUSTER_DIM_M; cii++) {
                 MD_FLOAT xtmp = cluster_x(ciptr, cii);
@@ -83,6 +84,7 @@ double computeForceLJ_ref(Parameter *param, Atom *atom, Neighbor *neighbor, Stat
                             fix += delx * force;
                             fiy += dely * force;
                             fiz += delz * force;
+                            any = 1;
                             addStat(stats->atoms_within_cutoff, 1);
                         } else {
                             addStat(stats->atoms_outside_cutoff, 1);
@@ -90,6 +92,12 @@ double computeForceLJ_ref(Parameter *param, Atom *atom, Neighbor *neighbor, Stat
                     }
                 }
 
+                if(any != 0) {
+                    addStat(stats->clusters_within_cutoff, 1);
+                } else {
+                    addStat(stats->clusters_outside_cutoff, 1);
+                }
+
                 cluster_x(cifptr, cii) += fix;
                 cluster_y(cifptr, cii) += fiy;
                 cluster_z(cifptr, cii) += fiz;
diff --git a/gromacs/includes/stats.h b/gromacs/includes/stats.h
index 4b08e3d..df022d0 100644
--- a/gromacs/includes/stats.h
+++ b/gromacs/includes/stats.h
@@ -31,6 +31,8 @@ typedef struct {
     long long int force_iters;
     long long int atoms_within_cutoff;
     long long int atoms_outside_cutoff;
+    long long int clusters_within_cutoff;
+    long long int clusters_outside_cutoff;
 } Stats;
 
 void initStats(Stats *s);
diff --git a/gromacs/stats.c b/gromacs/stats.c
index ccc6030..5f62918 100644
--- a/gromacs/stats.c
+++ b/gromacs/stats.c
@@ -11,6 +11,8 @@ void initStats(Stats *s) {
     s->force_iters = 0;
     s->atoms_within_cutoff = 0;
     s->atoms_outside_cutoff = 0;
+    s->clusters_within_cutoff = 0;
+    s->clusters_outside_cutoff = 0;
 }
 
 void displayStatistics(Atom *atom, Parameter *param, Stats *stats, double *timer) {
@@ -40,8 +42,10 @@ void displayStatistics(Atom *atom, Parameter *param, Stats *stats, double *timer
     printf("\tCycles/SIMD iteration: %.4f\n", timer[FORCE] * param->proc_freq * 1e9 / stats->force_iters);
 
     #ifdef USE_REFERENCE_VERSION
-    const double eff_pct = (double)stats->atoms_within_cutoff / (double)(stats->atoms_within_cutoff + stats->atoms_outside_cutoff) * 100.0;
-    printf("\tAtoms within/outside cutoff radius: %lld/%lld (%.2f%%)\n", stats->atoms_within_cutoff, stats->atoms_outside_cutoff, eff_pct);
+    const double atoms_eff = (double)stats->atoms_within_cutoff / (double)(stats->atoms_within_cutoff + stats->atoms_outside_cutoff) * 100.0;
+    printf("\tAtoms within/outside cutoff radius: %lld/%lld (%.2f%%)\n", stats->atoms_within_cutoff, stats->atoms_outside_cutoff, atoms_eff);
+    const double clusters_eff = (double)stats->clusters_within_cutoff / (double)(stats->clusters_within_cutoff + stats->clusters_outside_cutoff) * 100.0;
+    printf("\tClusters within/outside cutoff radius: %lld/%lld (%.2f%%)\n", stats->clusters_within_cutoff, stats->clusters_outside_cutoff, clusters_eff);
     #endif
 
 #endif