diff --git a/src/fix.cpp b/src/fix.cpp
index b88f3b70ae..1fb1c43f77 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -46,8 +46,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   style = new char[n];
   strcpy(style,arg[2]);
 
-  restart_global = 0;
-  restart_peratom = 0;
+  restart_global = restart_peratom = restart_file = 0;
   force_reneighbor = 0;
   box_change_size = box_change_shape = box_change_domain = 0;
   thermo_energy = 0;
diff --git a/src/fix.h b/src/fix.h
index c66552dfed..281edc6216 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -25,6 +25,7 @@ class Fix : protected Pointers {
 
   int restart_global;            // 1 if Fix saves global state, 0 if not
   int restart_peratom;           // 1 if Fix saves peratom state, 0 if not
+  int restart_file;              // 1 if Fix writes own restart file, 0 if not
   int force_reneighbor;          // 1 if Fix forces reneighboring, 0 if not
 
   int box_change_size;           // 1 if Fix changes box size, 0 if not
@@ -108,6 +109,7 @@ class Fix : protected Pointers {
   virtual void end_of_step() {}
   virtual void post_run() {}
   virtual void write_restart(FILE *) {}
+  virtual void write_restart_file() {}
   virtual void restart(char *) {}
 
   virtual void grow_arrays(int) {}
diff --git a/src/math_extra.h b/src/math_extra.h
index d1691f395b..8b0af2a460 100755
--- a/src/math_extra.h
+++ b/src/math_extra.h
@@ -41,29 +41,31 @@ namespace MathExtra {
 
   // 3x3 matrix operations
 
+  inline void col2mat(const double *ex, const double *ey, const double *ez,
+                      double m[3][3]);
   inline double det3(const double mat[3][3]);
-  inline void diag_times3(const double *diagonal, const double mat[3][3],
+  inline void diag_times3(const double *d, const double m[3][3],
+                          double ans[3][3]);
+  inline void times3_diag(const double m[3][3], const double *d,
                           double ans[3][3]);
   inline void plus3(const double m[3][3], const double m2[3][3],
                     double ans[3][3]);
   inline void times3(const double m[3][3], const double m2[3][3],
                      double ans[3][3]);
-  inline void transpose_times3(const double mat1[3][3],
-                               const double mat2[3][3],
+  inline void transpose_times3(const double m[3][3], const double m2[3][3],
                                double ans[3][3]);
-  inline void times3_transpose(const double mat1[3][3],
-                               const double mat2[3][3],
+  inline void times3_transpose(const double m[3][3], const double m2[3][3],
                                double ans[3][3]);
   inline void invert3(const double mat[3][3], double ans[3][3]);
   inline void matvec(const double mat[3][3], const double *vec, double *ans);
   inline void matvec(const double *ex, const double *ey, const double *ez,
                      const double *vec, double *ans);
-  inline void transpose_matvec(const double mat[3][3], const double*vec,
+  inline void transpose_matvec(const double mat[3][3], const double *vec,
                                double *ans);
   inline void transpose_matvec(const double *ex, const double *ey,
                                const double *ez, const double *v,
                                double *ans);
-  inline void transpose_diag3(const double mat[3][3], const double*vec,
+  inline void transpose_diag3(const double m[3][3], const double *d,
                               double ans[3][3]);
   inline void vecmat(const double *v, const double m[3][3], double *ans);
   inline void scalar_times3(const double f, double m[3][3]);
@@ -238,6 +240,24 @@ void MathExtra::cross3(const double *v1, const double *v2, double *ans)
   ans[2] = v1[0]*v2[1] - v1[1]*v2[0];
 }
 
+/* ----------------------------------------------------------------------
+   construct matrix from 3 column vectors
+------------------------------------------------------------------------- */
+
+void MathExtra::col2mat(const double *ex, const double *ey, const double *ez,
+                        double m[3][3])
+{
+  m[0][0] = ex[0];
+  m[1][0] = ex[1];
+  m[2][0] = ex[2];
+  m[0][1] = ey[0];
+  m[1][1] = ey[1];
+  m[2][1] = ey[2];
+  m[0][2] = ez[0];
+  m[1][2] = ez[1];
+  m[2][2] = ez[2];
+}
+
 /* ----------------------------------------------------------------------
    determinant of a matrix
 ------------------------------------------------------------------------- */
@@ -268,6 +288,24 @@ void MathExtra::diag_times3(const double *d, const double m[3][3],
   ans[2][2] = d[2]*m[2][2];
 }
 
+/* ----------------------------------------------------------------------
+   full matrix times a diagonal matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::times3_diag(const double m[3][3], const double *d,
+                            double ans[3][3])
+{
+  ans[0][0] = m[0][0]*d[0];
+  ans[0][1] = m[0][1]*d[1];
+  ans[0][2] = m[0][2]*d[2];
+  ans[1][0] = m[1][0]*d[0];
+  ans[1][1] = m[1][1]*d[1];
+  ans[1][2] = m[1][2]*d[2];
+  ans[2][0] = m[2][0]*d[0];
+  ans[2][1] = m[2][1]*d[1];
+  ans[2][2] = m[2][2]*d[2];
+}
+
 /* ----------------------------------------------------------------------
    add two matrices
 ------------------------------------------------------------------------- */
diff --git a/src/write_restart.cpp b/src/write_restart.cpp
index 77f2bed05f..6e1dd883a2 100644
--- a/src/write_restart.cpp
+++ b/src/write_restart.cpp
@@ -29,6 +29,7 @@
 #include "neighbor.h"
 #include "domain.h"
 #include "modify.h"
+#include "fix.h"
 #include "universe.h"
 #include "comm.h"
 #include "output.h"
@@ -292,6 +293,12 @@ void WriteRestart::write(char *file)
   }
 
   memory->destroy(buf);
+
+  // invoke any fixes that write their own restart file
+
+  for (int ifix = 0; ifix < modify->nfix; ifix++)
+    if (modify->fix[ifix]->restart_file)
+      modify->fix[ifix]->write_restart_file();
 }
 
 /* ----------------------------------------------------------------------