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

extend fix adapt to treat dihedrals

Author: akohlmey

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;
+}