--- /home/sjplimp/oldlammps/src/fix_spring.cpp	2005-08-04 13:53:37.000000000 -0600
+++ /home/sjplimp/lammps/src/fix_spring.cpp	2005-08-04 11:18:54.000000000 -0600
@@ -15,6 +15,7 @@
    Contributing author: Paul Crozier (SNL)
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "string.h"
 #include "fix_spring.h"
 #include "atom.h"
@@ -25,7 +26,8 @@
 #include "group.h"
 
 #define TETHER 0
-#define COUPLE 1
+#define VECTOR 1
+#define DISTANCE 2
 
 /* ---------------------------------------------------------------------- */
 
@@ -45,13 +47,13 @@
     if (strcmp(arg[7],"NULL") == 0) zflag = 0;
     else zc = atof(arg[7]);
 
-  } else if (strcmp(arg[3],"couple") == 0) {
+  } else if (strcmp(arg[3],"vector") == 0) {
     if (narg != 10) error->all("Illegal fix spring command");
-    styleflag = COUPLE;
+    styleflag = VECTOR;
     group1 = group->find(arg[4]);
     group2 = group->find(arg[5]);
     if (group1 == -1 || group2 == -1) 
-      error->all("Could not find fix spring couple group ID"); 
+      error->all("Could not find fix spring vector group ID"); 
     group1bit = group->bitmask[group1];
     group2bit = group->bitmask[group2];
     k_spring = atof(arg[6]);
@@ -62,6 +64,19 @@
     else yc = atof(arg[8]);
     if (strcmp(arg[9],"NULL") == 0) zflag = 0;
     else zc = atof(arg[9]);
+    
+  } else if (strcmp(arg[3],"distance") == 0) {
+    if (narg != 8) error->all("Illegal fix spring command");
+    styleflag = DISTANCE;
+    group1 = group->find(arg[4]);
+    group2 = group->find(arg[5]);
+    if (group1 == -1 || group2 == -1) 
+      error->all("Could not find fix spring distance group ID"); 
+    group1bit = group->bitmask[group1];
+    group2bit = group->bitmask[group2];
+    k_spring = atof(arg[6]);
+    r0 = atof(arg[7]);
+    if (r0 < 0) error->all("Cannot have r0 < 0 for fix spring distance");
 
   } else error->all("Illegal fix spring command");
 }
@@ -94,7 +109,7 @@
     MPI_Allreduce(&rmass,&masstotal,1,MPI_DOUBLE,MPI_SUM,world);
   }
 
-  if (styleflag == COUPLE) {
+  if (styleflag == VECTOR || styleflag == DISTANCE) {
     double rmass = 0.0;
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit && mask[i] & group1bit) rmass += mass[type[i]];
@@ -127,152 +142,242 @@
 
 void FixSpring::post_force(int vflag)
 {
-  double **x = atom->x;
+  if (styleflag == TETHER) spring_tether();
+  else if (styleflag == VECTOR) spring_vector();
+  else spring_distance();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSpring::spring_tether()
+{
+  double com[3];
+  com_group(com);
+  
+  // fx,fy,fz = components of k * (r-r0)
+  
+  double dx,dy,dz,fx,fy,fz;
+  
+  dx = com[0] - xc;
+  dy = com[1] - yc;
+  dz = com[2] - zc;
+  if (!xflag) dx = 0.0;
+  if (!yflag) dy = 0.0;
+  if (!zflag) dz = 0.0;
+  fx = k_spring * dx;
+  fy = k_spring * dy;
+  fz = k_spring * dz;
+  
+  // apply restoring force to atoms in group
+  // f = -k*(r-r0)*mass/masstotal
+
   double **f = atom->f;
   int *mask = atom->mask;
   int *type = atom->type;
-  int *image = atom->image;
   double *mass = atom->mass;
   int nlocal = atom->nlocal;
+  
+  double massfrac;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      massfrac = mass[type[i]]/masstotal;
+      f[i][0] -= fx*massfrac;
+      f[i][1] -= fy*massfrac;
+      f[i][2] -= fz*massfrac;
+    }
+}
 
-  if (styleflag == COUPLE) {
+/* ---------------------------------------------------------------------- */
+
+void FixSpring::spring_vector()
+{
+  double com[6];
+  com_group2(com);
   
-    // compute center-of-mass of each group
-    // use unwrapped coordinates via image
-    //   else answer is bogus if group straddles a periodic boundary
-
-    double com[6];
-    for (int i = 0; i < 6; i++) com[i] = 0.0;
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
-   
-    int xbox,ybox,zbox;
-    double massone;
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-        if (mask[i] & group1bit) {
-          xbox = image[i]%1000 - 500;
-          ybox = (image[i]/1000) % 1000 - 500;
-          zbox = image[i]/1000000 - 500;
-          massone = mass[type[i]];
-          com[0] += (x[i][0] + xbox*xprd) * massone;
-          com[1] += (x[i][1] + ybox*yprd) * massone;
-          com[2] += (x[i][2] + zbox*zprd) * massone;
-        }
-        if (mask[i] & group2bit) {
-          xbox = image[i]%1000 - 500;
-          ybox = (image[i]/1000) % 1000 - 500;
-          zbox = image[i]/1000000 - 500;
-          massone = mass[type[i]];
-          com[3] += (x[i][0] + xbox*xprd) * massone;
-          com[4] += (x[i][1] + ybox*yprd) * massone;
-          com[5] += (x[i][2] + zbox*zprd) * massone;
-        }
+  // fx,fy,fz = components of k * (r2-r-r0)
+  
+  double dx,dy,dz,fx,fy,fz;
+  
+  dx = com[3] - com[0] - xc;
+  dy = com[4] - com[1] - yc;
+  dz = com[5] - com[2] - zc;
+  if (!xflag) dx = 0.0;
+  if (!yflag) dy = 0.0;
+  if (!zflag) dz = 0.0;
+  fx = k_spring * dx;
+  fy = k_spring * dy;
+  fz = k_spring * dz;
+  
+  // apply restoring force to atoms in group
+  // f = -k*(r-r0)*mass/masstotal
+
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  int nlocal = atom->nlocal;
+
+  double massfrac;
+  for (int i = 0; i < nlocal; i++) {         
+    if (mask[i] & groupbit) {
+      if (mask[i] & group1bit) {
+	massfrac = mass[type[i]]/masstotal1;
+	f[i][0] += fx*massfrac;
+	f[i][1] += fy*massfrac;
+	f[i][2] += fz*massfrac;
+      }
+      if (mask[i] & group2bit) {
+	massfrac = mass[type[i]]/masstotal2;
+	f[i][0] -= fx*massfrac;
+	f[i][1] -= fy*massfrac;
+	f[i][2] -= fz*massfrac;
       }
     }
-   
-    double com_all[6];
-    MPI_Allreduce(com,com_all,6,MPI_DOUBLE,MPI_SUM,world);
-    com_all[0] /= masstotal1;
-    com_all[1] /= masstotal1;
-    com_all[2] /= masstotal1;
-    com_all[3] /= masstotal2;
-    com_all[4] /= masstotal2;
-    com_all[5] /= masstotal2;
-   
-    // fx,fy,fz = components of k * (r2-r-r0)
-   
-    double dx,dy,dz,fx,fy,fz;
-   
-    dx = com_all[3] - com_all[0] - xc;
-    dy = com_all[4] - com_all[1] - yc;
-    dz = com_all[5] - com_all[2] - zc;
-    if (!xflag) dx = 0.0;
-    if (!yflag) dy = 0.0;
-    if (!zflag) dz = 0.0;
-    fx = k_spring * dx;
-    fy = k_spring * dy;
-    fz = k_spring * dz;
-   
-    // apply restoring force to atoms in group
-    // f = -k*(r-r0)*mass/masstotal
-   
-    double massfrac;
-    for (int i = 0; i < nlocal; i++) {         
-      if (mask[i] & groupbit) {
-        if (mask[i] & group1bit) {
-          massfrac = mass[type[i]]/masstotal1;
-          f[i][0] += fx*massfrac;
-          f[i][1] += fy*massfrac;
-          f[i][2] += fz*massfrac;
-        }
-        if (mask[i] & group2bit) {
-          massfrac = mass[type[i]]/masstotal2;
-          f[i][0] -= fx*massfrac;
-          f[i][1] -= fy*massfrac;
-          f[i][2] -= fz*massfrac;
-        }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSpring::spring_distance()
+{
+  double com[6];
+  com_group2(com);
+  
+  // fx,fy,fz = components of k * (r-r0)
+  
+  double dx,dy,dz,fx,fy,fz,r,dr;
+  
+  dx = com[3] - com[0];
+  dy = com[4] - com[1];
+  dz = com[5] - com[2];
+  r = sqrt(dx*dx + dy*dy + dz*dz);
+  dr = r - r0;
+  fx = k_spring*dx*dr/r;
+  fy = k_spring*dy*dr/r;
+  fz = k_spring*dz*dr/r;
+  
+  // apply restoring force to atoms in group
+  // f = -k*(r-r0)*mass/masstotal
+
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  int nlocal = atom->nlocal;
+  
+  double massfrac;
+  for (int i = 0; i < nlocal; i++) {         
+    if (mask[i] & groupbit) {
+      if (mask[i] & group1bit) {
+	massfrac = mass[type[i]]/masstotal1;
+	f[i][0] += fx*massfrac;
+	f[i][1] += fy*massfrac;
+	f[i][2] += fz*massfrac;
+      }
+      if (mask[i] & group2bit) {
+	massfrac = mass[type[i]]/masstotal2;
+	f[i][0] -= fx*massfrac;
+	f[i][1] -= fy*massfrac;
+	f[i][2] -= fz*massfrac;
       }
     }
-    	  
-  } else {
-    
-    // compute center-of-mass of group
-    // use unwrapped coordinates via image
-    //   else answer is bogus if group straddles a periodic boundary
-   
-    double com[3];
-    com[0] = com[1] = com[2] = 0.0;
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
-   
-    int xbox,ybox,zbox;
-    double massone;
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-        xbox = image[i]%1000 - 500;
-        ybox = (image[i]/1000) % 1000 - 500;
-        zbox = image[i]/1000000 - 500;
-        massone = mass[type[i]];
-        com[0] += (x[i][0] + xbox*xprd) * massone;
-        com[1] += (x[i][1] + ybox*yprd) * massone;
-        com[2] += (x[i][2] + zbox*zprd) * massone;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSpring::com_group(double *com_all)
+{
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  int nlocal = atom->nlocal;
+
+  // compute center-of-mass of group
+  // use unwrapped coordinates via image
+  // else answer is bogus if group straddles a periodic boundary
+
+  double com[3];
+  com[0] = com[1] = com[2] = 0.0;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+   
+  int xbox,ybox,zbox;
+  double massone;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      xbox = image[i]%1000 - 500;
+      ybox = (image[i]/1000) % 1000 - 500;
+      zbox = image[i]/1000000 - 500;
+      massone = mass[type[i]];
+      com[0] += (x[i][0] + xbox*xprd) * massone;
+      com[1] += (x[i][1] + ybox*yprd) * massone;
+      com[2] += (x[i][2] + zbox*zprd) * massone;
+    }
+  
+  MPI_Allreduce(com,com_all,3,MPI_DOUBLE,MPI_SUM,world);
+  com_all[0] /= masstotal;
+  com_all[1] /= masstotal;
+  com_all[2] /= masstotal;
+}  
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixSpring::com_group2(double *com_all)
+{
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int *image = atom->image;
+  double *mass = atom->mass;
+  int nlocal = atom->nlocal;
+
+  // compute center-of-mass of each group
+  // use unwrapped coordinates via image
+  // else answer is bogus if group straddles a periodic boundary
+
+  double com[6];
+  for (int i = 0; i < 6; i++) com[i] = 0.0;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  
+  int xbox,ybox,zbox;
+  double massone;
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (mask[i] & group1bit) {
+	xbox = image[i]%1000 - 500;
+	ybox = (image[i]/1000) % 1000 - 500;
+	zbox = image[i]/1000000 - 500;
+	massone = mass[type[i]];
+	com[0] += (x[i][0] + xbox*xprd) * massone;
+	com[1] += (x[i][1] + ybox*yprd) * massone;
+	com[2] += (x[i][2] + zbox*zprd) * massone;
       }
-   
-    double com_all[3];
-    MPI_Allreduce(com,com_all,3,MPI_DOUBLE,MPI_SUM,world);
-    com_all[0] /= masstotal;
-    com_all[1] /= masstotal;
-    com_all[2] /= masstotal;
-   
-    // fx,fy,fz = components of k * (r-r0)
-   
-    double dx,dy,dz,fx,fy,fz;
-   
-    dx = com_all[0] - xc;
-    dy = com_all[1] - yc;
-    dz = com_all[2] - zc;
-    if (!xflag) dx = 0.0;
-    if (!yflag) dy = 0.0;
-    if (!zflag) dz = 0.0;
-    fx = k_spring * dx;
-    fy = k_spring * dy;
-    fz = k_spring * dz;
-   
-    // apply restoring force to atoms in group
-    // f = -k*(r-r0)*mass/masstotal
-   
-    double massfrac;
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-        massfrac = mass[type[i]]/masstotal;
-        f[i][0] -= fx*massfrac;
-        f[i][1] -= fy*massfrac;
-        f[i][2] -= fz*massfrac;
+      if (mask[i] & group2bit) {
+	xbox = image[i]%1000 - 500;
+	ybox = (image[i]/1000) % 1000 - 500;
+	zbox = image[i]/1000000 - 500;
+	massone = mass[type[i]];
+	com[3] += (x[i][0] + xbox*xprd) * massone;
+	com[4] += (x[i][1] + ybox*yprd) * massone;
+	com[5] += (x[i][2] + zbox*zprd) * massone;
       }
+    }
   }
+  
+  MPI_Allreduce(com,com_all,6,MPI_DOUBLE,MPI_SUM,world);
+  com_all[0] /= masstotal1;
+  com_all[1] /= masstotal1;
+  com_all[2] /= masstotal1;
+  com_all[3] /= masstotal2;
+  com_all[4] /= masstotal2;
+  com_all[5] /= masstotal2;
 }
 
 /* ---------------------------------------------------------------------- */
--- /home/sjplimp/oldlammps/src/fix_spring.h	2005-08-04 13:53:37.000000000 -0600
+++ /home/sjplimp/lammps/src/fix_spring.h	2005-08-04 10:48:51.000000000 -0600
@@ -27,13 +27,19 @@
   void post_force_respa(int, int, int);
 
  private:
-  double xc,yc,zc;
+  double xc,yc,zc,r0;
   double k_spring;
   int xflag,yflag,zflag;
   int styleflag;
   int group1,group2,group1bit,group2bit;
   double masstotal,masstotal1,masstotal2;
   int nlevels_respa;
+
+  void spring_tether();
+  void spring_vector();
+  void spring_distance();
+  void com_group(double *);
+  void com_group2(double *);
 };
 
 #endif
--- /home/sjplimp/oldlammps/doc/fix_spring.txt	2005-08-04 13:53:37.000000000 -0600
+++ /home/sjplimp/lammps/doc/fix_spring.txt	2005-08-04 10:24:46.000000000 -0600
@@ -14,20 +14,25 @@
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 spring = style name of this fix command :l
-keyword = {tether} or {couple} :l
+keyword = {tether} or {vector} or {distance}:l
   {tether} values = K x y z
     K = spring constant (force/distance units)
     x,y,z = point to which spring is tethered
-  {couple} values = group-ID1 group-ID2 K x y z
+  {vector} values = group-ID1 group-ID2 K x y z
     group-ID1,group-ID2 = two groups to couple together with a spring
     K = spring constant (force/distance units)
-    x,y,z = equilibrium size and direction of spring :pre
+    x,y,z = equilibrium size and direction of spring
+  {distance} values = group-ID1 group-ID2 K r0
+    group-ID1,group-ID2 = two groups to couple together with a spring
+    K = spring constant (force/distance units)
+    r0 = equilibrium size of spring :pre
 :ule
 
 [Examples:]
 
 fix pull ligand spring tether 50.0 0.0 1.0 1.0
-fix 5 lipids spring couple bilayer1 bilayer2 100.0 NULL NULL 10.0 :pre
+fix 5 lipids spring vector bilayer1 bilayer2 100.0 NULL NULL 10.0
+fix 2 all spring distance atom1 atom2 10.0 3.0 :pre
 
 [Description:]
 
@@ -49,16 +54,19 @@
 represents the total force on the group of atoms, not a per-atom
 force.
 
-The {couple} style links two groups of atoms together with a spring.
+The {vector} style links two groups of atoms together with a spring.
 Only atoms that are in both a coupling group and the fix group are
 treated as coupled.  The specified {x,y,z} determines the orientation
 and distance at which the spring is in equilibrium (exerts no force).
 When the centers of mass of the two groups are not in equilibrium
 positions with respect to each other, an equal and opposite force is
-applied to the atoms of each group to pull them towards an equilibrium
-orientation.
+applied to the atoms of each group to pull them towards equilibrium.
+
+The {distance} style also links two groups of atoms together with a
+spring, but only the distance between the two groups is specified; the
+orientation is not constrained.
 
-For both the {tether} and {couple} styles, any of the x,y,z values can
+For both the {tether} and {vector} styles, any of the x,y,z values can
 be specified as NULL which means do not include that dimension in the
 distance calculation or force application.
 
--- /home/sjplimp/oldlammps/doc/fix_spring.html	2005-08-04 13:53:37.000000000 -0600
+++ /home/sjplimp/lammps/doc/fix_spring.html	2005-08-04 13:53:06.000000000 -0600
@@ -19,22 +19,26 @@
 
 <LI>spring = style name of this fix command 
 
-<LI>keyword = <I>tether</I> or <I>couple</I> 
-
-<PRE>  <I>tether</I> values = K x y z
+<PRE>keyword = <I>tether</I> or <I>vector</I> or <I>distance</I>:l
+  <I>tether</I> values = K x y z
     K = spring constant (force/distance units)
     x,y,z = point to which spring is tethered
-  <I>couple</I> values = group-ID1 group-ID2 K x y z
+  <I>vector</I> values = group-ID1 group-ID2 K x y z
+    group-ID1,group-ID2 = two groups to couple together with a spring
+    K = spring constant (force/distance units)
+    x,y,z = equilibrium size and direction of spring
+  <I>distance</I> values = group-ID1 group-ID2 K r0
     group-ID1,group-ID2 = two groups to couple together with a spring
     K = spring constant (force/distance units)
-    x,y,z = equilibrium size and direction of spring 
+    r0 = equilibrium size of spring 
 </PRE>
 
 </UL>
 <P><B>Examples:</B>
 </P>
 <PRE>fix pull ligand spring tether 50.0 0.0 1.0 1.0
-fix 5 lipids spring couple bilayer1 bilayer2 100.0 NULL NULL 10.0 
+fix 5 lipids spring vector bilayer1 bilayer2 100.0 NULL NULL 10.0
+fix 2 all spring distance atom1 atom2 10.0 3.0 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -56,16 +60,19 @@
 represents the total force on the group of atoms, not a per-atom
 force.
 </P>
-<P>The <I>couple</I> style links two groups of atoms together with a spring.
+<P>The <I>vector</I> style links two groups of atoms together with a spring.
 Only atoms that are in both a coupling group and the fix group are
 treated as coupled.  The specified <I>x,y,z</I> determines the orientation
 and distance at which the spring is in equilibrium (exerts no force).
 When the centers of mass of the two groups are not in equilibrium
 positions with respect to each other, an equal and opposite force is
-applied to the atoms of each group to pull them towards an equilibrium
-orientation.
+applied to the atoms of each group to pull them towards equilibrium.
+</P>
+<P>The <I>distance</I> style also links two groups of atoms together with a
+spring, but only the distance between the two groups is specified; the
+orientation is not constrained.
 </P>
-<P>For both the <I>tether</I> and <I>couple</I> styles, any of the x,y,z values can
+<P>For both the <I>tether</I> and <I>vector</I> styles, any of the x,y,z values can
 be specified as NULL which means do not include that dimension in the
 distance calculation or force application.
 </P>