add support to fix adapt for angle coeffs

doc/src/fix_adapt.rst

@@ -14,7 +14,7 @@ Syntax
 * adapt = style name of this fix command
 * N = adapt simulation settings every this many timesteps
 * one or more attribute/arg pairs may be appended
-* attribute = *pair* or *bond* or *kspace* or *atom*
+* attribute = *pair* or *bond* or *angle* or *kspace* or *atom*
 
   .. parsed-literal::
 
@@ -28,11 +28,16 @@ Syntax
          bparam = parameter to adapt over time
          I = type bond to set parameter for
          v_name = variable with name that calculates value of bparam
+       *angle* args = astyle aparam I v_name
+         astyle = angle style name, e.g. harmonic
+         aparam = parameter to adapt over time
+         I = type angle to set parameter for
+         v_name = variable with name that calculates value of aparam
        *kspace* arg = v_name
          v_name = variable with name that calculates scale factor on K-space terms
-       *atom* args = aparam v_name
+       *atom* args = atomparam v_name
-         aparam = parameter to adapt over time
+         atomparam = parameter to adapt over time
-         v_name = variable with name that calculates value of aparam
+         v_name = variable with name that calculates value of atomparam
 
 * zero or more keyword/value pairs may be appended
 * keyword = *scale* or *reset* or *mass*
@@ -283,30 +288,60 @@ operates. The only difference is that now a bond coefficient for a
 given bond type is adapted.
 
 A wild-card asterisk can be used in place of or in conjunction with
-the bond type argument to set the coefficients for multiple bond types.
+the bond type argument to set the coefficients for multiple bond
-This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N = the number of
+types.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N =
-atom types, then an asterisk with no numeric values means all types
+the number of bond types, then an asterisk with no numeric values
-from 1 to N.  A leading asterisk means all types from 1 to n (inclusive).
+means all types from 1 to N.  A leading asterisk means all types from
-A trailing asterisk means all types from n to N (inclusive).  A middle
+1 to n (inclusive).  A trailing asterisk means all types from n to N
-asterisk means all types from m to n (inclusive).
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).
 
 Currently *bond* does not support bond_style hybrid nor bond_style
 hybrid/overlay as bond styles. The only bonds that currently are
 working with fix_adapt are
 
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`class2 <bond_class2>`        | r0    | type bonds |
+| :doc:`class2 <bond_class2>`        | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`fene <bond_fene>`            | k, r0 | type bonds |
+| :doc:`fene <bond_fene>`            | k, r0      | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`gromos <bond_gromos>`        | k, r0 | type bonds |
+| :doc:`fene/nm <bond_fene_nm>`      | k, r0      | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`harmonic <bond_harmonic>`    | k,r0  | type bonds |
+| :doc:`gromos <bond_gromos>`        | k, r0      | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`morse <bond_morse>`          | r0    | type bonds |
+| :doc:`harmonic <bond_harmonic>`    | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`nonlinear <bond_nonlinear>`  | r0    | type bonds |
+| :doc:`morse <bond_morse>`          | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
+| :doc:`nonlinear <bond_nonlinear>`  | epsilon,r0 | type bonds |
++------------------------------------+-------+-----------------+
+
+----------
+
+The *angle* keyword uses the specified variable to change the value of
+an angle coefficient over time, very similar to how the *pair* keyword
+operates. The only difference is that now an angle coefficient for a
+given angle type is adapted.
+
+A wild-card asterisk can be used in place of or in conjunction with
+the angle type argument to set the coefficients for multiple angle
+types.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N =
+the number of angle types, then an asterisk with no numeric values
+means all types from 1 to N.  A leading asterisk means all types from
+1 to n (inclusive).  A trailing asterisk means all types from n to N
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).
+
+Currently *angle* does not support angle_style hybrid nor angle_style
+hybrid/overlay as angle styles. The only angles that currently are
+working with fix_adapt are
+
++------------------------------------+-------+-----------------+
+| :doc:`harmonic <angle_harmonic>`    | k,theta0 | type angles |
++------------------------------------+-------+-----------------+
+
+Note that internally, theta0 is stored in radians, so the variable
+this fix uses to reset theta0 needs to generate values in radians.
 
 ----------
 

src/EXTRA-MOLECULE/bond_fene_nm.cpp

@@ -273,7 +273,7 @@ double BondFENENM::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
 void *BondFENENM::extract(const char *str, int &dim)
 {
   dim = 1;
-  if (strcmp(str, "kappa") == 0) return (void *) k;
+  if (strcmp(str, "k") == 0) return (void *) k;
   if (strcmp(str, "r0") == 0) return (void *) r0;
   return nullptr;
 }

src/EXTRA-MOLECULE/bond_gaussian.cpp

@@ -34,9 +34,7 @@ using namespace MathConst;
 BondGaussian::BondGaussian(LAMMPS *lmp) :
     Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr),
     r0(nullptr)
-{
+{}
-  reinitflag = 1;
-}
 
 /* ---------------------------------------------------------------------- */
 

src/MOLECULE/angle_harmonic.cpp

@@ -264,3 +264,15 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3)
   double tk = k[type] * dtheta;
   return tk * dtheta;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleHarmonic::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/MOLECULE/angle_harmonic.h

@@ -35,6 +35,7 @@ class AngleHarmonic : public Angle {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *theta0;

src/MOLECULE/bond_fene.cpp

@@ -265,7 +265,7 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/, double &ffor
 void *BondFENE::extract(const char *str, int &dim)
 {
   dim = 1;
-  if (strcmp(str, "kappa") == 0) return (void *) k;
+  if (strcmp(str, "k") == 0) return (void *) k;
   if (strcmp(str, "r0") == 0) return (void *) r0;
   return nullptr;
 }

src/MOLECULE/bond_gromos.cpp

@@ -30,10 +30,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp)
+BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) {}
-{
-  reinitflag = 1;
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -200,7 +197,7 @@ double BondGromos::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
 void *BondGromos::extract(const char *str, int &dim)
 {
   dim = 1;
-  if (strcmp(str, "kappa") == 0) return (void *) k;
+  if (strcmp(str, "k") == 0) return (void *) k;
   if (strcmp(str, "r0") == 0) return (void *) r0;
   return nullptr;
 }

src/MOLECULE/bond_harmonic.cpp

@@ -27,10 +27,7 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp)
+BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {}
-{
-  reinitflag = 1;
-}
 
 /* ---------------------------------------------------------------------- */
 
@@ -201,12 +198,13 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, double &
 }
 
 /* ----------------------------------------------------------------------
-    Return ptr to internal members upon request.
+   return ptr to internal members upon request
 ------------------------------------------------------------------------ */
+
 void *BondHarmonic::extract(const char *str, int &dim)
 {
   dim = 1;
-  if (strcmp(str, "kappa") == 0) return (void *) k;
+  if (strcmp(str, "k") == 0) return (void *) k;
   if (strcmp(str, "r0") == 0) return (void *) r0;
   return nullptr;
 }

src/angle.cpp

@@ -353,3 +353,15 @@ double Angle::memory_usage()
   bytes += (double) comm->nthreads * maxcvatom * 9 * sizeof(double);
   return bytes;
 }
+
+/* -----------------------------------------------------------------------
+   reset all type-based angle params via init()
+-------------------------------------------------------------------------- */
+
+void Angle::reinit()
+{
+  if (!reinitflag) 
+    error->all(FLERR, "Fix adapt interface to this angle style not supported");
+
+  init();
+}

src/angle.h

@@ -36,6 +36,9 @@ class Angle : protected Pointers {
                              // CENTROID_AVAIL = different and implemented
                              // CENTROID_NOTAVAIL = different, not yet implemented
 
+  int reinitflag;            // 0 if not compatible with fix adapt
+                             // extract() method may still need to be added
+
   // KOKKOS host/device flag and data masks
 
   ExecutionSpace execution_space;
@@ -57,6 +60,8 @@ class Angle : protected Pointers {
   virtual void write_data(FILE *) {}
   virtual double single(int, int, int, int) = 0;
   virtual double memory_usage();
+  virtual void *extract(const char *, int &) { return nullptr; }
+  void reinit();
 
  protected:
   int suffix_flag;    // suffix compatibility flag

src/bond.cpp

@@ -43,6 +43,7 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp)
   energy = 0.0;
   virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
   writedata = 1;
+  reinitflag = 1;
 
   comm_forward = comm_reverse = comm_reverse_off = 0;
 
@@ -336,11 +337,13 @@ double Bond::memory_usage()
 }
 
 /* -----------------------------------------------------------------------
-   Reset all type-based bond params via init.
+   reset all type-based bond params via init()
 -------------------------------------------------------------------------- */
+
 void Bond::reinit()
 {
-  if (!reinitflag) error->all(FLERR, "Fix adapt interface to this bond style not supported");
+  if (!reinitflag) 
+    error->all(FLERR, "Fix adapt interface to this bond style not supported");
 
   init();
 }

src/bond.h

@@ -37,7 +37,8 @@ class Bond : protected Pointers {
   int comm_reverse;          // size of reverse communication (0 if none)
   int comm_reverse_off;      // size of reverse comm even if newton off
 
-  int reinitflag;    // 1 if compatible with fix adapt and alike
+  int reinitflag;            // 0 if not compatible with fix adapt
+                             // extract() method may still need to be added
 
   // KOKKOS host/device flag and data masks
 
@@ -61,7 +62,8 @@ class Bond : protected Pointers {
   virtual double single(int, double, int, int, double &) = 0;
   virtual double memory_usage();
   virtual void *extract(const char *, int &) { return nullptr; }
-  virtual void reinit();
+  void reinit();
+
   virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;}
   virtual void unpack_forward_comm(int, int, double *) {}
   virtual int pack_reverse_comm(int, int, double *) {return 0;}

src/fix_adapt.cpp

@@ -14,6 +14,7 @@
 
 #include "fix_adapt.h"
 
+#include "angle.h"
 #include "atom.h"
 #include "bond.h"
 #include "domain.h"
@@ -38,7 +39,7 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{PAIR,KSPACE,ATOM,BOND};
+enum{PAIR,KSPACE,ATOM,BOND,ANGLE};
 enum{DIAMETER,CHARGE};
 
 /* ---------------------------------------------------------------------- */
@@ -75,6 +76,10 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
       if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command");
       nadapt++;
       iarg += 5;
+    } else if (strcmp(arg[iarg],"angle") == 0) {
+      if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command");
+      nadapt++;
+      iarg += 5;
     } else break;
   }
 
@@ -119,6 +124,20 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
       nadapt++;
       iarg += 5;
 
+    } else if (strcmp(arg[iarg],"angle") == 0) {
+      if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command");
+      adapt[nadapt].which = ANGLE;
+      adapt[nadapt].angle = nullptr;
+      adapt[nadapt].astyle = utils::strdup(arg[iarg+1]);
+      adapt[nadapt].aparam = utils::strdup(arg[iarg+2]);
+      utils::bounds(FLERR,arg[iarg+3],1,atom->nangletypes,
+                    adapt[nadapt].ilo,adapt[nadapt].ihi,error);
+      if (utils::strmatch(arg[iarg+4],"^v_")) {
+        adapt[nadapt].var = utils::strdup(arg[iarg+4]+2);
+      } else error->all(FLERR,"Illegal fix adapt command");
+      nadapt++;
+      iarg += 5;
+
     } else if (strcmp(arg[iarg],"kspace") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command");
       adapt[nadapt].which = KSPACE;
@@ -133,12 +152,12 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
       adapt[nadapt].which = ATOM;
       if (strcmp(arg[iarg+1],"diameter") == 0 ||
           strcmp(arg[iarg+1],"diameter/disc") == 0) {
-        adapt[nadapt].aparam = DIAMETER;
+        adapt[nadapt].atomparam = DIAMETER;
         diamflag = 1;
         discflag = 0;
         if (strcmp(arg[iarg+1],"diameter/disc") == 0) discflag = 1;
       } else if (strcmp(arg[iarg+1],"charge") == 0) {
-        adapt[nadapt].aparam = CHARGE;
+        adapt[nadapt].atomparam = CHARGE;
         chgflag = 1;
       } else error->all(FLERR,"Illegal fix adapt command");
       if (utils::strmatch(arg[iarg+2],"^v_")) {
@@ -191,6 +210,13 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
   for (int m = 0; m < nadapt; ++m)
     if (adapt[m].which == BOND)
       memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
+
+  // allocate angle style arrays:
+
+  n = atom->nbondtypes;
+  for (int m = 0; m < nadapt; ++m)
+    if (adapt[m].which == ANGLE)
+      memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -207,6 +233,10 @@ FixAdapt::~FixAdapt()
       delete [] adapt[m].bstyle;
       delete [] adapt[m].bparam;
       memory->destroy(adapt[m].vector_orig);
+    } else if (adapt[m].which == ANGLE) {
+      delete [] adapt[m].astyle;
+      delete [] adapt[m].aparam;
+      memory->destroy(adapt[m].vector_orig);
     }
   }
   delete [] adapt;
@@ -357,6 +387,7 @@ void FixAdapt::init()
       }
 
       delete [] pstyle;
+
     } else if (ad->which == BOND) {
       ad->bond = nullptr;
       anybond = 1;
@@ -383,13 +414,39 @@ void FixAdapt::init()
 
       delete [] bstyle;
 
+    } else if (ad->which == ANGLE) {
+      ad->angle = nullptr;
+      anyangle = 1;
+
+      char *astyle = utils::strdup(ad->astyle);
+      if (lmp->suffix_enable)
+        ad->angle = force->angle_match(fmt::format("{}/{}",astyle,lmp->suffix));
+
+      if (ad->angle == nullptr) ad->angle = force->angle_match(astyle);
+      if (ad->angle == nullptr )
+        error->all(FLERR,"Fix adapt angle style does not exist");
+
+      void *ptr = ad->angle->extract(ad->aparam,ad->bdim);
+
+      if (ptr == nullptr)
+        error->all(FLERR,"Fix adapt angle style param not supported");
+
+      // for angle styles, use a vector
+
+      if (ad->adim == 1) ad->vector = (double *) ptr;
+
+      if (utils::strmatch(force->angle_style,"^hybrid"))
+        error->all(FLERR,"Fix adapt does not support angle_style hybrid");
+
+      delete [] astyle;
+
     } else if (ad->which == KSPACE) {
       if (force->kspace == nullptr)
         error->all(FLERR,"Fix adapt kspace style does not exist");
       kspace_scale = (double *) force->kspace->extract("scale");
 
     } else if (ad->which == ATOM) {
-      if (ad->aparam == DIAMETER) {
+      if (ad->atomparam == DIAMETER) {
         if (!atom->radius_flag)
           error->all(FLERR,"Fix adapt requires atom attribute diameter");
         if (!atom->rmass_flag)
@@ -398,7 +455,7 @@ void FixAdapt::init()
           error->all(FLERR,"Fix adapt requires 2d simulation");
         if (!restart_reset) previous_diam_scale = 1.0;
       }
-      if (ad->aparam == CHARGE) {
+      if (ad->atomparam == CHARGE) {
         if (!atom->q_flag)
           error->all(FLERR,"Fix adapt requires atom attribute charge");
         if (!restart_reset) previous_chg_scale = 1.0;
@@ -408,7 +465,7 @@ void FixAdapt::init()
 
   if (restart_reset) restart_reset = 0;
 
-  // make copy of original pair/bond array values
+  // make copy of original pair/bond/angle array values
 
   for (int m = 0; m < nadapt; m++) {
     Adapt *ad = &adapt[m];
@@ -422,6 +479,10 @@ void FixAdapt::init()
     } else if (ad->which == BOND && ad->bdim == 1) {
       for (i = ad->ilo; i <= ad->ihi; ++i )
         ad->vector_orig[i] = ad->vector[i];
+
+    } else if (ad->which == ANGLE && ad->adim == 1) {
+      for (i = ad->ilo; i <= ad->ihi; ++i )
+        ad->vector_orig[i] = ad->vector[i];
     }
 
   }
@@ -483,7 +544,7 @@ void FixAdapt::post_run()
 }
 
 /* ----------------------------------------------------------------------
-   change pair,kspace,atom parameters based on variable evaluation
+   change pair,bond,angle,kspace,atom parameters based on variable evaluation
 ------------------------------------------------------------------------- */
 
 void FixAdapt::change_settings()
@@ -527,6 +588,18 @@ void FixAdapt::change_settings()
             ad->vector[i] = value;
       }
 
+    // set angle type array values:
+
+    } else if (ad->which == ANGLE) {
+      if (ad->adim == 1) {
+        if (scaleflag)
+          for (i = ad->ilo; i <= ad->ihi; ++i )
+            ad->vector[i] = value*ad->vector_orig[i];
+        else
+          for (i = ad->ilo; i <= ad->ihi; ++i )
+            ad->vector[i] = value;
+      }
+
     // set kspace scale factor
 
     } else if (ad->which == KSPACE) {
@@ -540,7 +613,7 @@ void FixAdapt::change_settings()
       // also reset rmass to new value assuming density remains constant
       // for scaleflag, previous_diam_scale is the scale factor on previous step
 
-      if (ad->aparam == DIAMETER) {
+      if (ad->atomparam == DIAMETER) {
         double scale;
         double *radius = atom->radius;
         double *rmass = atom->rmass;
@@ -567,7 +640,7 @@ void FixAdapt::change_settings()
       // reset charge to new value, for both owned and ghost atoms
       // for scaleflag, previous_chg_scale is the scale factor on previous step
 
-      } else if (ad->aparam == CHARGE) {
+      } else if (ad->atomparam == CHARGE) {
         double scale;
         double *q = atom->q;
         int *mask = atom->mask;
@@ -591,7 +664,7 @@ void FixAdapt::change_settings()
   modify->addstep_compute(update->ntimestep + nevery);
 
   // re-initialize pair styles if any PAIR settings were changed
-  // ditto for bond styles if any BOND settings were changed
+  // ditto for bond/angle styles if any BOND/ANGLE settings were changed
   // this resets other coeffs that may depend on changed values,
   //   and also offset and tail corrections
   // we must call force->pair->reinit() instead of the individual
@@ -601,6 +674,7 @@ void FixAdapt::change_settings()
 
   if (anypair) force->pair->reinit();
   if (anybond) force->bond->reinit();
+  if (anyangle) force->angle->reinit();
 
   // reset KSpace charges if charges have changed
 
@@ -624,7 +698,13 @@ void FixAdapt::restore_settings()
       }
 
     } else if (ad->which == BOND) {
-      if (ad->pdim == 1) {
+      if (ad->bdim == 1) {
+        for (int i = ad->ilo; i <= ad->ihi; i++)
+          ad->vector[i] = ad->vector_orig[i];
+      }
+
+    } else if (ad->which == ANGLE) {
+      if (ad->adim == 1) {
         for (int i = ad->ilo; i <= ad->ihi; i++)
           ad->vector[i] = ad->vector_orig[i];
       }

src/fix_adapt.h

@@ -45,7 +45,7 @@ class FixAdapt : public Fix {
 
  private:
   int nadapt, resetflag, scaleflag, massflag;
-  int anypair, anybond;
+  int anypair, anybond, anyangle;
   int nlevels_respa;
   char *id_fix_diam, *id_fix_chg;
   class FixStore *fix_diam, *fix_chg;
@@ -57,14 +57,16 @@ class FixAdapt : public Fix {
     char *var;
     char *pstyle, *pparam;
     char *bstyle, *bparam;
+    char *astyle, *aparam;
     int ilo, ihi, jlo, jhi;
-    int pdim, bdim;
+    int pdim, bdim, adim;
     double *scalar, scalar_orig;
     double *vector, *vector_orig;
     double **array, **array_orig;
-    int aparam;
+    int atomparam;
     class Pair *pair;
     class Bond *bond;
+    class Angle *angle;
   };
 
   Adapt *adapt;