Added two modes, one via velocity and one via quaternion.

src/USER-MISC/fix_propel_self.cpp

@@ -51,24 +51,34 @@ static constexpr const bool debug_out = false;
 FixPropelSelf::FixPropelSelf( LAMMPS *lmp, int narg, char **argv )
   : Fix(lmp, narg, argv)
 {
-  // if( lmp->citeme) lmp->citeme->add(cite_fix_active);
-  if( narg < 4 ) error->all(FLERR, "Illegal fix propel/self command");
-  
-  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);
-  if (!av) error->all(FLERR, "FixPropelSelf requires atom_style ellipsoid");
+  if( narg < 5 ) error->all(FLERR, "Illegal fix propel/self command");
 
-  if( debug_out && comm->me == 0 ){
-    fprintf(screen, "This is fix active, narg is %d\n", narg );
-    fprintf(screen, "args:");
-    for( int i = 0; i < narg; ++i ){
-      fprintf(screen, " %s", argv[i]);
+  // The fix is to support the following cases:
+  // 1. Simple atoms, in which case the force points along the velocity
+  // 2. Cases where the atoms have an internal ellipsoid. Currently those
+  //    styles are only body and aspherical particles.
+  // The first argument (mode) is used to differentiate between these.
+
+  // args: fix ID all propel/self mode magnitude
+  // Optional args are
+  const char *mode_str = argv[3];
+  if (strncmp(mode_str, "velocity", 8) == 0) {
+    mode = VELOCITY;
+  } else if (strncmp(mode_str, "quaternion", 10) == 0) {
+    // This mode should only be supported if the atom style has
+    // a quaternion (and if all atoms in the group have it)
+    if (verify_atoms_have_quaternion()) {
+      error->all(FLERR, "All atoms need a quaternion");
     }
-    fprintf(screen, "\n");
+    mode = QUATERNION;
+  } else {
+    char msg[2048];
+    sprintf(msg, "Illegal mode \"%s\" for fix propel/self", mode_str);
+    error->all(FLERR, msg);
   }
 
-  // args: fix ID all propel/self magnitude prop1 prop2 prop3
-  // Optional args are
-  magnitude = force->numeric( FLERR, argv[3] );
+  
+  magnitude = force->numeric( FLERR, argv[4] );
 }
 
 
@@ -92,17 +102,29 @@ double FixPropelSelf::memory_usage()
 }
 
 
-
-void FixPropelSelf::post_force(int /* vflag */ )
+void FixPropelSelf::post_force(int vflag )
 {
-  // Then do the rest:
-  double **f = atom->f;
+	switch(mode) {
+	case QUATERNION:
+		post_force_quaternion(vflag);
+		break;
+	case VELOCITY:
+		post_force_velocity(vflag);
+		break;
+	default:
+		error->all(FLERR, "reached statement that should be unreachable");
+	}
+}
 
+
+
+void FixPropelSelf::post_force_quaternion(int /* vflag */ )
+{
+  double **f = atom->f;
   AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);
 
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
-  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   AtomVecEllipsoid::Bonus *bonus = av->bonus;
   // Add the active force to the atom force:
@@ -141,7 +163,71 @@ void FixPropelSelf::post_force(int /* vflag */ )
       f[i][0] += magnitude * f_rot[0];
       f[i][1] += magnitude * f_rot[1];
       f[i][2] += magnitude * f_rot[2];
-
     }
   }
 }
+
+
+
+void FixPropelSelf::post_force_velocity(int /*vflag*/ )
+{
+  double **f = atom->f;
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  tagint *tag = atom->tag;
+
+  // Add the active force to the atom force:
+  for( int i = 0; i < nlocal; ++i ){
+    if( mask[i] & groupbit ){
+	    const double *vi = v[i];
+	    double f_act[3] = { vi[0], vi[1], vi[2] };
+	    double nv2 = vi[0]*vi[0] + vi[1]*vi[1] + vi[2]*vi[2];
+	    double fnorm = magnitude / sqrt(nv2);
+      
+      if (debug_out && comm->me == 0) {
+        // Magical reference particle:
+        if (tag[i] == 12) {
+          fprintf(screen, "Active force on atom %d: (%f %f %f)\n",
+                  tag[i], f_act[0], f_act[1], f_act[2]);
+          fprintf(screen, "  Total force before: (%f %f %f)\n",
+                  f[i][0], f[i][1], f[i][2]);
+        }
+      }
+
+      f[i][0] += fnorm * f_act[0];
+      f[i][1] += fnorm * f_act[1];
+      f[i][2] += fnorm * f_act[2];
+    }
+  }
+}
+
+
+int FixPropelSelf::verify_atoms_have_quaternion()
+{
+	int ellipsoid_flag = atom->ellipsoid_flag;
+	int body_flag = atom->body_flag;
+	int *mask = atom->mask;
+  if (! (ellipsoid_flag || body_flag) ){
+    error->all(FLERR, "mode quaternion requires body or ellipsoid atoms");
+  }
+  
+  // Make sure all atoms have ellipsoid or body set:
+  for (int i = 0; i < atom->nlocal; ++i) {
+    if (mask[i] & groupbit) {
+	    if (ellipsoid_flag && atom->ellipsoid[i] < 0) {
+		    error->all(FLERR, "Got atom without ellipsoid set");
+		    // Kind-of pointless but silences some compiler warnings:
+		    return 1;
+	    }
+	    if (body_flag && atom->body[i] < 0) {
+		    error->all(FLERR, "Got atom without body set");
+		    // Kind-of pointless but silences some compiler warnings:
+		    return 1;
+	    }
+    }
+  }
+  
+  return 0;
+}
+

src/USER-MISC/fix_propel_self.h

@@ -26,6 +26,12 @@ namespace LAMMPS_NS {
 
 class FixPropelSelf : public Fix {
  public:
+
+  enum operation_modes {
+    VELOCITY = 0,
+    QUATERNION = 1
+  };
+	
   FixPropelSelf(class LAMMPS *, int, char **);
   virtual ~FixPropelSelf();
   virtual int setmask();
@@ -34,9 +40,16 @@ class FixPropelSelf : public Fix {
 
   double memory_usage();
 
-protected:
+private:
   double magnitude;
   int thermostat_orient;
+  int mode;
+
+  int verify_atoms_have_quaternion();
+  void post_force_velocity(int);
+  void post_force_quaternion(int);
+	
+	
 };
 
 }