git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11054 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/MOLECULE/angle_cosine_periodic.cpp

@@ -269,8 +269,10 @@ void AngleCosinePeriodic::read_restart(FILE *fp)
 
 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);
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-MISC/angle_quartic.cpp

@@ -120,7 +120,7 @@ void AngleQuartic::compute(int eflag, int vflag)
       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;
@@ -253,7 +253,7 @@ void AngleQuartic::read_restart(FILE *fp)
 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]);
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-MISC/bond_harmonic_shift.cpp

@@ -134,6 +134,8 @@ void BondHarmonicShift::coeff(int narg, char **arg)
   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++) {
@@ -187,6 +189,16 @@ void BondHarmonicShift::read_restart(FILE *fp)
   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,

src/USER-MISC/bond_harmonic_shift.h

@@ -34,6 +34,7 @@ class BondHarmonicShift : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
+  void write_data(FILE *);
   double single(int, double, int, int, double &);
 
  protected:

src/USER-MISC/bond_harmonic_shift_cut.cpp

@@ -135,6 +135,8 @@ void BondHarmonicShiftCut::coeff(int narg, char **arg)
   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++) {
@@ -188,6 +190,16 @@ void BondHarmonicShiftCut::read_restart(FILE *fp)
   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,

src/USER-MISC/bond_harmonic_shift_cut.h

@@ -34,6 +34,7 @@ class BondHarmonicShiftCut : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
+  void write_data(FILE *);
   double single(int, double, int, int, double &);
 
  protected:

src/USER-OMP/angle_quartic_omp.cpp

@@ -143,7 +143,7 @@ void AngleQuarticOMP::eval(int nfrom, int nto, ThrData * const thr)
       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;