Some code clean-up, added safety check in post_force_velocity.

examples/USER/misc/fix_propel_self/2d_velocity/in.2d_langevin

@@ -0,0 +1,37 @@
+dimension 2
+boundary p p p
+
+variable L equal 20
+region total block -$L $L -$L $L -0.5 0.5
+lattice hex 0.3
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+mass * 1.0
+
+fix step all nve
+fix temp all langevin 1.0 1.0 1.0 13
+fix twod all enforce2d
+
+neighbor 0.6 bin
+
+dump traj all custom 250 2d_active.dump.bin id type x y z
+
+thermo_style custom time step pe ke etotal temp
+thermo 1000
+run 25000
+
+fix active all propel/self velocity 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 25000

examples/USER/misc/fix_propel_self/2d_velocity/in.2d_viscous

@@ -0,0 +1,37 @@
+dimension 2
+boundary p p p
+
+variable L equal 20
+region total block -$L $L -$L $L -0.5 0.5
+lattice hex 0.3
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+mass * 1.0
+
+fix step all nve
+fix twod all enforce2d
+
+neighbor 0.6 bin
+
+dump traj all custom 250 2d_active.dump.bin id type x y z
+
+thermo_style custom time step pe ke etotal temp
+thermo 1000
+run 25000
+
+fix active all propel/self velocity 1.0
+fix fric all viscous 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 25000

examples/USER/misc/fix_propel_self/3d_quaternion/in.3d_quaternion

@@ -0,0 +1,40 @@
+dimension 3
+boundary p p p
+
+atom_style ellipsoid
+variable L equal 20
+region total block -$L $L -$L $L -$L $L
+lattice sc 0.1
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+# mass * 1.0
+set type * density 1.0
+set type * shape 1.0 1.0 1.0
+set type * quat 0 0 1 0
+
+fix step all nve
+fix temp all langevin 1.0 1.0 1.0 13
+
+neighbor 0.6 bin
+
+dump traj all custom 100 3d_active.dump.bin id type x y z fx fy fz
+
+thermo_style custom time step pe ke etotal temp
+thermo 100
+run 5000
+
+fix active all propel/self quaternion 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 5000

src/USER-MISC/fix_propel_self.cpp

@@ -103,16 +103,16 @@ double FixPropelSelf::memory_usage()
 
 void FixPropelSelf::post_force(int vflag )
 {
-	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");
-	}
+  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");
+  }
 }
 
 
@@ -179,10 +179,16 @@ void FixPropelSelf::post_force_velocity(int /*vflag*/ )
   // 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);
+      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 = 0.0;
+      constexpr const double TOL = 1e-14;
+      if (nv2 > TOL) {
+        // Without this check you can run into numerical issues
+        // because fnorm will blow up.
+        fnorm = magnitude / sqrt(nv2);
+      }
       
       if (debug_out && comm->me == 0) {
         // Magical reference particle:
@@ -214,16 +220,16 @@ int FixPropelSelf::verify_atoms_have_quaternion()
   // 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 return but silences compiler warnings:
-		    return 1;
-	    }
-	    if (body_flag && atom->body[i] < 0) {
-		    error->all(FLERR, "Got atom without body set");
-		    // Kind-of pointless return silences compiler warnings:
-		    return 1;
-	    }
+      if (ellipsoid_flag && atom->ellipsoid[i] < 0) {
+        error->all(FLERR, "Got atom without ellipsoid set");
+        // Kind-of pointless return but silences compiler warnings:
+        return 1;
+      }
+      if (body_flag && atom->body[i] < 0) {
+        error->all(FLERR, "Got atom without body set");
+        // Kind-of pointless return silences compiler warnings:
+        return 1;
+      }
     }
   }