Commit 6e945f72 authored by Charlles Abreu's avatar Charlles Abreu
Browse files

Keyword 'bodyforces' added to fix_modify command

- Applies to classes FixRigid, FixRigidSmall, and FixPOEMS
parent 251bc882
Loading
Loading
Loading
Loading
+64 −19
Original line number Diff line number Diff line
@@ -280,6 +280,10 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
      fprintf(logfile,"%d clusters, %d bodies, %d joints, %d atoms\n",
              ncluster,nbody,njoint,nsum);
  }

  // compute per body forces and torques inside final_integrate() by default

  earlyflag = 0;
}

/* ----------------------------------------------------------------------
@@ -330,6 +334,25 @@ FixPOEMS::~FixPOEMS()

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

int FixPOEMS::modify_param(int narg, char **arg)
{
  if (strcmp(arg[0],"bodyforces") == 0) {
    if (narg < 2)
      error->all(FLERR,"Illegal fix_modify command");
    if (strcmp(arg[1],"early") == 0)
      earlyflag = 1;
    else if (strcmp(arg[1],"late") == 0)
      earlyflag = 0;
    else
      error->all(FLERR,"Illegal fix_modify command");
    return 2;
  }
  else
    return 0;
}

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

int FixPOEMS::setmask()
{
  int mask = 0;
@@ -352,9 +375,17 @@ void FixPOEMS::init()
  // warn if more than one POEMS fix

  int count = 0;
  for (int i = 0; i < modify->nfix; i++)
    if (strcmp(modify->fix[i]->style,"poems") == 0) count++;
  if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix poems");
  int myindex, myposition;
  for (i = 0; i < modify->nfix; i++)
    if (strcmp(modify->fix[i]->style,"poems") == 0) {
      count++;
      if (strcmp(modify->fix[i]->id,id) != 0) {
        myindex = i;
        myposition = count;
      }
    }
  if (count > 1 && myposition == 1 && comm->me == 0)
    error->warning(FLERR,"More than one fix poems");

  // error if npt,nph fix comes before rigid fix

@@ -368,20 +399,25 @@ void FixPOEMS::init()
        error->all(FLERR,"POEMS fix must come before NPT/NPH fix");
  }

  // warn if any non-rigid fix with post_force succeeds any fix rigid:
  int first_rigid = 0;
  while (!modify->fix[first_rigid]->rigid_flag)
    first_rigid++;
  // warn if fix poems preceeds non-rigid fixes with post-force tasks
  // when computing body forces and torques in post_force() instead of final_integrate()

  if (earlyflag) {
    int has_post_force[modify->nfix - myindex];
    count = 0;
  for (i = first_rigid + 1; i < modify->nfix; i++) {
    Fix *ifix = modify->fix[i];
    if ( (modify->fmask[i] & POST_FORCE) && (!ifix->rigid_flag) ) {
      count++;
      if (comm->me == 0) {
      if (count == 1)
        error->warning(FLERR,"The following fixes must preceed any rigid-body fixes");
        if (screen) fprintf(screen,"> fix %s %s\n",ifix->id,ifix->style);
        if (logfile) fprintf(logfile,"> fix %s %s\n",ifix->id,ifix->style);
    for (i = myindex + 1; i < modify->nfix; i++)
      if ( (modify->fmask[i] & POST_FORCE) && (!modify->fix[i]->rigid_flag) )
        has_post_force[count++] = i;
    if (count) {
      FILE *p[2] = {screen, logfile};
      for (int j = 0; j < 2; j++)
        if (p[j]) {
          fprintf(p[j],"WARNING: fix %s %s",id,style);
          fprintf(p[j]," will add up forces before they are handled by:\n");
          for (int k = 0; k < count; k++) {
            Fix *fix = modify->fix[has_post_force[k]];
            fprintf(p[j],"         => fix %s %s\n",fix->id,fix->style);
          }
        }
    }
  }
@@ -747,7 +783,7 @@ void FixPOEMS::initial_integrate(int vflag)
   only count joint atoms in 1st body
------------------------------------------------------------------------- */

void FixPOEMS::post_force(int vflag)
void FixPOEMS::compute_forces_and_torques()
{
  int i,ibody;
  int xbox,ybox,zbox;
@@ -798,6 +834,13 @@ void FixPOEMS::post_force(int vflag)
  }
}

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

void FixPOEMS::post_force(int vflag)
{
  if (earlyflag) compute_forces_and_torques();
}

/* ----------------------------------------------------------------------
   update vcm,omega by last 1/2 step
   set v of body atoms accordingly
@@ -805,6 +848,8 @@ void FixPOEMS::post_force(int vflag)

void FixPOEMS::final_integrate()
{
  if (!earlyflag) compute_forces_and_torques();

  // perform POEMS integration

  poems->LobattoTwo(vcm,omega,torque,fcm);
+4 −0
Original line number Diff line number Diff line
@@ -28,6 +28,7 @@ class FixPOEMS : public Fix {
 public:
  FixPOEMS(class LAMMPS *, int narg, char **arg);
  ~FixPOEMS();
  int modify_param(int, char **);
  int setmask();
  void init();
  void setup(int);
@@ -56,6 +57,8 @@ class FixPOEMS : public Fix {
  int nlevels_respa;
  double total_ke;

  int earlyflag;    // 1 if forces and torques are computed at post_force()

  // atom assignment to rigid bodies
  // double count joint atoms as being in multiple bodies

@@ -96,6 +99,7 @@ class FixPOEMS : public Fix {

  // internal class functions

  void compute_forces_and_torques();
  void readfile(char *);
  int readline(FILE *, char **, int *);
  void jointbuild();
+52 −34
Original line number Diff line number Diff line
@@ -639,6 +639,25 @@ FixRigid::~FixRigid()

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

int FixRigid::modify_param(int narg, char **arg)
{
  if (strcmp(arg[0],"bodyforces") == 0) {
    if (narg < 2)
      error->all(FLERR,"Illegal fix_modify command");
    if (strcmp(arg[1],"early") == 0)
      earlyflag = 1;
    else if (strcmp(arg[1],"late") == 0)
      earlyflag = 0;
    else
      error->all(FLERR,"Illegal fix_modify command");
    return 2;
  }
  else
    return 0;
}

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

int FixRigid::setmask()
{
  int mask = 0;
@@ -668,9 +687,17 @@ void FixRigid::init()
  // warn if more than one rigid fix

  int count = 0;
  int myindex, myposition;
  for (i = 0; i < modify->nfix; i++)
    if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
  if (count > 1 && me == 0) error->warning(FLERR,"More than one fix rigid");
    if (strcmp(modify->fix[i]->style,"rigid") == 0) {
      count++;
      if (strcmp(modify->fix[i]->id,id) != 0) {
        myindex = i;
        myposition = count;
      }
    }
  if (count > 1 && myposition == 1 && comm->me == 0)
    error->warning(FLERR,"More than one fix rigid");

  // error if npt,nph fix comes before rigid fix

@@ -684,23 +711,24 @@ void FixRigid::init()
        error->all(FLERR,"Rigid fix must come before NPT/NPH fix");
  }

  // warn if any non-rigid fix with post_force succeeds any fix rigid instance
  // warn if fix rigid preceeds non-rigid fixes with post-force tasks
  // when computing body forces and torques in post_force() instead of final_integrate()

  if (earlyflag) {
    int first_rigid = 0;
    while (!modify->fix[first_rigid]->rigid_flag)
      first_rigid++;
    int has_post_force[modify->nfix - myindex];
    count = 0;
    for (i = first_rigid + 1; i < modify->nfix; i++) {
      Fix *ifix = modify->fix[i];
      if ( (modify->fmask[i] & POST_FORCE) && (!ifix->rigid_flag) ) {
        count++;
        if (comm->me == 0) {
          if (count == 1)
            error->warning(FLERR,"The following fixes must preceed any rigid-body fixes");
          if (screen) fprintf(screen,"> fix %s %s\n",ifix->id,ifix->style);
          if (logfile) fprintf(logfile,"> fix %s %s\n",ifix->id,ifix->style);
    for (i = myindex + 1; i < modify->nfix; i++)
      if ( (modify->fmask[i] & POST_FORCE) && (!modify->fix[i]->rigid_flag) )
        has_post_force[count++] = i;
    if (count) {
      FILE *p[2] = {screen, logfile};
      for (int j = 0; j < 2; j++)
        if (p[j]) {
          fprintf(p[j],"WARNING: fix %s %s",id,style);
          fprintf(p[j]," will add up forces before they are handled by:\n");
          for (int k = 0; k < count; k++) {
            Fix *fix = modify->fix[has_post_force[k]];
            fprintf(p[j],"         => fix %s %s\n",fix->id,fix->style);
          }
        }
    }
@@ -1054,8 +1082,10 @@ void FixRigid::compute_forces_and_torques()

void FixRigid::post_force(int vflag)
{
  if (earlyflag) {
    if (langflag) apply_langevin_thermostat();
  if (earlyflag) compute_forces_and_torques();
    compute_forces_and_torques();
  }
}

/* ---------------------------------------------------------------------- */
@@ -1065,7 +1095,10 @@ void FixRigid::final_integrate()
  int ibody;
  double dtfm;

  if (!earlyflag) compute_forces_and_torques();
  if (!earlyflag) {
    if (langflag) apply_langevin_thermostat();
    compute_forces_and_torques();
  }

  // update vcm and angmom
  // fflag,tflag = 0 for some dimensions in 2d
@@ -2701,18 +2734,3 @@ double FixRigid::compute_array(int i, int j)
  if (j == 13) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX;
  return (imagebody[i] >> IMG2BITS) - IMGMAX;
}

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

int FixRigid::modify_param(int narg, char **arg)
{
  if (strcmp(arg[0],"bodyforces") == 0) {
    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
    if (strcmp(arg[1],"early") == 0) earlyflag=1;
    else if (strcmp(arg[1],"late") == 0) earlyflag=0;
    else error->all(FLERR,"Illegal fix_modify command");

    return 2;
  }
  return 0;
}
+2 −2
Original line number Diff line number Diff line
@@ -28,6 +28,7 @@ class FixRigid : public Fix {
 public:
  FixRigid(class LAMMPS *, int, char **);
  virtual ~FixRigid();
  int modify_param(int, char **);
  virtual int setmask();
  virtual void init();
  virtual void setup(int);
@@ -38,7 +39,6 @@ class FixRigid : public Fix {
  void final_integrate_respa(int, int);
  void write_restart_file(char *);
  virtual double compute_scalar();
  virtual int modify_param(int, char **);

  double memory_usage();
  void grow_arrays(int);
@@ -72,6 +72,7 @@ class FixRigid : public Fix {
  char *infile;             // file to read rigid body attributes from
  int rstyle;               // SINGLE,MOLECULE,GROUP
  int setupflag;            // 1 if body properties are setup, else 0
  int earlyflag;     // 1 if forces and torques are computed at post_force()

  int dimension;            // # of dimensions
  int nbody;                // # of rigid bodies
@@ -107,7 +108,6 @@ class FixRigid : public Fix {
  int orientflag;           // 1 if particles store spatial orientation
  int dorientflag;          // 1 if particles store dipole orientation
  int reinitflag;           // 1 if re-initialize rigid bodies between runs
  int earlyflag;            // 1 if compute body forces and torques early, i.e. in post_force()

  imageint *xcmimage;       // internal image flags for atoms in rigid bodies
                            // set relative to in-box xcm of each body
+2 −3
Original line number Diff line number Diff line
@@ -1438,9 +1438,8 @@ int FixRigidNHSmall::modify_param(int narg, char **arg)
    if (pressure->pressflag == 0)
      error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
    return 2;
  }

  return 0;
  } else
    return FixRigidSmall::modify_param(narg,arg);
}

/* ---------------------------------------------------------------------- */
Loading