Unverified Commit f4d9715c authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

make code follow LAMMPS conventions more closely and do some cleanups 

- remove tabs and trailing whitespace
- remove references to atom style body, since code only works with ellipsoid
- adjust function names and tests for requirements to be more obvious and work correctly in parallel
- remove rather specific debug code
- remove non-essential c++11 features
- refactor, correct, and simplify parsing of types keyword arguments
parent d739e017

src/USER-MISC/fix_propel_self.cpp

+94 −127
Original line number Diff line number Diff line
@@ -17,12 +17,12 @@ 
   Thanks to Liesbeth Janssen @ Eindhoven University for useful discussions!

----------------------------------------------------------------------- */





#include <stdio.h>

#include <string.h>



#include "fix_propel_self.h"



#include <cstdio>

#include <cstring>

#include <string>



#include "atom.h"

#include "atom_vec_ellipsoid.h"

#include "citeme.h"
@@ -39,39 +39,43 @@ 
#include "respa.h"

#include "update.h"





using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;



#define PRINT_DEBUG_OUTPUT 0



static constexpr const bool debug_out = false;



/* ---------------------------------------------------------------------- */



FixPropelSelf::FixPropelSelf( LAMMPS *lmp, int narg, char **argv )

	: Fix(lmp, narg, argv), magnitude(1.0), thermostat_orient(0),

  : Fix(lmp, narg, argv), magnitude(0.0),

    mode(VELOCITY), n_types_filter(0), apply_to_type(NULL)

{

  if (narg < 5) error->all(FLERR, "Illegal fix propel/self command");



  // 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.

  // 2. Aspherical particles with an orientation.

  // 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");



    if (!atoms_have_quaternion()) {

      error->all(FLERR, "All fix atoms need to be extended particles");

    }

    mode = QUATERNION;



  } else {

    char msg[2048];

    sprintf(msg, "Illegal mode \"%s\" for fix propel/self", mode_str);
@@ -81,56 +85,47 @@ FixPropelSelf::FixPropelSelf( LAMMPS *lmp, int narg, char **argv ) 
  magnitude = force->numeric( FLERR, argv[4] );



  // Handle rest of args:



  int iarg = 5;

  while (iarg < narg) {

    const char *arg = argv[iarg];

    if (strncmp(arg, "types", 5) == 0) {

      // In this case we need to allocate the type list.

      // First find the largest given type:

      int max_type_in_list = 0;

      ++iarg;

      while ( (iarg + n_types_filter < narg) &&

              (!isalpha(argv[iarg + n_types_filter][0]))) {

        int curr_type = force->numeric(FLERR, argv[iarg + n_types_filter]);

        if (curr_type > max_type_in_list) max_type_in_list = curr_type;

        ++n_types_filter;

      }

      if (n_types_filter == 0) {

        error->all(FLERR, "\"types\" option requires at least one type");

      }



      // sanity check:

      if (max_type_in_list < n_types_filter) {

        error->warning(FLERR, "Found duplicate type in \"types\" option");

      }

    if (strcmp(argv[iarg],"types") == 0) {



      apply_to_type = new int[max_type_in_list+1];

      if (!apply_to_type) error->all(FLERR, "Failed to allocate type storage!");

      apply_to_type = new int[atom->ntypes+1];

      memset(apply_to_type, 0, atom->ntypes * sizeof(int));



      for (int i = 0; i < max_type_in_list+1; ++i) {

        apply_to_type[i] = 0;

      }

      // consume all following numerical arguments as types



      iarg++;

      int flag=0;

      while (iarg < narg) {

        if (isdigit(argv[iarg][0])) {

          int thistype = force->inumeric(FLERR,argv[iarg]);

          if ((thistype < 1) || (thistype > atom->ntypes))

            error->all(FLERR,"Illegal atom type to types keyword");

          apply_to_type[thistype] = 1;

          flag = 1;

          iarg++;

        } else break;

      }

      if (!flag)

        error->all(FLERR,"'types' keyword requires at least one type");

      else

        n_types_filter = 1;



      for (int i = 0; i < n_types_filter; ++i) {

        int curr_type = force->numeric(FLERR, argv[iarg + i]);

        apply_to_type[curr_type] = 1;

      }

      // Done handling types argument.

      iarg += n_types_filter; // -1 because we incremented by 1 previously.

    } else {

      char msg[2048];

      sprintf(msg, "Illegal fix propel/self command: "

              "Unrecognized argument %s", arg);

      error->all(FLERR, msg);

      error->all(FLERR,"Illegal fix propel/self command.");

    }

  }



}



/* ---------------------------------------------------------------------- */



FixPropelSelf::~FixPropelSelf()

{}



{

  delete[] apply_to_type;

}

/* ---------------------------------------------------------------------- */



int FixPropelSelf::setmask()

{
@@ -140,16 +135,18 @@ int FixPropelSelf::setmask() 
  return mask;

}



/* ---------------------------------------------------------------------- */



double FixPropelSelf::memory_usage()

{

  // magnitude + thermostat_orient + mode + n_types_filter + apply_to_type

  double bytes = sizeof(double) + 3*sizeof(int) + sizeof(int*);

  bytes += sizeof(int)*n_types_filter;

  bytes += sizeof(int)*atom->ntypes*n_types_filter;



  return bytes;

}



/* ---------------------------------------------------------------------- */



void FixPropelSelf::post_force(int vflag )

{
@@ -163,11 +160,11 @@ void FixPropelSelf::post_force(int vflag ) 
    else                post_force_velocity<0>(vflag);

    break;

  default:

    error->all(FLERR, "reached statement that should be unreachable");

    ;

  }

}





/* ---------------------------------------------------------------------- */



template <int filter_by_type>

void FixPropelSelf::post_force_quaternion(int /* vflag */ )
@@ -181,7 +178,9 @@ void FixPropelSelf::post_force_quaternion(int /* vflag */ ) 
  int* ellipsoid = atom->ellipsoid;



  AtomVecEllipsoid::Bonus *bonus = av->bonus;



  // Add the active force to the atom force:



  for( int i = 0; i < nlocal; ++i ){

    if( mask[i] & groupbit ){

      if (filter_by_type && !apply_to_type[type[i]]) {
@@ -192,31 +191,11 @@ void FixPropelSelf::post_force_quaternion(int /* vflag */ ) 
      double f_rot[3];



      double *quat  = bonus[ellipsoid[i]].quat;

      tagint *tag = atom->tag;



      double Q[3][3];

      MathExtra::quat_to_mat( quat, Q );

      MathExtra::matvec( Q, f_act, f_rot );



      if (debug_out && comm->me == 0) {

        // Magical reference particle:

        if (tag[i] == 12) {

          fprintf(screen, "rotation quaternion: (%f %f %f %f)\n",

                  quat[0], quat[1], quat[2], quat[3]);

          fprintf(screen, "rotation matrix: / %f %f %f \\\n",

                  Q[0][0] ,Q[0][1], Q[0][2]);

          fprintf(screen, "                 | %f %f %f |\n",

                  Q[1][0] ,Q[1][1], Q[1][2]);

          fprintf(screen, "                 \\ %f %f %f /\n",

                  Q[2][0] ,Q[2][1], Q[2][2]);



          fprintf(screen, "Active force on atom %d: (%f %f %f)\n",

                  tag[i], f_rot[0], f_rot[1], f_rot[2]);

          fprintf(screen, "  Total force before: (%f %f %f)\n",

                  f[i][0], f[i][1], f[i][2]);

        }

      }



      f[i][0] += magnitude * f_rot[0];

      f[i][1] += magnitude * f_rot[1];

      f[i][2] += magnitude * f_rot[2];
@@ -224,6 +203,7 @@ void FixPropelSelf::post_force_quaternion(int /* vflag */ ) 
  }

}



/* ---------------------------------------------------------------------- */



template <int filter_by_type>

void FixPropelSelf::post_force_velocity(int /*vflag*/)
@@ -232,10 +212,10 @@ void FixPropelSelf::post_force_velocity(int /*vflag*/ ) 
  double **v = atom->v;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  tagint *tag = atom->tag;

  int *type = atom->type;



  // Add the active force to the atom force:



  for(int i = 0; i < nlocal; ++i) {

    if( mask[i] & groupbit ){

      if (filter_by_type && !apply_to_type[type[i]]) {
@@ -246,21 +226,14 @@ void FixPropelSelf::post_force_velocity(int /*vflag*/ ) 
      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;

      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:

        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]);

        }

        // Without this check you can run into numerical

        // issues because fnorm will blow up.



        fnorm = magnitude / sqrt(nv2);

      }



      f[i][0] += fnorm * f_act[0];
@@ -270,32 +243,26 @@ void FixPropelSelf::post_force_velocity(int /*vflag*/ ) 
  }

}



/* ---------------------------------------------------------------------- */



int FixPropelSelf::verify_atoms_have_quaternion()

int FixPropelSelf::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");

  if (!atom->ellipsoid_flag) {

    error->all(FLERR, "Mode 'quaternion' requires atom style ellipsoid");

    return 0;

  }



  // 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;

      }

    }

  }

  int *mask = atom->mask;

  int flag=0,flagall=0;



  return 0;

}

  // Make sure all atoms have ellipsoid data:



  for (int i = 0; i < atom->nlocal; ++i)

    if (mask[i] & groupbit)

      if (atom->ellipsoid[i] < 0) ++flag;



  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);

  if (flagall > 0) return 0;

  

  return 1;

}

src/USER-MISC/fix_propel_self.h

+8 −13
Original line number Diff line number Diff line
@@ -27,39 +27,34 @@ 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();

  virtual void post_force(int);

  // virtual void post_force_respa(int, int, int);



  double memory_usage();



 protected:

  enum operation_modes {

    VELOCITY = 0,

    QUATERNION = 1

  };



private:

  double magnitude;

  int thermostat_orient;

  int mode;





  // If 0, apply fix to everything in group. If > 0, apply only to those

  // types i for which i <= n_types_filter _and_ apply_to_type[i] == 1:

  int n_types_filter;

  int *apply_to_type; //< Specifies, per type, if the fix applies to it or not.





  int verify_atoms_have_quaternion();

  int atoms_have_quaternion();



  template <int filter_by_type> void post_force_velocity(int);

  template <int filter_by_type> void post_force_quaternion(int);

	

	

};



}



#endif

CRA Git | Maintained and supported by SUSTech CRA and CCSE