Merge pull request #4540 from evoyiatzis/fix-adapt-dihedral

extend fix adapt to treat dihedrals

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 *angle* or *improper* or *kspace* or *atom*
+* attribute = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *atom*
 
   .. parsed-literal::
 
@@ -33,6 +33,11 @@ Syntax
          aparam = parameter to adapt over time
          I = type angle to set parameter for (integer or type label)
          v_name = variable with name that calculates value of aparam
+       *dihedral* args = dstyle dparam I v_name
+         dstyle = dihedral style name (e.g., quadratic)
+         dparam = parameter to adapt over time
+         I = type dihedral to set parameter for (integer or type label)
+         v_name = variable with name that calculates value of iparam
        *improper* args = istyle iparam I v_name
          istyle = improper style name (e.g., cvff)
          iparam = parameter to adapt over time
@@ -433,6 +438,48 @@ this fix uses to reset theta0 needs to generate values in radians.
 
 ----------
 
+.. versionadded:: TBD
+
+The *dihedral* keyword uses the specified variable to change the value of
+a dihedral coefficient over time, very similar to how the *angle* keyword
+operates. The only difference is that now a dihedral coefficient for a
+given dihedral type is adapted.
+
+A wild-card asterisk can be used in place of or in conjunction with the
+dihedral type argument to set the coefficients for multiple dihedral types.
+This takes the form "\*" or "\*n" or "m\*" or "m\*n".  If :math:`N` is
+the number of dihedral types, then an asterisk with no numeric values means
+all types from 1 to :math:`N`.  A leading asterisk means all types from
+1 to n (inclusive).  A trailing asterisk means all types from m to
+:math:`N` (inclusive).  A middle asterisk means all types from m to n
+(inclusive).
+
+If :doc:`dihedral_style hybrid <dihedral_hybrid>` is used, *dstyle* should be a
+sub-style name. The dihedral styles that currently work with fix adapt are:
+
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`charmm  <dihedral_charmm>`                                       | k,n,d                   | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`charmmfsw <dihedral_charmm>`                                     | k,n,d                   | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`class2 <dihedral_class2>`                                        | k1,k2,k3,phi1,phi2,phi3 | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`cosine/squared/restricted <dihedral_cosine_squared_restricted>`  | k,phi0                  | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`helix <dihedral_helix>`                                          | a,b,c                   | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`multi/harmonic <dihedral_multi_harmonic>`                        | a1,a2,a3,a4,a5          | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`opls <dihedral_opls>`                                            | k1,k2,k3,k4             | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+| :doc:`quadratic <dihedral_quadratic>`                                  | k,phi0                  | type dihedrals |
++------------------------------------------------------------------------+-------------------------+----------------+
+
+Note that internally, phi0 is stored in radians, so the variable
+this fix use to reset phi0 needs to generate values in radians.
+
+----------
+
 .. versionadded:: 2Apr2025
 
 The *improper* keyword uses the specified variable to change the value of

src/CLASS2/dihedral_class2.cpp

@@ -946,3 +946,18 @@ void DihedralClass2::write_data(FILE *fp)
             at_theta0_1[i]*180.0/MY_PI,at_theta0_2[i]*180.0/MY_PI);
 }
 
+/* ----------------------------------------------------------------------
+    return ptr to internal members upon request
+ ------------------------------------------------------------------------ */
+
+ void *DihedralClass2::extract(const char *str, int &dim)
+ {
+   dim = 1;
+   if (strcmp(str, "k1") == 0) return (void *) k1;
+   if (strcmp(str, "k2") == 0) return (void *) k2;
+   if (strcmp(str, "k3") == 0) return (void *) k3;
+   if (strcmp(str, "phi1") == 0) return (void *) phi1;
+   if (strcmp(str, "phi2") == 0) return (void *) phi2;
+   if (strcmp(str, "phi3") == 0) return (void *) phi3;
+   return nullptr;
+ }

src/CLASS2/dihedral_class2.h

@@ -33,6 +33,7 @@ class DihedralClass2 : public Dihedral {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k1, *k2, *k3;

src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.cpp

@@ -393,3 +393,15 @@ void DihedralCosineSquaredRestricted::born_matrix(int nd, int i1, int i2, int i3
 
   du2 = 2 * k[type] * numerator / denominator;
 }
+
+/* ----------------------------------------------------------------------
+    return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralCosineSquaredRestricted::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "phi0") == 0) return (void *) phi0;
+  return nullptr;
+}

src/EXTRA-MOLECULE/dihedral_cosine_squared_restricted.h

@@ -34,6 +34,7 @@ class DihedralCosineSquaredRestricted : public Dihedral {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   void born_matrix(int, int, int, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *phi0;

src/EXTRA-MOLECULE/dihedral_helix.cpp

@@ -430,3 +430,16 @@ void DihedralHelix::born_matrix(int nd, int i1, int i2, int i3, int i4,
   du2 = -(9.0*bphi[type]*cos(3.0*phi) + cphi[type]*cos(phi + MY_PI4))*siinv*siinv +
           (3.0*bphi[type]*sin(3.0*phi) + cphi[type]*sin(phi + MY_PI4))*c*siinv*siinv*siinv;
 }
+
+/* ----------------------------------------------------------------------
+    return ptr to internal members upon request
+ ------------------------------------------------------------------------ */
+
+ void *DihedralHelix::extract(const char *str, int &dim)
+ {
+   dim = 1;
+   if (strcmp(str, "a") == 0) return (void *) aphi;
+   if (strcmp(str, "b") == 0) return (void *) bphi;
+   if (strcmp(str, "c") == 0) return (void *) cphi;
+   return nullptr;
+ }

src/EXTRA-MOLECULE/dihedral_helix.h

@@ -34,6 +34,7 @@ class DihedralHelix : public Dihedral {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   void born_matrix(int, int, int, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *aphi, *bphi, *cphi;

src/EXTRA-MOLECULE/dihedral_quadratic.cpp

@@ -435,3 +435,15 @@ void DihedralQuadratic::born_matrix(int nd, int i1, int i2, int i3, int i4,
   du = - 2.0 * k[type] * dphi * siinv;
   du2 = 2.0 * k[type] * siinv * siinv * ( 1.0 - dphi * c * siinv) ;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralQuadratic::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "phi0") == 0) return (void *) phi0;
+  return nullptr;
+}

src/EXTRA-MOLECULE/dihedral_quadratic.h

@@ -34,6 +34,7 @@ class DihedralQuadratic : public Dihedral {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   void born_matrix(int, int, int, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *phi0;

src/MOLECULE/dihedral_charmm.cpp

@@ -432,3 +432,16 @@ void DihedralCharmm::write_data(FILE *fp)
   for (int i = 1; i <= atom->ndihedraltypes; i++)
     fprintf(fp, "%d %g %d %d %g\n", i, k[i], multiplicity[i], shift[i], weight[i]);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralCharmm::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "n") == 0) return (void *) multiplicity;
+  if (strcmp(str, "d") == 0) return (void *) shift;
+  return nullptr;
+}

src/MOLECULE/dihedral_charmm.h

@@ -34,6 +34,7 @@ class DihedralCharmm : public Dihedral {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *weight, *cos_shift, *sin_shift;

src/MOLECULE/dihedral_charmmfsw.cpp

@@ -475,3 +475,16 @@ void DihedralCharmmfsw::write_data(FILE *fp)
   for (int i = 1; i <= atom->ndihedraltypes; i++)
     fprintf(fp, "%d %g %d %d %g\n", i, k[i], multiplicity[i], shift[i], weight[i]);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralCharmmfsw::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "n") == 0) return (void *) multiplicity;
+  if (strcmp(str, "d") == 0) return (void *) shift;
+  return nullptr;
+}

src/MOLECULE/dihedral_charmmfsw.h

@@ -34,6 +34,7 @@ class DihedralCharmmfsw : public Dihedral {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
+  void *extract(const char *, int &) override;
 
  protected:
   int implicit, weightflag, dihedflag;

src/MOLECULE/dihedral_multi_harmonic.cpp

@@ -407,3 +407,19 @@ void DihedralMultiHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4,
   du = a2[type] + c * (2.0 * a3[type] + c * (3.0 * a4[type] + c * 4.0 * a5[type]));
   du2 = 2.0 * a3[type] + 6.0 * c * (a4[type] + 2.0 * a5[type] * c);
 }
+
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralMultiHarmonic::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "a1") == 0) return (void *) a1;
+  if (strcmp(str, "a2") == 0) return (void *) a2;
+  if (strcmp(str, "a3") == 0) return (void *) a3;
+  if (strcmp(str, "a4") == 0) return (void *) a4;
+  if (strcmp(str, "a5") == 0) return (void *) a5;
+  return nullptr;
+}

src/MOLECULE/dihedral_multi_harmonic.h

@@ -34,6 +34,7 @@ class DihedralMultiHarmonic : public Dihedral {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   void born_matrix(int, int, int, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *a1, *a2, *a3, *a4, *a5;

src/MOLECULE/dihedral_opls.cpp

@@ -431,3 +431,17 @@ void DihedralOPLS::born_matrix(int nd, int i1, int i2, int i3, int i4, double &d
          16.0 * k4[type] * si * cos(4.0 * phi) - 4.0 * k4[type] * sin(4.0 * phi)) /
       (si * si * si);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *DihedralOPLS::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k1") == 0) return (void *) k1;
+  if (strcmp(str, "k2") == 0) return (void *) k2;
+  if (strcmp(str, "k3") == 0) return (void *) k3;
+  if (strcmp(str, "k4") == 0) return (void *) k4;
+  return nullptr;
+}

src/MOLECULE/dihedral_opls.h

@@ -34,6 +34,7 @@ class DihedralOPLS : public Dihedral {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   void born_matrix(int, int, int, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k1, *k2, *k3, *k4;

src/dihedral.cpp

@@ -33,6 +33,7 @@ Dihedral::Dihedral(LAMMPS *_lmp) : Pointers(_lmp)
 {
   energy = 0.0;
   writedata = 0;
+  reinitflag = 1;
 
   allocated = 0;
   suffix_flag = Suffix::NONE;
@@ -428,3 +429,15 @@ double Dihedral::memory_usage()
   bytes += (double) comm->nthreads * maxcvatom * 9 * sizeof(double);
   return bytes;
 }
+
+/* -----------------------------------------------------------------------
+   reset all type-based dihedral params via init()
+-------------------------------------------------------------------------- */
+
+void Dihedral::reinit()
+{
+  if (!reinitflag)
+    error->all(FLERR, "Fix adapt interface to this dihedral style not supported");
+
+  init();
+}

src/dihedral.h

@@ -37,6 +37,9 @@ class Dihedral : 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;
@@ -62,6 +65,8 @@ class Dihedral : protected Pointers {
     du = 0.0;
     du2 = 0.0;
   }
+  virtual void *extract(const char *, int &) { return nullptr; }
+  void reinit();
 
  protected:
   int suffix_flag;    // suffix compatibility flag

src/dihedral_hybrid.cpp

@@ -326,6 +326,14 @@ void DihedralHybrid::init_style()
     if (styles[m]) styles[m]->init_style();
 }
 
+/* ---------------------------------------------------------------------- */
+
+int DihedralHybrid::check_itype(int itype, char *substyle)
+{
+  if (strcmp(keywords[map[itype]], substyle) == 0) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    proc 0 writes to restart file
 ------------------------------------------------------------------------- */

src/dihedral_hybrid.h

@@ -40,6 +40,8 @@ class DihedralHybrid : public Dihedral {
   void read_restart(FILE *) override;
   double memory_usage() override;
 
+  int check_itype(int, char *);
+
  protected:
   int *map;    // which style each dihedral type points to
 

src/fix_adapt.cpp

@@ -19,6 +19,8 @@
 #include "atom.h"
 #include "bond.h"
 #include "bond_hybrid.h"
+#include "dihedral.h"
+#include "dihedral_hybrid.h"
 #include "domain.h"
 #include "error.h"
 #include "fix_store_atom.h"
@@ -43,14 +45,14 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{PAIR, KSPACE, ATOM, BOND, ANGLE, IMPROPER};
+enum{PAIR, KSPACE, ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER};
 enum{DIAMETER, CHARGE};
 
 /* ---------------------------------------------------------------------- */
 
 FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0),
-  anyimproper(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
+  anydihedral(0), anyimproper(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr)
 {
   if (narg < 5) utils::missing_cmd_args(FLERR,"fix adapt", error);
   nevery = utils::inumeric(FLERR,arg[3],false,lmp);
@@ -85,6 +87,10 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
       if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt angle", error);
       nadapt++;
       iarg += 5;
+    } else if (strcmp(arg[iarg],"dihedral") == 0) {
+      if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt dihedral", error);
+      nadapt++;
+      iarg += 5;
     } else if (strcmp(arg[iarg],"improper") == 0) {
       if (iarg+5 > narg) utils::missing_cmd_args(FLERR,"fix adapt improper", error);
       nadapt++;
@@ -155,6 +161,19 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
       nadapt++;
       iarg += 5;
 
+    } else if (strcmp(arg[iarg],"dihedral") == 0) {
+      adapt[nadapt].which = DIHEDRAL;
+      adapt[nadapt].dihedral = nullptr;
+      adapt[nadapt].dstyle = utils::strdup(arg[iarg+1]);
+      adapt[nadapt].dparam = utils::strdup(arg[iarg+2]);
+      utils::bounds_typelabel(FLERR, arg[iarg+3], 1, atom->ndihedraltypes,
+                              adapt[nadapt].ilo, adapt[nadapt].ihi, lmp, Atom::DIHEDRAL);
+      if (utils::strmatch(arg[iarg+4],"^v_")) {
+        adapt[nadapt].var = utils::strdup(arg[iarg+4]+2);
+      } else error->all(FLERR,"Argument #{} must be variable not {}", iarg+5, arg[iarg+4]);
+      nadapt++;
+      iarg += 5;
+
     } else if (strcmp(arg[iarg],"improper") == 0) {
       adapt[nadapt].which = IMPROPER;
       adapt[nadapt].improper = nullptr;
@@ -246,6 +265,12 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) :
 
   // allocate improper style arrays:
 
+  n = atom->ndihedraltypes;
+  for (int m = 0; m < nadapt; ++m)
+    if (adapt[m].which == DIHEDRAL) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
+
+  // allocate improper style arrays:
+
   n = atom->nimpropertypes;
   for (int m = 0; m < nadapt; ++m)
     if (adapt[m].which == IMPROPER) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
@@ -270,6 +295,10 @@ FixAdapt::~FixAdapt()
       delete[] adapt[m].astyle;
       delete[] adapt[m].aparam;
       memory->destroy(adapt[m].vector_orig);
+    } else if (adapt[m].which == DIHEDRAL) {
+      delete[] adapt[m].dstyle;
+      delete[] adapt[m].dparam;
+      memory->destroy(adapt[m].vector_orig);
     } else if (adapt[m].which == IMPROPER) {
       delete[] adapt[m].istyle;
       delete[] adapt[m].iparam;
@@ -368,6 +397,7 @@ void FixAdapt::init()
   anypair = 0;
   anybond = 0;
   anyangle = 0;
+  anydihedral = 0;
   anyimproper = 0;
 
   for (int m = 0; m < nadapt; m++) {
@@ -509,6 +539,40 @@ void FixAdapt::init()
 
       delete[] astyle;
 
+    } else if (ad->which == DIHEDRAL) {
+      ad->dihedral = nullptr;
+      anydihedral = 1;
+
+      char *dstyle = utils::strdup(ad->dstyle);
+      if (lmp->suffix_enable)
+        ad->dihedral = force->dihedral_match(fmt::format("{}/{}",dstyle,lmp->suffix));
+
+      if (ad->dihedral == nullptr) ad->dihedral = force->dihedral_match(dstyle);
+      if (ad->dihedral == nullptr )
+        error->all(FLERR,"Fix adapt dihedral style {} does not exist", dstyle);
+
+      void *ptr = ad->dihedral->extract(ad->dparam,ad->ddim);
+
+      if (ptr == nullptr)
+        error->all(FLERR,"Fix adapt dihedral style parameter {} not supported", ad->dparam);
+
+      // for dihedral styles, use a vector
+
+      if (ad->ddim == 1) ad->vector = (double *) ptr;
+
+      if (utils::strmatch(force->dihedral_style,"^hybrid")) {
+        auto dihedral = dynamic_cast<DihedralHybrid *>(force->dihedral);
+        if (dihedral) {
+          for (i = ad->ilo; i <= ad->ihi; i++) {
+            if (!dihedral->check_itype(i,dstyle))
+              error->all(FLERR,"Fix adapt type dihedral range is not valid "
+                         "for dihedral hybrid sub-style {}", dstyle);
+          }
+        }
+      }
+
+      delete[] dstyle;
+
     } else if (ad->which == IMPROPER) {
       ad->improper = nullptr;
       anyimproper = 1;
@@ -569,7 +633,7 @@ void FixAdapt::init()
 
   if (restart_reset) restart_reset = 0;
 
-  // make copy of original pair/bond/angle/improper array values
+  // make copy of original pair/bond/angle/dihedral/improper array values
 
   for (int m = 0; m < nadapt; m++) {
     Adapt *ad = &adapt[m];
@@ -588,6 +652,10 @@ void FixAdapt::init()
       for (i = ad->ilo; i <= ad->ihi; ++i )
         ad->vector_orig[i] = ad->vector[i];
 
+    } else if (ad->which == DIHEDRAL && ad->ddim == 1) {
+      for (i = ad->ilo; i <= ad->ihi; ++i )
+        ad->vector_orig[i] = ad->vector[i];
+
     } else if (ad->which == IMPROPER && ad->idim == 1) {
       for (i = ad->ilo; i <= ad->ihi; ++i )
         ad->vector_orig[i] = ad->vector[i];
@@ -710,6 +778,18 @@ void FixAdapt::change_settings()
             ad->vector[i] = value;
       }
 
+    // set dihedral type array values:
+
+  } else if (ad->which == DIHEDRAL) {
+    if (ad->ddim == 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 improper type array values:
 
     } else if (ad->which == IMPROPER) {
@@ -797,6 +877,7 @@ void FixAdapt::change_settings()
   if (anypair) force->pair->reinit();
   if (anybond) force->bond->reinit();
   if (anyangle) force->angle->reinit();
+  if (anydihedral) force->dihedral->reinit();
   if (anyimproper) force->improper->reinit();
 
   // reset KSpace charges if charges have changed
@@ -832,6 +913,12 @@ void FixAdapt::restore_settings()
           ad->vector[i] = ad->vector_orig[i];
       }
 
+    } else if (ad->which == DIHEDRAL) {
+      if (ad->ddim == 1) {
+        for (int i = ad->ilo; i <= ad->ihi; i++)
+          ad->vector[i] = ad->vector_orig[i];
+      }
+
     } else if (ad->which == IMPROPER) {
       if (ad->idim == 1) {
         for (int i = ad->ilo; i <= ad->ihi; i++)
@@ -878,6 +965,7 @@ void FixAdapt::restore_settings()
   if (anypair) force->pair->reinit();
   if (anybond) force->bond->reinit();
   if (anyangle) force->angle->reinit();
+  if (anydihedral) force->dihedral->reinit();
   if (anyimproper) force->improper->reinit();
   if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }

src/fix_adapt.h

@@ -44,7 +44,7 @@ class FixAdapt : public Fix {
 
  private:
   int nadapt, resetflag, scaleflag, massflag;
-  int anypair, anybond, anyangle, anyimproper;
+  int anypair, anybond, anyangle, anydihedral, anyimproper;
   int nlevels_respa;
   char *id_fix_diam, *id_fix_chg;
   class FixStoreAtom *fix_diam, *fix_chg;
@@ -57,9 +57,10 @@ class FixAdapt : public Fix {
     char *pstyle, *pparam;
     char *bstyle, *bparam;
     char *astyle, *aparam;
+    char *dstyle, *dparam;
     char *istyle, *iparam;
     int ilo, ihi, jlo, jhi;
-    int pdim, bdim, adim, idim;
+    int pdim, bdim, adim, ddim, idim;
     double *scalar, scalar_orig;
     double *vector, *vector_orig;
     double **array, **array_orig;
@@ -67,6 +68,7 @@ class FixAdapt : public Fix {
     class Pair *pair;
     class Bond *bond;
     class Angle *angle;
+    class Dihedral *dihedral;
     class Improper *improper;
   };
 

unittest/force-styles/test_dihedral_style.cpp

@@ -746,3 +746,44 @@ TEST(DihedralStyle, numdiff)
     cleanup_lammps(lmp, test_config);
     if (!verbose) ::testing::internal::GetCapturedStdout();
 }
+
+TEST(DihedralStyle, extract)
+{
+    if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
+
+    LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"};
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    if (!lmp) {
+        std::cerr << "One or more prerequisite styles are not available "
+                     "in this LAMMPS configuration:\n";
+        for (auto prerequisite : test_config.prerequisites) {
+            std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+        }
+        GTEST_SKIP();
+    }
+
+    auto *dihedral = lmp->force->dihedral;
+    void *ptr   = nullptr;
+    int dim     = 0;
+    for (auto extract : test_config.extract) {
+        ptr = dihedral->extract(extract.first.c_str(), dim);
+        EXPECT_NE(ptr, nullptr);
+        EXPECT_EQ(dim, extract.second);
+    }
+    ptr = dihedral->extract("does_not_exist", dim);
+    EXPECT_EQ(ptr, nullptr);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    cleanup_lammps(lmp, test_config);
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+}

unittest/force-styles/tests/dihedral-charmm.yaml

@@ -19,7 +19,10 @@ dihedral_coeff: ! |
   3  56.0 0 110 0.0
   4  23.0 1 180 0.5
   5  19.0 3  90 1.0
-extract: ! ""
+extract: ! |
+  k 1
+  n 1
+  d 1
 natoms: 29
 init_energy: 1317.959844120986
 init_stress: ! |2-

unittest/force-styles/tests/dihedral-charmmfsw.yaml

@@ -19,7 +19,10 @@ dihedral_coeff: ! |
   3  56.0 0 110 0.0
   4  23.0 1 180 0.5
   5  19.0 3  90 1.0
-extract: ! ""
+extract: ! |
+  k 1
+  n 1
+  d 1
 natoms: 29
 init_energy: 1317.959844120986
 init_stress: ! |2-

unittest/force-styles/tests/dihedral-class2.yaml

@@ -17,7 +17,13 @@ dihedral_coeff: ! |
   *  at   75 42 31 75 42 31 120 50
   *  aat  75 120 160
   *  bb13 75 1.4 1.4
-extract: ! ""
+extract: ! |
+  k1 1
+  k2 1
+  k3 1
+  phi1 1
+  phi2 1
+  phi3 1
 natoms: 29
 init_energy: 3355.0074717375933
 init_stress: ! |2-

unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml

@@ -16,7 +16,9 @@ dihedral_coeff: ! |
   3  15.0 -10.0
   4  12.0 0.0
   5  11.0 45.0
-extract: ! ""
+extract: ! |
+  k 1
+  phi0 1
 natoms: 29
 init_energy: 10643.96352037142
 init_stress: ! |2-

unittest/force-styles/tests/dihedral-helix.yaml

@@ -16,7 +16,10 @@ dihedral_coeff: ! |
   3  56 10 32
   4  23 80 61
   5  19 90 13
-extract: ! ""
+extract: ! |
+  a 1
+  b 1
+  c 1
 natoms: 29
 init_energy: 4634.436545707672
 init_stress: ! |-

unittest/force-styles/tests/dihedral-multi_harmonic.yaml

@@ -16,7 +16,12 @@ dihedral_coeff: ! |
   3  56 10 32 84 52
   4  23 80 61 83 74
   5  19 90 13 15 58
-extract: ! ""
+extract: ! |
+  a1 1
+  a2 1
+  a3 1
+  a4 1
+  a5 1
 natoms: 29
 init_energy: 2854.9857316566695
 init_stress: ! |2-

unittest/force-styles/tests/dihedral-opls.yaml

@@ -16,7 +16,11 @@ dihedral_coeff: ! |
   3  56 10 32 84
   4  23 80 61 83
   5  19 90 13 15
-extract: ! ""
+extract: ! |
+  k1 1
+  k2 1
+  k3 1
+  k4 1
 natoms: 29
 init_energy: 2260.6834525285753
 init_stress: ! |-

unittest/force-styles/tests/dihedral-quadratic.yaml

@@ -17,7 +17,9 @@ dihedral_coeff: ! |
   3  25 100
   4  35 110
   5  85 120
-extract: ! ""
+extract: ! |
+  k 1
+  phi0 1
 natoms: 29
 init_energy: 6216.314347066598
 init_stress: ! |2-