Merge pull request #4333 from evoyiatzis/patch-4

add extract() function to several angle styles

doc/src/fix_adapt.rst

@@ -322,25 +322,29 @@ all types from 1 to :math:`N`.  A leading asterisk means all types from
 If :doc:`bond_style hybrid <bond_hybrid>` is used, *bstyle* should be a
 sub-style name. The bond styles that currently work with fix adapt are:
 
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`class2 <bond_class2>`                       | r0                        | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`fene <bond_fene>`                           | k,r0                      | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
+| :doc:`fene/expand <bond_fene_expand>`             | k,r0,epsilon,sigma,shift  | type bonds |
++---------------------------------------------------+---------------------------+------------+
 | :doc:`fene/nm <bond_fene>`                        | k,r0                      | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`gromos <bond_gromos>`                       | k,r0                      | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`harmonic <bond_harmonic>`                   | k,r0                      | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`harmonic/shift <bond_harmonic_shift>`       | k,r0,r1                   | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`harmonic/restrain <bond_harmonic_restrain>` | k                         | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
+| :doc:`mm3 <bond_mm3>`                             | k,r0                      | type bonds |
++---------------------------------------------------+---------------------------+------------+
 | :doc:`morse <bond_morse>`                         | r0                        | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 | :doc:`nonlinear <bond_nonlinear>`                 | epsilon,r0                | type bonds |
-+---------------------------------------------------+------------+------------+
++---------------------------------------------------+---------------------------+------------+
 
 ----------
 
@@ -363,11 +367,31 @@ all types from 1 to :math:`N`.  A leading asterisk means all types from
 If :doc:`angle_style hybrid <angle_hybrid>` is used, *astyle* should be a
 sub-style name. The angle styles that currently work with fix adapt are:
 
-+------------------------------------+----------+-------------+
++--------------------------------------------------------------------+-----------------+-------------+
 | :doc:`harmonic <angle_harmonic>`                                   | k,theta0        | type angles |
-+------------------------------------+----------+-------------+
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`charmm <angle_charmm>`                                       | k,theta0        | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`class2 <angle_class2>`                                       | k2,k3,k4,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
 | :doc:`cosine <angle_cosine>`                                       | k               | type angles |
-+------------------------------------+----------+-------------+
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`cosine/periodic <angle_cosine_periodic>`                     | k,b,n           | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`cosine/squared/restricted <angle_cosine_squared_restricted>` | k,theta0        | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`dipole <angle_dipole>`                                       | k,gamma0        | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`fourier <angle_fourier>`                                     | k,c0,c1,c2      | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`fourier/simple <angle_fourier_simple>`                       | k,c,n           | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`mm3 <angle_mm3>`                                             | k,theta0        | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`quartic <angle_quartic>`                                     | k2,k3,k4,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`spica <angle_spica>`                                         | 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/CG-SPICA/angle_spica.cpp

@@ -522,3 +522,15 @@ double AngleSPICA::single(int type, int i1, int i2, int i3)
   double tk = k[type] * dtheta;
   return tk*dtheta + e13;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleSPICA::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/CG-SPICA/angle_spica.h

@@ -37,6 +37,7 @@ class AngleSPICA : 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/CLASS2/angle_class2.cpp

@@ -467,3 +467,17 @@ double AngleClass2::single(int type, int i1, int i2, int i3)
 
   return energy;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleClass2::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k2") == 0) return (void *) k2;
+  if (strcmp(str, "k3") == 0) return (void *) k3;
+  if (strcmp(str, "k4") == 0) return (void *) k4;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/CLASS2/angle_class2.h

@@ -35,6 +35,7 @@ class AngleClass2 : 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 *theta0, *k2, *k3, *k4;

src/DIPOLE/angle_dipole.cpp

@@ -263,3 +263,15 @@ double AngleDipole::single(int type, int iRef, int iDip, int /*iDummy*/)
 
   return kdg * deltaGamma;    // energy
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleDipole::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "gamma0") == 0) return (void *) gamma0;
+  return nullptr;
+}

src/DIPOLE/angle_dipole.h

@@ -36,6 +36,7 @@ class AngleDipole : 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, *gamma0;

src/EXTRA-MOLECULE/angle_cosine_periodic.cpp

@@ -336,3 +336,16 @@ void AngleCosinePeriodic::born_matrix(int type, int i1, int i2, int i3, double &
   du = prefactor * sin(m_angle) / s;
   du2 = prefactor * (c * sin(m_angle) - s * cos(m_angle) * multiplicity[type]) / (s * s * s);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCosinePeriodic::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "b") == 0) return (void *) b;
+  if (strcmp(str, "multiplicity") == 0) return (void *) multiplicity;
+  return nullptr;
+}

src/EXTRA-MOLECULE/angle_cosine_periodic.h

@@ -36,6 +36,7 @@ class AngleCosinePeriodic : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k;

src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp

@@ -296,3 +296,15 @@ void AngleCosineSquaredRestricted::born_matrix(int type, int i1, int i2, int i3,
 
   du2 = 2 * k[type] * numerator / denominator;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCosineSquaredRestricted::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/EXTRA-MOLECULE/angle_cosine_squared_restricted.h

@@ -36,6 +36,7 @@ class AngleCosineSquaredRestricted : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *theta0;

src/EXTRA-MOLECULE/angle_fourier.cpp

@@ -309,3 +309,16 @@ void AngleFourier::born_matrix(int type, int i1, int i2, int i3, double &du, dou
   du = k[type] * (C1[type] + 4 * C2[type] * c);
 }
 
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleFourier::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "C0") == 0) return (void *) C0;
+  if (strcmp(str, "C1") == 0) return (void *) C1;
+  if (strcmp(str, "C2") == 0) return (void *) C2;
+  return nullptr;
+}

src/EXTRA-MOLECULE/angle_fourier.h

@@ -36,6 +36,7 @@ class AngleFourier : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *C0, *C1, *C2;

src/EXTRA-MOLECULE/angle_fourier_simple.cpp

@@ -316,3 +316,16 @@ void AngleFourierSimple::born_matrix(int type, int i1, int i2, int i3, double &d
   du2 = k[type] * C[type] * N[type] * (cos(theta) * sin(N[type] * theta)
                   - N[type] * sin(theta) * cos(N[type] * theta)) / pow(sin(theta),3);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleFourierSimple::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "C") == 0) return (void *) C;
+  if (strcmp(str, "N") == 0) return (void *) N;
+  return nullptr;
+}

src/EXTRA-MOLECULE/angle_fourier_simple.h

@@ -36,6 +36,7 @@ class AngleFourierSimple : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *C, *N;

src/EXTRA-MOLECULE/angle_quartic.cpp

@@ -325,3 +325,17 @@ void AngleQuartic::born_matrix(int type, int i1, int i2, int i3, double &du, dou
   du2 = (2.0 * k2[type] + 6.0 * k3[type] * dtheta + 12.0 * k4[type] * dtheta2) / (s*s) -
           (2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) * c / (s*s*s);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleQuartic::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k2") == 0) return (void *) k2;
+  if (strcmp(str, "k3") == 0) return (void *) k3;
+  if (strcmp(str, "k4") == 0) return (void *) k4;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/EXTRA-MOLECULE/angle_quartic.h

@@ -36,6 +36,7 @@ class AngleQuartic : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k2, *k3, *k4, *theta0;

src/MOLECULE/angle_charmm.cpp

@@ -309,3 +309,15 @@ double AngleCharmm::single(int type, int i1, int i2, int i3)
 
   return (tk * dtheta + rk * dr);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCharmm::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_charmm.h

@@ -35,6 +35,7 @@ class AngleCharmm : 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, *k_ub, *r_ub;

src/MOLECULE/bond_fene_expand.cpp

@@ -273,3 +273,18 @@ double BondFENEExpand::single(int type, double rsq, int /*i*/, int /*j*/, double
 
   return eng;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *BondFENEExpand::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "r0") == 0) return (void *) r0;
+  if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
+  if (strcmp(str, "sigma") == 0) return (void *) sigma;
+  if (strcmp(str, "shift") == 0) return (void *) shift;
+  return nullptr;
+}

src/MOLECULE/bond_fene_expand.h

@@ -36,6 +36,7 @@ class BondFENEExpand : public Bond {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   double single(int, double, int, int, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *k, *r0, *epsilon, *sigma, *shift;

src/YAFF/angle_mm3.cpp

@@ -327,3 +327,15 @@ void AngleMM3::born_matrix(int type, int i1, int i2, int i3, double &du, double
   du = -k2[type] * df / s;
   du2 = k2[type] * (d2f - df  * c / s) / (s * s) ;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleMM3::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k2") == 0) return (void *) k2;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/YAFF/angle_mm3.h

@@ -36,6 +36,7 @@ class AngleMM3 : public Angle {
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
   void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *theta0, *k2;

src/YAFF/bond_mm3.cpp

@@ -238,3 +238,15 @@ void BondMM3::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du
   du = 2.0 * k2[type] * dr + 3.0 * K3 * dr2 + 4.0 * K4 * dr3;
   du2 = 2.0 * k2[type] + 6.0 * K3 * dr + 12.0 * K4 * dr2;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *BondMM3::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k2") == 0) return (void *) k2;
+  if (strcmp(str, "r0") == 0) return (void *) r0;
+  return nullptr;
+}

src/YAFF/bond_mm3.h

@@ -36,6 +36,7 @@ class BondMM3 : public Bond {
   void write_data(FILE *) override;
   double single(int, double, int, int, double &) override;
   void born_matrix(int, double, int, int, double &, double &) override;
+  void *extract(const char *, int &) override;
 
  protected:
   double *r0, *k2;

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

@@ -15,7 +15,9 @@ angle_coeff: ! |
   3  40.0 120.0  35.0 2.410
   4  33.0 108.5  30.0 2.163
 equilibrium: 4 1.9216075064457567 1.9425514574696887 2.0943951023931953 1.8936822384138476
-extract: ! ""
+extract: ! |
+  k 1
+  theta0 1
 natoms: 29
 init_energy: 85.42486388459771
 init_stress: ! |2-

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

@@ -23,7 +23,11 @@ angle_coeff: ! |
   3  ba  10.0 10.0 1.5 1.5
   4  ba   0.0 20.0 1.5 1.5
 equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k2 1
+  k3 1
+  k4 1
+  theta0 1
 natoms: 29
 init_energy: 46.44089683774903
 init_stress: ! |2-

unittest/force-styles/tests/angle-cosine_periodic.yaml

@@ -15,7 +15,10 @@ angle_coeff: ! |
   3  50.0 -1 3
   4 100.0 -1 4
 equilibrium: 4 3.141592653589793 3.141592653589793 2.0943951023931957 2.356194490192345
-extract: ! ""
+extract: ! |
+  k 1
+  b 1
+  multiplicity 1
 natoms: 29
 init_energy: 1178.5476942873006
 init_stress: ! |2-

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

@@ -17,7 +17,9 @@ angle_coeff: ! |
   3  50.0 120.0
   4 100.0 108.5
 equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476
-extract: ! ""
+extract: ! |
+  k 1
+  theta0 1
 natoms: 29
 init_energy: 43.16721849625078
 init_stress: ! |2-

unittest/force-styles/tests/angle-dipole.yaml

@@ -20,7 +20,9 @@ angle_coeff: ! |
   3  50.0 120.0
   4 100.0 108.5
 equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k 1
+  gamma0 1
 natoms: 29
 init_energy: 1003.6681304854917
 init_stress: ! |2-

unittest/force-styles/tests/angle-fourier.yaml

@@ -16,7 +16,11 @@ angle_coeff: ! |
   3  50.0 0.0 0.0 1.0
   4 100.0 0.3 0.3 0.3
 equilibrium: 4 3.141592653589793 1.5707963267948966 1.5707963267948966 1.8234765819369754
-extract: ! ""
+extract: ! |
+  k 1
+  C0 1
+  C1 1
+  C2 1
 natoms: 29
 init_energy: 400.84036632010225
 init_stress: ! |-

unittest/force-styles/tests/angle-fourier_simple.yaml

@@ -16,7 +16,10 @@ angle_coeff: ! |
   3  50.0  1.0 3.0
   4 100.0 -0.5 1.5
 equilibrium: 4 3.141592653589793 1.5707963267948966 1.0471975511965976 2.0943951023931953
-extract: ! ""
+extract: ! |
+  k 1
+  C 1
+  N 1
 natoms: 29
 init_energy: 2474.0748013590646
 init_stress: ! |-

unittest/force-styles/tests/angle-mm3.yaml

@@ -16,7 +16,9 @@ angle_coeff: ! |
   3  50.0 120.0
   4 100.0 108.5
 equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k2 1
+  theta0 1
 natoms: 29
 init_energy: 44.72461548562619
 init_stress: ! |2-

unittest/force-styles/tests/angle-quartic.yaml

@@ -16,7 +16,11 @@ angle_coeff: ! |
   3 120.0  50.0  -9.5 -1.5
   4 108.5 100.0   5.0 -2.0
 equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k2 1
+  k3 1
+  k4 1
+  theta0 1
 natoms: 29
 init_energy: 41.0458477552901
 init_stress: ! |2-

unittest/force-styles/tests/angle-spica.yaml

@@ -20,7 +20,9 @@ angle_coeff: ! |
   3  40.0 120.0
   4  33.0 108.5
 equilibrium: 4 1.9216075064457565 1.9425514574696887 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k 1
+  theta0 1
 natoms: 29
 init_energy: 38.36438529349082
 init_stress: ! |2-

unittest/force-styles/tests/bond-fene_expand.yaml

@@ -17,7 +17,12 @@ bond_coeff: ! |
   4 650 2.4 0.015 1.2 0.15
   5 450 2   0.018 1   0.09
 equilibrium: 5 1.5550000000000002 1.117 1.321 1.3139999999999998 1.06
-extract: ! ""
+extract: ! |
+  k 1
+  r0 1
+  epsilon 1
+  sigma 1
+  shift 1
 natoms: 29
 init_energy: 5926.020859124294
 init_stress: ! |-

unittest/force-styles/tests/bond-mm3.yaml

@@ -17,7 +17,9 @@ bond_coeff: ! |
   4 650.0 1.2
   5 450.0 1.0
 equilibrium: 5 1.5 1.1 1.3 1.2 1
-extract: ! ""
+extract: ! |
+  k2 1
+  r0 1
 natoms: 29
 init_energy: 4.247265008273143
 init_stress: ! |-