diff --git a/doc/src/fix_ave_correlate.rst b/doc/src/fix_ave_correlate.rst
index 659e15105c..fea927a073 100644
--- a/doc/src/fix_ave_correlate.rst
+++ b/doc/src/fix_ave_correlate.rst
@@ -32,13 +32,14 @@ Syntax
 
   .. parsed-literal::
 
-       *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full*
+       *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full* or *first*
          auto = correlate each value with itself
          upper = correlate each value with each succeeding value
          lower = correlate each value with each preceding value
          auto/upper = auto + upper
          auto/lower = auto + lower
          full = correlate each value with every other value, including itself = auto + upper + lower
+         first = correlate each value with the first value
        *ave* args = *one* or *running*
          one = zero the correlation accumulation every Nfreq steps
          running = accumulate correlations continuously
@@ -257,6 +258,9 @@ time.
 * If *type* is set to *full* then each input value is correlated with
   itself and every other value (i.e., :math:`C_{ij} = V_i V_j` for
   :math:`\{i,j\} = \{1,N\}`, so :math:`N_\text{pair} = N^2`).
+* If *type* is set to *first* then each input value is correlated with
+  the first input value (i.e., :math:`C_{ij} = V_1 V_j` for
+  :math:`\{i,j\} = \{1,N\}`, so :math:`N_\text{pair} = N`).
 
 The *ave* keyword determines what happens to the accumulation of correlation
 samples every :math:`N_\text{freq}` timesteps.  If the *ave* setting is *one*,
@@ -369,6 +373,8 @@ above.
 * For *type* = *full*, the :math:`N_\text{pair} = N^2` columns are ordered:
   :math:`C_{11}, C_{12}, \dotsc, C_{1N}, C_{21}, C_{22}, \dotsc, C_{2N},
   C_{31}, \dotsc, C_{3N}, \dotsc, C_{N1}, \dotsc, C_{N,N-1}, C_{NN}`
+* For *type* = *first*, the :math:`N_\text{pair} = N` columns are ordered:
+  :math:`C_{11}, C_{12}, \dotsc, C_{1N}`
 
 The array values calculated by this fix are treated as extensive.  If
 you need to divide them by the number of atoms, you must do this in a
diff --git a/src/fix_ave_correlate.cpp b/src/fix_ave_correlate.cpp
index 10280682df..a48788b777 100644
--- a/src/fix_ave_correlate.cpp
+++ b/src/fix_ave_correlate.cpp
@@ -36,7 +36,7 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 
 enum { ONE, RUNNING };
-enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL };
+enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL, FIRST };
 
 /* ---------------------------------------------------------------------- */
 
@@ -108,6 +108,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER;
       else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER;
       else if (strcmp(arg[iarg+1],"full") == 0) type = FULL;
+      else if (strcmp(arg[iarg+1], "first") == 0) type = FIRST;
       else error->all(FLERR, iarg+1, "Unknown fix ave/correlate type: {}");
       iarg += 2;
     } else if (strcmp(arg[iarg],"ave") == 0) {
@@ -203,7 +204,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) :
 
   // npair = # of correlation pairs to calculate
 
-  if (type == AUTO) npair = nvalues;
+  if (type == AUTO || type == FIRST) npair = nvalues;
   if (type == UPPER || type == LOWER) npair = nvalues*(nvalues-1)/2;
   if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues*(nvalues+1)/2;
   if (type == FULL) npair = nvalues*nvalues;
@@ -242,6 +243,9 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) :
         for (int i = 0; i < nvalues; i++)
           for (int j = 0; j < nvalues; j++)
             fprintf(fp," %s*%s",earg[i],earg[j]);
+      else if (type == FIRST)
+        for (int i = 0; i < nvalues; i++)
+          fprintf(fp," %s*%s",earg[0],earg[i]);
       fprintf(fp,"\n");
     }
     if (ferror(fp))
@@ -579,6 +583,16 @@ void FixAveCorrelate::accumulate()
       m--;
       if (m < 0) m = nrepeat-1;
     }
+  } else if (type == FIRST) {
+    m = n = lastindex;
+    for (k = 0; k < nsample; k++) {
+      ipair = 0;
+      for (i = 0; i < nvalues; i++) {
+        corr[k][ipair++] += cvalues[m][0]*cvalues[n][i];
+      }
+      m--;
+      if (m < 0) m = nrepeat-1;
+    }
   }
 }