diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp
index 034cfc4c1..b288c4702 100644
--- a/src/MOLECULE/angle_cosine_periodic.cpp
+++ b/src/MOLECULE/angle_cosine_periodic.cpp
@@ -1,301 +1,303 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Tod A Pascal (Caltech)
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "angle_cosine_periodic.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
 #include "math_const.h"
 #include "math_special.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 using namespace MathSpecial;
 
 #define SMALL 0.001
 
 /* ---------------------------------------------------------------------- */
 
 AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
 AngleCosinePeriodic::~AngleCosinePeriodic()
 {
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(k);
     memory->destroy(b);
     memory->destroy(multiplicity);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::compute(int eflag, int vflag)
 {
   int i,i1,i2,i3,n,m,type,b_factor;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double eangle,f1[3],f3[3];
   double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
   double tn,tn_1,tn_2,un,un_1,un_2;
 
   eangle = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = 0;
 
   double **x = atom->x;
   double **f = atom->f;
   int **anglelist = neighbor->anglelist;
   int nanglelist = neighbor->nanglelist;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   for (n = 0; n < nanglelist; n++) {
     i1 = anglelist[n][0];
     i2 = anglelist[n][1];
     i3 = anglelist[n][2];
     type = anglelist[n][3];
 
     // 1st bond
 
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
 
     // 2nd bond
 
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
 
     // c = cosine of angle
 
     c = delx1*delx2 + dely1*dely2 + delz1*delz2;
     c /= r1*r2;
     if (c > 1.0) c = 1.0;
     if (c < -1.0) c = -1.0;
 
     m = multiplicity[type];
     b_factor = b[type];
 
     // cos(n*x) = Tn(cos(x))
     // Tn(x) = Chebyshev polynomials of the first kind: T_0 = 1, T_1 = x, ...
     // recurrence relationship:
     // Tn(x) = 2*x*T[n-1](x) - T[n-2](x) where T[-1](x) = 0
     // also, dTn(x)/dx = n*U[n-1](x)
     // where Un(x) = 2*x*U[n-1](x) - U[n-2](x) and U[-1](x) = 0
     // finally need to handle special case for n = 1
 
     tn = 1.0;
     tn_1 = 1.0;
     tn_2 = 0.0;
     un = 1.0;
     un_1 = 2.0;
     un_2 = 0.0;
 
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;
     s = 1.0/s;
 
     // force & energy
 
     tn_2 = c;
     for (i = 1; i <= m; i++) {
       tn = 2*c*tn_1 - tn_2;
       tn_2 = tn_1;
       tn_1 = tn;
     }
 
     for (i = 2; i <= m; i++) {
       un = 2*c*un_1 - un_2;
       un_2 = un_1;
       un_1 = un;
     }
     tn = b_factor*powsign(m)*tn;
     un = b_factor*powsign(m)*m*un;
 
     if (eflag) eangle = 2*k[type]*(1.0 - tn);
 
     a = -k[type]*un;
     a11 = a*c / rsq1;
     a12 = -a / (r1*r2);
     a22 = a*c / rsq2;
 
     f1[0] = a11*delx1 + a12*delx2;
     f1[1] = a11*dely1 + a12*dely2;
     f1[2] = a11*delz1 + a12*delz2;
     f3[0] = a22*delx2 + a12*delx1;
     f3[1] = a22*dely2 + a12*dely1;
     f3[2] = a22*delz2 + a12*delz1;
 
     // apply force to each of 3 atoms
 
     if (newton_bond || i1 < nlocal) {
       f[i1][0] += f1[0];
       f[i1][1] += f1[1];
       f[i1][2] += f1[2];
     }
 
     if (newton_bond || i2 < nlocal) {
       f[i2][0] -= f1[0] + f3[0];
       f[i2][1] -= f1[1] + f3[1];
       f[i2][2] -= f1[2] + f3[2];
     }
 
     if (newton_bond || i3 < nlocal) {
       f[i3][0] += f3[0];
       f[i3][1] += f3[1];
       f[i3][2] += f3[2];
     }
 
     if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
                          delx1,dely1,delz1,delx2,dely2,delz2);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::allocate()
 {
   allocated = 1;
   int n = atom->nangletypes;
 
   memory->create(k,n+1,"angle:k");
   memory->create(multiplicity,n+1,"angle:multiplicity");
   memory->create(b,n+1,"angle:b");
 
   memory->create(setflag,n+1,"angle:setflag");
   for (int i = 1; i <= n; i++) setflag[i] = 0;
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::coeff(int narg, char **arg)
 {
   if (narg != 4) error->all(FLERR,"Incorrect args for angle coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi;
   force->bounds(arg[0],atom->nangletypes,ilo,ihi);
 
   double c_one = force->numeric(FLERR,arg[1]);
   int b_one = force->inumeric(FLERR,arg[2]);
   int n_one = force->inumeric(FLERR,arg[3]);
   if (n_one <= 0) error->all(FLERR,"Incorrect args for angle coefficients");
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = c_one/(n_one*n_one);
     b[i] = b_one;
     multiplicity[i] = n_one;
     setflag[i] = 1;
     count++;
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
 }
 
 /* ---------------------------------------------------------------------- */
 
 double AngleCosinePeriodic::equilibrium_angle(int i)
 {
   return MY_PI;
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes out coeffs to restart file
 ------------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::write_restart(FILE *fp)
 {
   fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&b[1],sizeof(int),atom->nangletypes,fp);
   fwrite(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads coeffs from restart file, bcasts them
 ------------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::read_restart(FILE *fp)
 {
   allocate();
 
   if (comm->me == 0) {
     fread(&k[1],sizeof(double),atom->nangletypes,fp);
     fread(&b[1],sizeof(int),atom->nangletypes,fp);
     fread(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
   }
 
   MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&b[1],atom->nangletypes,MPI_INT,0,world);
   MPI_Bcast(&multiplicity[1],atom->nangletypes,MPI_INT,0,world);
   for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
 }
 
 
 /* ----------------------------------------------------------------------
    proc 0 writes to data file
 ------------------------------------------------------------------------- */
 
 void AngleCosinePeriodic::write_data(FILE *fp)
 {
-  for (int i = 1; i <= atom->nangletypes; i++)
-    fprintf(fp,"%d %g %d %d\n",i,k[i],b[i],multiplicity[i]);
+  for (int i = 1; i <= atom->nangletypes; i++) {
+    int m = multiplicity[i];
+    fprintf(fp,"%d %g %d %d\n",i,k[i]*m*m,b[i],m);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 double AngleCosinePeriodic::single(int type, int i1, int i2, int i3)
 {
   double **x = atom->x;
 
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
   domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
 
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
   domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
   c /= r1*r2;
   if (c > 1.0) c = 1.0;
   if (c < -1.0) c = -1.0;
 
   c = cos(acos(c)*multiplicity[type]);
   return k[type]*(1.0-b[type]*powsign(multiplicity[type])*c);
 }
diff --git a/src/USER-MISC/angle_quartic.cpp b/src/USER-MISC/angle_quartic.cpp
index 8c007143f..c6ba466ef 100644
--- a/src/USER-MISC/angle_quartic.cpp
+++ b/src/USER-MISC/angle_quartic.cpp
@@ -1,289 +1,289 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under 
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Loukas D. Peristeras (Scienomics SARL)
    [ based on angle_harmonic.cpp]
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "angle_quartic.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
 #define SMALL 0.001
 
 /* ---------------------------------------------------------------------- */
 
 AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
 AngleQuartic::~AngleQuartic()
 {
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(k2);
     memory->destroy(k3);
     memory->destroy(k4);
     memory->destroy(theta0);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AngleQuartic::compute(int eflag, int vflag)
 {
   int i1,i2,i3,n,type;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double eangle,f1[3],f3[3];
   double dtheta,dtheta2,dtheta3,dtheta4,tk;
   double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
 
   eangle = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = 0;
 
   double **x = atom->x;
   double **f = atom->f;
   int **anglelist = neighbor->anglelist;
   int nanglelist = neighbor->nanglelist;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   for (n = 0; n < nanglelist; n++) {
     i1 = anglelist[n][0];
     i2 = anglelist[n][1];
     i3 = anglelist[n][2];
     type = anglelist[n][3];
 
     // 1st bond
 
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
 
     // 2nd bond
 
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
 
     // angle (cos and sin)
 
     c = delx1*delx2 + dely1*dely2 + delz1*delz2;
     c /= r1*r2;
 
     if (c > 1.0) c = 1.0;
     if (c < -1.0) c = -1.0;
 
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;
     s = 1.0/s;
 
     // force & energy
 
     dtheta = acos(c) - theta0[type];
     dtheta2 = dtheta * dtheta;
     dtheta3 = dtheta2 * dtheta;
     tk =  2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 
       4.0 * k4[type] * dtheta3;
 
     if (eflag) {
       dtheta4 = dtheta3 * dtheta;
       eangle = k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4;
     }
 
-    a = -2.0 * tk * s;
+    a = -tk * s;
     a11 = a*c / rsq1;
     a12 = -a / (r1*r2);
     a22 = a*c / rsq2;
 
     f1[0] = a11*delx1 + a12*delx2;
     f1[1] = a11*dely1 + a12*dely2;
     f1[2] = a11*delz1 + a12*delz2;
     f3[0] = a22*delx2 + a12*delx1;
     f3[1] = a22*dely2 + a12*dely1;
     f3[2] = a22*delz2 + a12*delz1;
 
     // apply force to each of 3 atoms
 
     if (newton_bond || i1 < nlocal) {
       f[i1][0] += f1[0];
       f[i1][1] += f1[1];
       f[i1][2] += f1[2];
     }
 
     if (newton_bond || i2 < nlocal) {
       f[i2][0] -= f1[0] + f3[0];
       f[i2][1] -= f1[1] + f3[1];
       f[i2][2] -= f1[2] + f3[2];
     }
 
     if (newton_bond || i3 < nlocal) {
       f[i3][0] += f3[0];
       f[i3][1] += f3[1];
       f[i3][2] += f3[2];
     }
 
     if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
 			 delx1,dely1,delz1,delx2,dely2,delz2);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AngleQuartic::allocate()
 {
   allocated = 1;
   int n = atom->nangletypes;
 
   memory->create(k2,n+1,"angle:k2");
   memory->create(k3,n+1,"angle:k3");
   memory->create(k4,n+1,"angle:k4");
   memory->create(theta0,n+1,"angle:theta0");
 
   memory->create(setflag,n+1,"angle:setflag");
   for (int i = 1; i <= n; i++) setflag[i] = 0;
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
 
 void AngleQuartic::coeff(int narg, char **arg)
 {
   if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi;
   force->bounds(arg[0],atom->nangletypes,ilo,ihi);
 
   double theta0_one = force->numeric(FLERR,arg[1]);
   double k2_one = force->numeric(FLERR,arg[2]);
   double k3_one = force->numeric(FLERR,arg[3]);
   double k4_one = force->numeric(FLERR,arg[4]);
 
   // convert theta0 from degrees to radians
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k2[i] = k2_one;
     k3[i] = k3_one;
     k4[i] = k4_one;
     theta0[i] = theta0_one/180.0 * MY_PI;
     setflag[i] = 1;
     count++;
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
 }
 
 /* ---------------------------------------------------------------------- */
 
 double AngleQuartic::equilibrium_angle(int i)
 {
   return theta0[i];
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes out coeffs to restart file
 ------------------------------------------------------------------------- */
 
 void AngleQuartic::write_restart(FILE *fp)
 {
   fwrite(&k2[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&k3[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&k4[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads coeffs from restart file, bcasts them 
 ------------------------------------------------------------------------- */
 
 void AngleQuartic::read_restart(FILE *fp)
 {
   allocate();
 
   if (comm->me == 0) {
     fread(&k2[1],sizeof(double),atom->nangletypes,fp);
     fread(&k3[1],sizeof(double),atom->nangletypes,fp);
     fread(&k4[1],sizeof(double),atom->nangletypes,fp);
     fread(&theta0[1],sizeof(double),atom->nangletypes,fp);
   }
   MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&k3[1],atom->nangletypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&k4[1],atom->nangletypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
 
   for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes to data file
 ------------------------------------------------------------------------- */
 
 void AngleQuartic::write_data(FILE *fp)
 {
   for (int i = 1; i <= atom->nangletypes; i++)
-    fprintf(fp,"%d %g %g %g %g\n",i,k2[i],k3[i],k4[i],theta0[i]/MY_PI*180.0);
+    fprintf(fp,"%d %g %g %g %g\n",i,theta0[i]/MY_PI*180.0,k2[i],k3[i],k4[i]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 double AngleQuartic::single(int type, int i1, int i2, int i3)
 {
   double **x = atom->x;
 
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
   domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
 
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
   domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
   c /= r1*r2;
   if (c > 1.0) c = 1.0;
   if (c < -1.0) c = -1.0;
 
   double dtheta = acos(c) - theta0[type];
   double dtheta2 = dtheta * dtheta;
   double dtheta3 = dtheta2 * dtheta;
   double dtheta4 = dtheta3 * dtheta;
   double tk = 2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 
     4.0 * k4[type] * dtheta3;
   return k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4;
 }
diff --git a/src/USER-MISC/bond_harmonic_shift.cpp b/src/USER-MISC/bond_harmonic_shift.cpp
index e8c2b7d78..3b809a859 100644
--- a/src/USER-MISC/bond_harmonic_shift.cpp
+++ b/src/USER-MISC/bond_harmonic_shift.cpp
@@ -1,201 +1,213 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Carsten Svaneborg, science@zqex.dk
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "bond_harmonic_shift.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 BondHarmonicShift::BondHarmonicShift(LAMMPS *lmp) : Bond(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
 BondHarmonicShift::~BondHarmonicShift()
 {
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(k);
     memory->destroy(r0);
     memory->destroy(r1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void BondHarmonicShift::compute(int eflag, int vflag)
 {
   int i1,i2,n,type;
   double delx,dely,delz,ebond,fbond;
   double rsq,r,dr,rk;
 
   ebond = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = 0;
 
   double **x = atom->x;
   double **f = atom->f;
   int **bondlist = neighbor->bondlist;
   int nbondlist = neighbor->nbondlist;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   for (n = 0; n < nbondlist; n++) {
     i1 = bondlist[n][0];
     i2 = bondlist[n][1];
     type = bondlist[n][2];
 
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
 
     rsq = delx*delx + dely*dely + delz*delz;
     r = sqrt(rsq);
 
     dr = r - r0[type];
     rk = k[type] * dr;
 
     // force & energy
 
     if (r > 0.0) fbond = -2.0*rk/r;
     else fbond = 0.0;
 
     if (eflag)
       ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]) );
 
     // apply force to each of 2 atoms
 
     if (newton_bond || i1 < nlocal) {
       f[i1][0] += delx*fbond;
       f[i1][1] += dely*fbond;
       f[i1][2] += delz*fbond;
     }
 
     if (newton_bond || i2 < nlocal) {
       f[i2][0] -= delx*fbond;
       f[i2][1] -= dely*fbond;
       f[i2][2] -= delz*fbond;
     }
 
     if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void BondHarmonicShift::allocate()
 {
   allocated = 1;
   int n = atom->nbondtypes;
 
   memory->create(k ,    n+1,"bond:k");
   memory->create(r0,    n+1,"bond:r0");
   memory->create(r1,    n+1,"bond:r1");
   memory->create(setflag,n+1,"bond:setflag");
 
   for (int i = 1; i <= n; i++) setflag[i] = 0;
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShift::coeff(int narg, char **arg)
 {
   if (narg != 4) error->all(FLERR,"Incorrect args for bond coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi;
   force->bounds(arg[0],atom->nbondtypes,ilo,ihi);
 
   double Umin = force->numeric(FLERR,arg[1]);   // energy at minimum
   double r0_one = force->numeric(FLERR,arg[2]); // position of minimum
   double r1_one = force->numeric(FLERR,arg[3]);  // position where energy = 0
+  if (r0_one == r1_one)
+    error->all(FLERR,"Bond harmonic/shift r0 and r1 must be different");
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = Umin/((r0_one-r1_one)*(r0_one-r1_one));
     r0[i] = r0_one;
     r1[i] = r1_one;
     setflag[i] = 1;
     count++;
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
 }
 
 /* ----------------------------------------------------------------------
    return an equilbrium bond length
 ------------------------------------------------------------------------- */
 
 double BondHarmonicShift::equilibrium_distance(int i)
 {
   return r0[i];
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes out coeffs to restart file
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShift::write_restart(FILE *fp)
 {
   fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);
   fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
   fwrite(&r1[1],sizeof(double),atom->nbondtypes,fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads coeffs from restart file, bcasts them
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShift::read_restart(FILE *fp)
 {
   allocate();
 
   if (comm->me == 0) {
     fread(&k[1],sizeof(double),atom->nbondtypes,fp);
     fread(&r0[1],sizeof(double),atom->nbondtypes,fp);
     fread(&r1[1],sizeof(double),atom->nbondtypes,fp);
   }
   MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&r1[1],atom->nbondtypes,MPI_DOUBLE,0,world);
 
   for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
 }
 
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void BondHarmonicShift::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+    fprintf(fp,"%d %g %g %g\n",i,k[i],r0[i],r1[i]);
+}
+
 /* ---------------------------------------------------------------------- */
 
 double BondHarmonicShift::single(int type, double rsq, int i, int j,
 				 double &fforce)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double dr2=r0[type]-r1[type];
 
   fforce =  -2.0*k[type]*dr/r;
   return k[type]*(dr*dr - dr2*dr2);
 }
diff --git a/src/USER-MISC/bond_harmonic_shift.h b/src/USER-MISC/bond_harmonic_shift.h
index 0a5751cb0..a9aabf96e 100644
--- a/src/USER-MISC/bond_harmonic_shift.h
+++ b/src/USER-MISC/bond_harmonic_shift.h
@@ -1,48 +1,49 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef BOND_CLASS
 
 BondStyle(harmonic/shift,BondHarmonicShift)
 
 #else
 
 #ifndef LMP_BOND_HARMONIC_SHIFT_H
 #define LMP_BOND_HARMONIC_SHIFT_H
 
 #include "stdio.h"
 #include "bond.h"
 
 namespace LAMMPS_NS {
 
 class BondHarmonicShift : public Bond {
  public:
   BondHarmonicShift(class LAMMPS *);
   virtual ~BondHarmonicShift();
   virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
+  void write_data(FILE *);
   double single(int, double, int, int, double &);
 
  protected:
   double *k,*r0,*r1;
 
   void allocate();
 };
 
 }
 
 #endif
 #endif
diff --git a/src/USER-MISC/bond_harmonic_shift_cut.cpp b/src/USER-MISC/bond_harmonic_shift_cut.cpp
index a6d9425cb..f75a50ceb 100644
--- a/src/USER-MISC/bond_harmonic_shift_cut.cpp
+++ b/src/USER-MISC/bond_harmonic_shift_cut.cpp
@@ -1,206 +1,218 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Carsten Svaneborg, science@zqex.dk
 ------------------------------------------------------------------------- */
 
 #include "math.h"
 #include "stdlib.h"
 #include "bond_harmonic_shift_cut.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "domain.h"
 #include "comm.h"
 #include "force.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 BondHarmonicShiftCut::BondHarmonicShiftCut(LAMMPS *lmp) : Bond(lmp) {}
 
 /* ---------------------------------------------------------------------- */
 
 BondHarmonicShiftCut::~BondHarmonicShiftCut()
 {
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(k);
     memory->destroy(r0);
     memory->destroy(r1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void BondHarmonicShiftCut::compute(int eflag, int vflag)
 {
   int i1,i2,n,type;
   double delx,dely,delz,ebond,fbond;
   double rsq,r,dr,rk;
 
   ebond = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = 0;
 
   double **x = atom->x;
   double **f = atom->f;
   int **bondlist = neighbor->bondlist;
   int nbondlist = neighbor->nbondlist;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
   for (n = 0; n < nbondlist; n++) {
     i1 = bondlist[n][0];
     i2 = bondlist[n][1];
     type = bondlist[n][2];
 
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
 
     rsq = delx*delx + dely*dely + delz*delz;
     r = sqrt(rsq);
 
     if (r>r1[type]) continue;
 
     dr = r - r0[type];
     rk = k[type] * dr;
 
     // force & energy
 
     if (r > 0.0) fbond = -2.0*rk/r;
     else fbond = 0.0;
 
     if (eflag)
       ebond = k[type]*(dr*dr -(r0[type]-r1[type])*(r0[type]-r1[type]));
 
     // apply force to each of 2 atoms
 
     if (newton_bond || i1 < nlocal) {
       f[i1][0] += delx*fbond;
       f[i1][1] += dely*fbond;
       f[i1][2] += delz*fbond;
     }
 
     if (newton_bond || i2 < nlocal) {
       f[i2][0] -= delx*fbond;
       f[i2][1] -= dely*fbond;
       f[i2][2] -= delz*fbond;
     }
 
     if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void BondHarmonicShiftCut::allocate()
 {
   allocated = 1;
   int n = atom->nbondtypes;
 
   memory->create(k ,    n+1,"bond:k");
   memory->create(r0,    n+1,"bond:r0");
   memory->create(r1,    n+1,"bond:r1");
   memory->create(setflag,n+1,"bond:setflag");
   for (int i = 1; i <= n; i++) setflag[i] = 0;
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShiftCut::coeff(int narg, char **arg)
 {
   if (narg != 4) error->all(FLERR,"Incorrect args for bond coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi;
   force->bounds(arg[0],atom->nbondtypes,ilo,ihi);
 
   double Umin = force->numeric(FLERR,arg[1]);   // energy at minimum
   double r0_one = force->numeric(FLERR,arg[2]); // position of minimum
   double r1_one = force->numeric(FLERR,arg[3]);  // position where energy = 0 = cutoff
+  if (r0_one == r1_one)
+    error->all(FLERR,"Bond harmonic/shift/cut r0 and r1 must be different");
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     k[i] = Umin/((r0_one-r1_one)*(r0_one-r1_one));
     r0[i] = r0_one;
     r1[i] = r1_one;
     setflag[i] = 1;
     count++;
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
 }
 
 /* ----------------------------------------------------------------------
    return an equilbrium bond length
 ------------------------------------------------------------------------- */
 
 double BondHarmonicShiftCut::equilibrium_distance(int i)
 {
   return r0[i];
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes out coeffs to restart file
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShiftCut::write_restart(FILE *fp)
 {
   fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);
   fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
   fwrite(&r1[1],sizeof(double),atom->nbondtypes,fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads coeffs from restart file, bcasts them
 ------------------------------------------------------------------------- */
 
 void BondHarmonicShiftCut::read_restart(FILE *fp)
 {
   allocate();
 
   if (comm->me == 0) {
     fread(&k[1],sizeof(double),atom->nbondtypes,fp);
     fread(&r0[1],sizeof(double),atom->nbondtypes,fp);
     fread(&r1[1],sizeof(double),atom->nbondtypes,fp);
   }
   MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&r1[1],atom->nbondtypes,MPI_DOUBLE,0,world);
 
   for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
 }
 
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void BondHarmonicShiftCut::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+    fprintf(fp,"%d %g %g %g\n",i,k[i],r0[i],r1[i]);
+}
+
 /* ---------------------------------------------------------------------- */
 
 double BondHarmonicShiftCut::single(int type, double rsq, int i, int j,
 				    double &fforce)
 {
   fforce = 0.0;
   double r = sqrt(rsq);
 
   if (r>r1[type]) return 0.0;
 
   double dr = r - r0[type];
   double dr2=r0[type]-r1[type];
 
   fforce = -2.0*k[type]*dr/r;
   return k[type]*(dr*dr - dr2*dr2);
 }
diff --git a/src/USER-MISC/bond_harmonic_shift_cut.h b/src/USER-MISC/bond_harmonic_shift_cut.h
index 5059b75ea..e74364160 100644
--- a/src/USER-MISC/bond_harmonic_shift_cut.h
+++ b/src/USER-MISC/bond_harmonic_shift_cut.h
@@ -1,48 +1,49 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #ifdef BOND_CLASS
 
 BondStyle(harmonic/shift/cut,BondHarmonicShiftCut)
 
 #else
 
 #ifndef LMP_BOND_HARMONIC_SHIFT_CUT_H
 #define LMP_BOND_HARMONIC_SHIFT_CUT_H
 
 #include "stdio.h"
 #include "bond.h"
 
 namespace LAMMPS_NS {
 
 class BondHarmonicShiftCut : public Bond {
  public:
   BondHarmonicShiftCut(class LAMMPS *);
   virtual ~BondHarmonicShiftCut();
   virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
+  void write_data(FILE *);
   double single(int, double, int, int, double &);
 
  protected:
   double *k,*r0,*r1;
 
   void allocate();
 };
 
 }
 
 #endif
 #endif
diff --git a/src/USER-OMP/angle_quartic_omp.cpp b/src/USER-OMP/angle_quartic_omp.cpp
index d2e8b1a71..a8119f51b 100644
--- a/src/USER-OMP/angle_quartic_omp.cpp
+++ b/src/USER-OMP/angle_quartic_omp.cpp
@@ -1,181 +1,181 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing author: Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
 #include "angle_quartic_omp.h"
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 #include "domain.h"
 
 #include "math_const.h"
 
 #include <math.h>
 
 #include "suffix.h"
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
 #define SMALL 0.001
 
 /* ---------------------------------------------------------------------- */
 
 AngleQuarticOMP::AngleQuarticOMP(class LAMMPS *lmp)
   : AngleQuartic(lmp), ThrOMP(lmp,THR_ANGLE)
 {
   suffix_flag |= Suffix::OMP;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void AngleQuarticOMP::compute(int eflag, int vflag)
 {
 
   if (eflag || vflag) {
     ev_setup(eflag,vflag);
   } else evflag = 0;
 
   const int nall = atom->nlocal + atom->nghost;
   const int nthreads = comm->nthreads;
   const int inum = neighbor->nanglelist;
 
 #if defined(_OPENMP)
 #pragma omp parallel default(none) shared(eflag,vflag)
 #endif
   {
     int ifrom, ito, tid;
 
     loop_setup_thr(ifrom, ito, tid, inum, nthreads);
     ThrData *thr = fix->get_thr(tid);
     ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
 
     if (inum > 0) {
       if (evflag) {
         if (eflag) {
           if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr);
           else eval<1,1,0>(ifrom, ito, thr);
         } else {
           if (force->newton_bond) eval<1,0,1>(ifrom, ito, thr);
           else eval<1,0,0>(ifrom, ito, thr);
         }
       } else {
         if (force->newton_bond) eval<0,0,1>(ifrom, ito, thr);
         else eval<0,0,0>(ifrom, ito, thr);
       }
     }
     reduce_thr(this, eflag, vflag, thr);
   } // end of omp parallel region
 }
 
 template <int EVFLAG, int EFLAG, int NEWTON_BOND>
 void AngleQuarticOMP::eval(int nfrom, int nto, ThrData * const thr)
 {
   int i1,i2,i3,n,type;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double eangle,f1[3],f3[3];
   double dtheta,dtheta2,dtheta3,dtheta4,tk;
   double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
 
   const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
   dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
   const int4_t * _noalias const anglelist = (int4_t *) neighbor->anglelist[0];
   const int nlocal = atom->nlocal;
 
   for (n = nfrom; n < nto; n++) {
     i1 = anglelist[n].a;
     i2 = anglelist[n].b;
     i3 = anglelist[n].c;
     type = anglelist[n].t;
 
     // 1st bond
 
     delx1 = x[i1].x - x[i2].x;
     dely1 = x[i1].y - x[i2].y;
     delz1 = x[i1].z - x[i2].z;
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
 
     // 2nd bond
 
     delx2 = x[i3].x - x[i2].x;
     dely2 = x[i3].y - x[i2].y;
     delz2 = x[i3].z - x[i2].z;
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
 
     // angle (cos and sin)
 
     c = delx1*delx2 + dely1*dely2 + delz1*delz2;
     c /= r1*r2;
 
     if (c > 1.0) c = 1.0;
     if (c < -1.0) c = -1.0;
 
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;
     s = 1.0/s;
 
     // force & energy
 
     dtheta = acos(c) - theta0[type];
     dtheta2 = dtheta * dtheta;
     dtheta3 = dtheta2 * dtheta;
     tk =  2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3;
 
     if (EFLAG) {
       dtheta4 = dtheta3 * dtheta;
       eangle = k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4;
     }
 
-    a = -2.0 * tk * s;
+    a = -tk * s;
     a11 = a*c / rsq1;
     a12 = -a / (r1*r2);
     a22 = a*c / rsq2;
 
     f1[0] = a11*delx1 + a12*delx2;
     f1[1] = a11*dely1 + a12*dely2;
     f1[2] = a11*delz1 + a12*delz2;
     f3[0] = a22*delx2 + a12*delx1;
     f3[1] = a22*dely2 + a12*dely1;
     f3[2] = a22*delz2 + a12*delz1;
 
     // apply force to each of 3 atoms
 
     if (NEWTON_BOND || i1 < nlocal) {
       f[i1].x += f1[0];
       f[i1].y += f1[1];
       f[i1].z += f1[2];
     }
 
     if (NEWTON_BOND || i2 < nlocal) {
       f[i2].x -= f1[0] + f3[0];
       f[i2].y -= f1[1] + f3[1];
       f[i2].z -= f1[2] + f3[2];
     }
 
     if (NEWTON_BOND || i3 < nlocal) {
       f[i3].x += f3[0];
       f[i3].y += f3[1];
       f[i3].z += f3[2];
     }
 
     if (EVFLAG) ev_tally_thr(this,i1,i2,i3,nlocal,NEWTON_BOND,eangle,f1,f3,
                              delx1,dely1,delz1,delx2,dely2,delz2,thr);
   }
 }