diff --git a/doc/src/PDF/pair_gayberne_extra.pdf b/doc/src/PDF/pair_gayberne_extra.pdf
index b794ac3382..5ac8bd051b 100644
Binary files a/doc/src/PDF/pair_gayberne_extra.pdf and b/doc/src/PDF/pair_gayberne_extra.pdf differ
diff --git a/doc/src/PDF/pair_gayberne_extra.tex b/doc/src/PDF/pair_gayberne_extra.tex
index a63ab2b1ce..0568b54dcf 100644
--- a/doc/src/PDF/pair_gayberne_extra.tex
+++ b/doc/src/PDF/pair_gayberne_extra.tex
@@ -131,7 +131,7 @@ and
 
 $$ \frac{ \partial \chi_{12} }{ \partial \mathbf{q}_i } = 4.0 \cdot
 r^{-2} \cdot \mathbf{A}_i (- \mathbf{\iota}^T \cdot \mathbf{B}_i
-\times \mathbf{\iota} ). $$
+\times \mathbf{\iota} ) \cdot \mu \cdot \chi_{12}^{ ( \mu -1 ) / \mu}. $$
 
 For the derivative of the $\eta$ term, we were unable to find a matrix
 expression due to the determinant. Let $a_{mi}$ be the mth row of the
diff --git a/lib/gpu/lal_gayberne.cu b/lib/gpu/lal_gayberne.cu
index dc6e00ec82..cd1ee59fc6 100644
--- a/lib/gpu/lal_gayberne.cu
+++ b/lib/gpu/lal_gayberne.cu
@@ -316,10 +316,9 @@ __kernel void k_gayberne(const __global numtyp4 *restrict x_,
         numtyp tempv[3];
         gpu_row_times3(iota,b1,tempv);
         gpu_cross3(tempv,iota,tchi);
-        temp1 = (numtyp)-4.0*ir*ir;
-        tchi[0] *= temp1;
-        tchi[1] *= temp1;
-        tchi[2] *= temp1;
+        tchi[0] *= temp2;
+        tchi[1] *= temp2;
+        tchi[2] *= temp2;
       }
 
       numtyp temp2 = factor_lj*eta*chi;
diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp
index ea5a7e9aaa..7dd5270309 100644
--- a/src/ASPHERE/pair_gayberne.cpp
+++ b/src/ASPHERE/pair_gayberne.cpp
@@ -28,7 +28,6 @@
 #include "citeme.h"
 #include "memory.h"
 #include "error.h"
-#include "utils.h"
 
 using namespace LAMMPS_NS;
 
@@ -461,20 +460,20 @@ void PairGayBerne::read_restart(FILE *fp)
   int i,j;
   int me = comm->me;
   for (i = 1; i <= atom->ntypes; i++) {
-    if (me == 0) utils::sfread(FLERR,&setwell[i],sizeof(int),1,fp,NULL,error);
+    if (me == 0) fread(&setwell[i],sizeof(int),1,fp);
     MPI_Bcast(&setwell[i],1,MPI_INT,0,world);
     if (setwell[i]) {
-      if (me == 0) utils::sfread(FLERR,&well[i][0],sizeof(double),3,fp,NULL,error);
+      if (me == 0) fread(&well[i][0],sizeof(double),3,fp);
       MPI_Bcast(&well[i][0],3,MPI_DOUBLE,0,world);
     }
     for (j = i; j <= atom->ntypes; j++) {
-      if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error);
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
       if (setflag[i][j]) {
         if (me == 0) {
-          utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,NULL,error);
-          utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,NULL,error);
-          utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,NULL,error);
+          fread(&epsilon[i][j],sizeof(double),1,fp);
+          fread(&sigma[i][j],sizeof(double),1,fp);
+          fread(&cut[i][j],sizeof(double),1,fp);
         }
         MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
         MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
@@ -506,12 +505,12 @@ void PairGayBerne::read_restart_settings(FILE *fp)
 {
   int me = comm->me;
   if (me == 0) {
-    utils::sfread(FLERR,&gamma,sizeof(double),1,fp,NULL,error);
-    utils::sfread(FLERR,&upsilon,sizeof(double),1,fp,NULL,error);
-    utils::sfread(FLERR,&mu,sizeof(double),1,fp,NULL,error);
-    utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,NULL,error);
-    utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,NULL,error);
-    utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,NULL,error);
+    fread(&gamma,sizeof(double),1,fp);
+    fread(&upsilon,sizeof(double),1,fp);
+    fread(&mu,sizeof(double),1,fp);
+    fread(&cut_global,sizeof(double),1,fp);
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
   }
   MPI_Bcast(&gamma,1,MPI_DOUBLE,0,world);
   MPI_Bcast(&upsilon,1,MPI_DOUBLE,0,world);
@@ -644,10 +643,10 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
   dchi[2] = temp2*(iota[2]-temp1*r12hat[2]);
 
   temp1 = -eta*u_r;
-  temp2 = eta*chi;
-  fforce[0] = temp1*dchi[0]-temp2*dUr[0];
-  fforce[1] = temp1*dchi[1]-temp2*dUr[1];
-  fforce[2] = temp1*dchi[2]-temp2*dUr[2];
+  temp3 = eta*chi;
+  fforce[0] = temp1*dchi[0]-temp3*dUr[0];
+  fforce[1] = temp1*dchi[1]-temp3*dUr[1];
+  fforce[2] = temp1*dchi[2]-temp3*dUr[2];
 
   // torque for particle 1 and 2
   // compute dUr
@@ -668,18 +667,17 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
 
   MathExtra::vecmat(iota,b1,tempv);
   MathExtra::cross3(tempv,iota,dchi);
-  temp1 = -4.0/rsq;
-  dchi[0] *= temp1;
-  dchi[1] *= temp1;
-  dchi[2] *= temp1;
+  dchi[0] *= temp2;
+  dchi[1] *= temp2;
+  dchi[2] *= temp2;
   double dchi2[3];
 
   if (newton_pair || j < nlocal) {
     MathExtra::vecmat(iota,b2,tempv);
     MathExtra::cross3(tempv,iota,dchi2);
-    dchi2[0] *= temp1;
-    dchi2[1] *= temp1;
-    dchi2[2] *= temp1;
+    dchi2[0] *= temp2;
+    dchi2[1] *= temp2;
+    dchi2[2] *= temp2;
   }
 
   // compute d_eta
diff --git a/src/USER-INTEL/pair_gayberne_intel.cpp b/src/USER-INTEL/pair_gayberne_intel.cpp
index 862dee2287..30e61f67b9 100644
--- a/src/USER-INTEL/pair_gayberne_intel.cpp
+++ b/src/USER-INTEL/pair_gayberne_intel.cpp
@@ -555,10 +555,10 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,
           dchi_2 = temp2 * (iota_2 - temp1 * r12hat_2);
 
           temp1 = -eta * u_r;
-          temp2 = eta * chi;
-          fforce_0 = temp1 * dchi_0 - temp2 * dUr_0;
-          fforce_1 = temp1 * dchi_1 - temp2 * dUr_1;
-          fforce_2 = temp1 * dchi_2 - temp2 * dUr_2;
+          temp3 = eta * chi;
+          fforce_0 = temp1 * dchi_0 - temp3 * dUr_0;
+          fforce_1 = temp1 * dchi_1 - temp3 * dUr_1;
+          fforce_2 = temp1 * dchi_2 - temp3 * dUr_2;
 
           // torque for particle 1 and 2
           // compute dUr
@@ -579,18 +579,17 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,
 
           ME_vecmat(iota, b1, tempv);
           ME_cross3(tempv, iota, dchi);
-          temp1 = (flt_t)-4.0 / rsq_form[jj];
-          dchi_0 *= temp1;
-          dchi_1 *= temp1;
-          dchi_2 *= temp1;
+          dchi_0 *= temp2;
+          dchi_1 *= temp2;
+          dchi_2 *= temp2;
           flt_t dchi2_0, dchi2_1, dchi2_2;
 
           if (NEWTON_PAIR) {
             ME_vecmat(iota, b2, tempv);
             ME_cross3(tempv, iota, dchi2);
-            dchi2_0 *= temp1;
-            dchi2_1 *= temp1;
-            dchi2_2 *= temp1;
+            dchi2_0 *= temp2;
+            dchi2_1 *= temp2;
+            dchi2_2 *= temp2;
           }
 
           // compute d_eta