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

src/RIGID/fix_rigid_small.cpp

@@ -2510,7 +2510,6 @@ void FixRigidSmall::write_restart_file(char *file)
     buf[i][8] = ispace[0][1];
     buf[i][9] = ispace[0][2];
     buf[i][10] = ispace[1][2];
-    buf[i][10] = ispace[1][2];
     buf[i][11] = body[i].vcm[0];
     buf[i][12] = body[i].vcm[1];
     buf[i][13] = body[i].vcm[2];
@@ -2546,7 +2545,8 @@ void FixRigidSmall::write_restart_file(char *file)
                 "%-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
                 static_cast<int> (buf[i][0]),buf[i][1],
                 buf[i][2],buf[i][3],buf[i][4],
-                buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],buf[i][10],
+                buf[i][5],buf[i][6],buf[i][7],
+                buf[i][8],buf[i][9],buf[i][10],
                 buf[i][11],buf[i][12],buf[i][13],
                 buf[i][14],buf[i][15],buf[i][16],
                 static_cast<int> (buf[i][17]),

src/fix_momentum.cpp

@@ -36,7 +36,7 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) :
   nevery = force->inumeric(FLERR,arg[3]);
   if (nevery <= 0) error->all(FLERR,"Illegal fix momentum command");
 
-  linear = angular = 0;
+  linear = angular = rescale = 0;
 
   int iarg = 4;
   while (iarg < narg) {
@@ -50,6 +50,9 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"angular") == 0) {
       angular = 1;
       iarg += 1;
+    } else if (strcmp(arg[iarg],"rescale") == 0) {
+      rescale = 1;
+      iarg += 1;
     } else error->all(FLERR,"Illegal fix momentum command");
   }
 
@@ -87,6 +90,35 @@ void FixMomentum::init()
 
 void FixMomentum::end_of_step()
 {
+  double **v = atom->v;
+  int *mask = atom->mask;
+  const int nlocal = atom->nlocal;
+  double ekin_old,ekin_new;
+  ekin_old = ekin_new = 0.0;
+
+  // compute kinetic energy before momentum removal, if needed
+
+  if (rescale) {
+
+    double *rmass = atom->rmass;
+    double *mass = atom->mass;
+    int *type = atom->type;
+    double ke=0.0;
+
+    if (rmass) {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit)
+          ke += rmass[i] *
+            (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
+    } else {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit)
+          ke +=  mass[type[i]] *
+            (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
+    }
+    MPI_Allreduce(&ke,&ekin_old,1,MPI_DOUBLE,MPI_SUM,world);
+  }
+
   if (linear) {
     double vcm[3];
     group->vcm(igroup,masstotal,vcm);
@@ -94,10 +126,6 @@ void FixMomentum::end_of_step()
     // adjust velocities by vcm to zero linear momentum
     // only adjust a component if flag is set
 
-    double **v = atom->v;
-    int *mask = atom->mask;
-    int nlocal = atom->nlocal;
-
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
         if (xflag) v[i][0] -= vcm[0];
@@ -137,4 +165,36 @@ void FixMomentum::end_of_step()
         v[i][2] -= omega[0]*dy - omega[1]*dx;
       }
   }
+
+  // compute kinetic energy after momentum removal, if needed
+
+  if (rescale) {
+
+    double ke=0.0, factor=1.0;
+    double *rmass = atom->rmass;
+    double *mass = atom->mass;
+    int *type = atom->type;
+
+    if (rmass) {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit)
+          ke += rmass[i] *
+            (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
+    } else {
+      for (int i = 0; i < nlocal; i++)
+        if (mask[i] & groupbit)
+          ke +=  mass[type[i]] *
+            (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
+    }
+    MPI_Allreduce(&ke,&ekin_new,1,MPI_DOUBLE,MPI_SUM,world);
+
+    if (ekin_new != 0.0) factor = sqrt(ekin_old/ekin_new);
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+        v[i][0] *= factor;
+        v[i][1] *= factor;
+        v[i][2] *= factor;
+      }
+    }
+  }
 }

src/fix_momentum.h

@@ -32,7 +32,7 @@ class FixMomentum : public Fix {
   void end_of_step();
 
  private:
-  int linear,angular;
+  int linear,angular,rescale;
   int xflag,yflag,zflag;
   double masstotal;
 };