Commit f56a345a authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

implement "earlyflag" to select when to compute body torques/forces

parent 76cc545d

src/RIGID/fix_rigid.cpp

+40 −15
Original line number Diff line number Diff line
@@ -575,6 +575,10 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : 


  setupflag = 0;



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



  earlyflag = 0;



  // print statistics



  int nsum = 0;
@@ -680,7 +684,10 @@ void FixRigid::init() 
        error->all(FLERR,"Rigid fix must come before NPT/NPH fix");

  }



  // error if any non-rigid fix with post_force succeeds any fix rigid:

  // warn if any non-rigid fix with post_force succeeds any fix rigid instance

  // 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++;
@@ -690,13 +697,14 @@ void FixRigid::init() 
      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);

        }

      }

    }

  if (count > 0)

    error->all(FLERR,"the fixes listed above must preceed all rigid-body fixes");

  }



  // timestep info
@@ -1047,7 +1055,7 @@ void FixRigid::compute_forces_and_torques() 
void FixRigid::post_force(int vflag)

{

  if (langflag) apply_langevin_thermostat();

  compute_forces_and_torques();

  if (earlyflag) compute_forces_and_torques();

}



/* ---------------------------------------------------------------------- */
@@ -1057,6 +1065,8 @@ void FixRigid::final_integrate() 
  int ibody;

  double dtfm;



  if (!earlyflag) compute_forces_and_torques();



  // update vcm and angmom

  // fflag,tflag = 0 for some dimensions in 2d
@@ -2691,3 +2701,18 @@ 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;

}

src/RIGID/fix_rigid.h

+2 −1
Original line number Diff line number Diff line
@@ -38,7 +38,7 @@ 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 **) {return 0;}

  virtual int modify_param(int, char **);



  double memory_usage();

  void grow_arrays(int);
@@ -107,6 +107,7 @@ 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

src/RIGID/fix_rigid_small.cpp

+44 −19
Original line number Diff line number Diff line
@@ -403,6 +403,10 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : 
  avec_line = (AtomVecLine *) atom->style_match("line");

  avec_tri = (AtomVecTri *) atom->style_match("tri");



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



  earlyflag = 0;



  // print statistics



  int one = 0;
@@ -516,7 +520,10 @@ void FixRigidSmall::init() 
        error->all(FLERR,"Rigid fix must come before NPT/NPH fix");

  }



  // error if any non-rigid fix with post_force succeeds any fix rigid:

  // warn if any non-rigid fix with post_force succeeds any fix rigid instance

  // 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++;
@@ -526,13 +533,14 @@ void FixRigidSmall::init() 
      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);

        }

       }

     }

  if (count > 0)

    error->all(FLERR,"the fixes listed above must preceed all rigid-body fixes");

   }



  // timestep info
@@ -901,7 +909,7 @@ void FixRigidSmall::compute_forces_and_torques() 
void FixRigidSmall::post_force(int vflag)

{

  if (langflag) apply_langevin_thermostat();

  compute_forces_and_torques();

  if (earlyflag) compute_forces_and_torques();

}



/* ---------------------------------------------------------------------- */
@@ -910,7 +918,9 @@ void FixRigidSmall::final_integrate() 
{

  double dtfm;



  //check(3);

  //check(3)



  if (!earlyflag) compute_forces_and_torques();



  // update vcm and angmom, recompute omega
@@ -3500,6 +3510,21 @@ double FixRigidSmall::compute_scalar() 
  return tall;

}



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



int FixRigidSmall::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;

}



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

   memory usage of local atom-based arrays

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

src/RIGID/fix_rigid_small.h

+3 −1
Original line number Diff line number Diff line
@@ -37,7 +37,7 @@ class FixRigidSmall : public Fix { 
  virtual void init();

  virtual void setup(int);

  virtual void initial_integrate(int);

  void post_force(int);

  virtual void post_force(int);

  virtual void final_integrate();

  void initial_integrate_respa(int, int, int);

  void final_integrate_respa(int, int);
@@ -69,6 +69,7 @@ class FixRigidSmall : public Fix { 
  double extract_ke();

  double extract_erotational();

  double compute_scalar();

  virtual int modify_param(int, char **);

  double memory_usage();



 protected:
@@ -130,6 +131,7 @@ class FixRigidSmall : 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()



  int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE;   // bitmasks for eflags

  int OMEGA,ANGMOM,TORQUE;

src/fix.h

+0 −1
Original line number Diff line number Diff line
@@ -102,7 +102,6 @@ class Fix : protected Pointers { 


  Fix(class LAMMPS *, int, char **);

  virtual ~Fix();

  void modify_params(int, char **);



  virtual int setmask() = 0;

CRA Git | Maintained and supported by SUSTech CRA and CCSE