Merge pull request #4356 from evoyiatzis/patch-1

add extract() function to few angle and bond styles

doc/src/fix_adapt.rst

@@ -323,7 +323,7 @@ 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:`class2 <bond_class2>`                         | k2,k3,k4,r0               | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
 | :doc:`fene <bond_fene>`                             | k,r0                      | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
@@ -331,6 +331,8 @@ sub-style name. The bond styles that currently work with fix adapt are:
 +-----------------------------------------------------+---------------------------+------------+
 | :doc:`fene/nm <bond_fene>`                          | k,r0                      | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
+| :doc:`gaussian <bond_gaussian>`                     | alpha,width,r0            | type bonds |
++-----------------------------------------------------+---------------------------+------------+
 | :doc:`gromos <bond_gromos>`                         | k,r0                      | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
 | :doc:`harmonic <bond_harmonic>`                     | k,r0                      | type bonds |
@@ -343,9 +345,9 @@ sub-style name. The bond styles that currently work with fix adapt are:
 +-----------------------------------------------------+---------------------------+------------+
 | :doc:`mm3 <bond_mm3>`                               | k,r0                      | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
-| :doc:`morse <bond_morse>`                           | r0                        | type bonds |
+| :doc:`morse <bond_morse>`                           | d0,alpha,r0               | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
-| :doc:`nonlinear <bond_nonlinear>`                   | epsilon,r0                | type bonds |
+| :doc:`nonlinear <bond_nonlinear>`                   | lamda,epsilon,r0          | type bonds |
 +-----------------------------------------------------+---------------------------+------------+
 
 ----------
@@ -369,31 +371,37 @@ 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:`harmonic <angle_harmonic>`                                   | k,theta0           | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`charmm <angle_charmm>`                                       | k,theta0        | type angles |
+| :doc:`charmm <angle_charmm>`                                       | k,theta0           | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`class2 <angle_class2>`                                       | k2,k3,k4,theta0 | type angles |
+| :doc:`class2 <angle_class2>`                                       | k2,k3,k4,theta0    | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`cosine <angle_cosine>`                                       | k               | type angles |
+| :doc:`cosine <angle_cosine>`                                       | k                  | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`cosine/periodic <angle_cosine_periodic>`                     | k,b,n           | type angles |
+| :doc:`cosine/delta <angle_cosine_delta>`                           | k                  | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`cosine/squared/restricted <angle_cosine_squared_restricted>` | k,theta0        | type angles |
+| :doc:`cosine/periodic <angle_cosine_periodic>`                     | k,b,n              | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`dipole <angle_dipole>`                                       | k,gamma0        | type angles |
+| :doc:`cosine/squared <angle_cosine_squared>`                       | k,theta0           | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`fourier <angle_fourier>`                                     | k,c0,c1,c2      | type angles |
+| :doc:`cosine/squared/restricted <angle_cosine_squared_restricted>` | k,theta0           | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`fourier/simple <angle_fourier_simple>`                       | k,c,n           | type angles |
+| :doc:`dipole <angle_dipole>`                                       | k,gamma0           | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`mm3 <angle_mm3>`                                             | k,theta0        | type angles |
+| :doc:`fourier <angle_fourier>`                                     | k,c0,c1,c2         | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`quartic <angle_quartic>`                                     | k2,k3,k4,theta0 | type angles |
+| :doc:`fourier/simple <angle_fourier_simple>`                       | k,c,n              | type angles |
-+--------------------------------------------------------------------+-----------------+-------------+
++--------------------------------------------------------------------+--------------------+-------------+
-| :doc:`spica <angle_spica>`                                         | k,theta0        | type angles |
+| :doc:`gaussian <angle_gaussian>`                                   | alpha,width,theta0 | 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/AMOEBA/angle_amoeba.cpp

@@ -868,3 +868,19 @@ double AngleAmoeba::single(int type, int i1, int i2, int i3)
 
   return energy;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleAmoeba::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, "k5") == 0) return (void *) k5;
+  if (strcmp(str, "k6") == 0) return (void *) k6;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/AMOEBA/angle_amoeba.h

@@ -36,6 +36,7 @@ class AngleAmoeba : 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:
   int *pflag, *ubflag;

src/CLASS2/bond_class2.cpp

@@ -242,6 +242,9 @@ void BondClass2::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double
 void *BondClass2::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,"r0")==0) return (void*) r0;
   return nullptr;
 }

src/EXTRA-MOLECULE/angle_gaussian.cpp

@@ -352,3 +352,14 @@ double AngleGaussian::single(int type, int i1, int i2, int i3)
   if (sum_g_i < SMALL) sum_g_i = SMALL;
   return -(force->boltz * angle_temperature[type]) * log(sum_g_i);
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *AngleGaussian::extract(const char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"alpha") == 0) return (void *) alpha;
+  if (strcmp(str,"width") == 0) return (void *) width;
+  if (strcmp(str,"theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/EXTRA-MOLECULE/angle_gaussian.h

@@ -35,6 +35,7 @@ class AngleGaussian : 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:
   int *nterms;

src/EXTRA-MOLECULE/bond_gaussian.cpp

@@ -337,3 +337,14 @@ void BondGaussian::born_matrix(int type, double rsq, int /*i*/, int /*j*/, doubl
 
   du2 = - (force->boltz * bond_temperature[type]) * numerator / denominator;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *BondGaussian::extract(const char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"alpha") == 0) return (void *) alpha;
+  if (strcmp(str,"width") == 0) return (void *) width;
+  if (strcmp(str,"r0") == 0) return (void *) r0;
+  return nullptr;
+}

src/EXTRA-MOLECULE/bond_gaussian.h

@@ -36,6 +36,7 @@ class BondGaussian : 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:
   int *nterms;

src/EXTRA-MOLECULE/bond_nonlinear.cpp

@@ -228,6 +228,7 @@ void BondNonlinear::born_matrix(int type, double rsq, int /*i*/, int /*j*/, doub
 void *BondNonlinear::extract(const char *str, int &dim)
 {
   dim = 1;
+  if (strcmp(str,"lamda")==0) return (void*) lamda;
   if (strcmp(str,"epsilon")==0) return (void*) epsilon;
   if (strcmp(str,"r0")==0) return (void*) r0;
   return nullptr;

src/MESONT/angle_mesocnt.cpp

@@ -396,3 +396,14 @@ double AngleMesoCNT::single(int type, int i1, int i2, int i3)
   else
     return kb[type] * dtheta + thetab[type] * (kh[type] * thetab[type] - kb[type]);
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleMesoCNT::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/MESONT/angle_mesocnt.h

@@ -41,6 +41,7 @@ class AngleMesoCNT : 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:
   int *buckling;

src/MOFFF/angle_class2_p6.cpp

@@ -487,3 +487,19 @@ double AngleClass2P6::single(int type, int i1, int i2, int i3)
 
   return energy;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleClass2P6::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, "k5") == 0) return (void *) k5;
+  if (strcmp(str, "k6") == 0) return (void *) k6;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/MOFFF/angle_class2_p6.h

@@ -35,6 +35,7 @@ class AngleClass2P6 : 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, *k5, *k6;

src/MOFFF/angle_cosine_buck6d.cpp

@@ -383,3 +383,16 @@ double AngleCosineBuck6d::single(int type, int i1, int i2, int i3)
 
   return k[type]*(1.0+cos(tk));
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCosineBuck6d::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "k") == 0) return (void *) k;
+  if (strcmp(str, "multiplicity") == 0) return (void *) multiplicity;
+  if (strcmp(str, "th0") == 0) return (void *) th0;
+  return nullptr;
+}

src/MOFFF/angle_cosine_buck6d.h

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

src/MOLECULE/bond_morse.cpp

@@ -227,6 +227,8 @@ void BondMorse::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &
 void *BondMorse::extract(const char *str, int &dim)
 {
   dim = 1;
+  if (strcmp(str, "d0") == 0) return (void *) d0;
+  if (strcmp(str, "alpha") == 0) return (void *) alpha;
   if (strcmp(str, "r0") == 0) return (void *) r0;
   return nullptr;
 }

src/YAFF/angle_cross.cpp

@@ -342,3 +342,19 @@ double AngleCross::single(int type, int i1, int i2, int i3)
   double energy = kss[type]*dr1*dr2+kbs0[type]*dr1*dtheta + kbs1[type]*dr2*dtheta;
   return energy;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCross::extract(const char *str, int &dim)
+{
+  dim = 1;
+  if (strcmp(str, "r00") == 0) return (void *) r00;
+  if (strcmp(str, "r01") == 0) return (void *) r01;
+  if (strcmp(str, "kss") == 0) return (void *) kss;
+  if (strcmp(str, "kbs0") == 0) return (void *) kbs0;
+  if (strcmp(str, "kbs1") == 0) return (void *) kbs1;
+  if (strcmp(str, "theta0") == 0) return (void *) theta0;
+  return nullptr;
+}

src/YAFF/angle_cross.h

@@ -35,6 +35,7 @@ class AngleCross : 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 *kss, *kbs0, *kbs1, *r00, *r01, *theta0;

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

@@ -20,7 +20,13 @@ angle_coeff: ! |
   * ub 0.0 0.0
   4 ub -7.6 1.5537
 equilibrium: 4 1.9320794819577227 1.9687313962496038 2.1676989309769574 1.8936822384138474
-extract: ! ""
+extract: ! |
+  k2 1
+  k3 1
+  k4 1
+  k5 1
+  k6 1
+  theta0 1
 natoms: 29
 init_energy: 22.644137017730763
 init_stress: ! |2-

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

@@ -23,7 +23,13 @@ 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
+  k5 1
+  k6 1
+  theta0 1
 natoms: 29
 init_energy: 46.44080287964778
 init_stress: ! |2-

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

@@ -16,7 +16,13 @@ angle_coeff: ! |
   3 100.0 150.0 120.0 1.45 1.35  98.2
   4 250.0  90.0 110.0 1.05 1.20  99.9
 equilibrium: 4 1.8675022996339325 1.911135530933791 1.7139133254584316 1.7435839227423353
-extract: ! ""
+extract: ! |
+  r00 1
+  r01 1
+  kss 1
+  kbs0 1
+  kbs1 1
+  theta0 1
 natoms: 29
 init_energy: -138.93227715361755
 init_stress: ! |-

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

@@ -16,7 +16,10 @@ angle_coeff: ! |
   3  300.0 2 0.2189  8.66 88.1 0.5439  9.94 142.7
   4  300.0 2 0.0214 14.29 85.3 0.3934 18.22 118.1
 equilibrium: 4 1.5376350710070041 1.4887658519511628 1.5376350710070041 1.4887658519511628
-extract: ! ""
+extract: ! |
+  alpha 2
+  width 2
+  theta0 2
 natoms: 29
 init_energy: 57.794009158943744
 init_stress: ! |2-

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

@@ -18,6 +18,9 @@ bond_coeff: ! |
   5 0.97   532.50 -1282.90  2004.76
 equilibrium: 5 1.42 1.1 1.3 1.2 0.97
 extract: ! |
+  k2 1
+  k3 1
+  k4 1
   r0 1
 natoms: 29
 init_energy: 26.429859642132705

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

@@ -17,7 +17,10 @@ bond_coeff: ! |
   4 300.0 1 0.0100 0.098 2.45
   5 300.0 1 0.0100 0.098 2.85
 equilibrium: 5 1.45 1.37 1.61 2.45 2.85
-extract: ! ""
+extract: ! |
+  alpha 2
+  width 2
+  r0 2
 natoms: 29
 init_energy: 4638.6541482649545
 init_stress: ! |2-

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

@@ -18,6 +18,8 @@ bond_coeff: ! |
   5 7000.0 0.3 1.0
 equilibrium: 5 1.5 1.1 1.3 1.2 1
 extract: ! |
+  d0 1
+  alpha 1
   r0 1
 natoms: 29
 init_energy: 5.607154540709207

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

@@ -17,7 +17,10 @@ bond_coeff: ! |
   4 650.0 1.2 1.4
   5 450.0 1.0 1.1
 equilibrium: 5 1.5 1.1 1.3 1.2 1
-extract: ! ""
+extract: ! |
+  lamda 1
+  epsilon 1
+  r0 1
 natoms: 29
 init_energy: 1.9451728032820943
 init_stress: ! |-