Made fix adapt work with bond_harmonic.

doc/src/fix_adapt.txt

@@ -22,6 +22,11 @@ attribute = {pair} or {kspace} or {atom} :l
     pparam = parameter to adapt over time
     I,J = type pair(s) to set parameter for
     v_name = variable with name that calculates value of pparam
+  {bond} args = bstyle bparam I v_name
+    bstyle = bond style name, e.g. harmonic
+    bparam = parameter to adapt over time
+    I = type bond to set parameter for
+    v_name = variable with name that calculates value of bparam
   {kspace} arg = v_name
     v_name = variable with name that calculates scale factor on K-space terms
   {atom} args = aparam v_name
@@ -42,7 +47,10 @@ keyword = {scale} or {reset} :l
 fix 1 all adapt 1 pair soft a 1 1 v_prefactor
 fix 1 all adapt 1 pair soft a 2* 3 v_prefactor
 fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
-fix 1 all adapt 10 atom diameter v_size :pre
+fix 1 all adapt 10 atom diameter v_size
+
+variable ramp_up equal "ramp(0.01,0.5)"
+fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up :pre
 
 [Description:]
 
@@ -192,6 +200,19 @@ fix 1 all adapt 1 pair soft a * * v_prefactor :pre
 
 :line
 
+The {bond} keyword uses the specified variable to change the value of
+a bond coefficient over time, very similar to how the {pair} keyword
+operates. The only difference is that now a bond coefficient for a
+given bond type is adapted.
+
+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
+
+"harmonic"_bond_harmonic.html: k,r0: type bonds :tb(c=3,s=:)
+
+:line
+
 The {kspace} keyword used the specified variable as a scale factor on
 the energy, forces, virial calculated by whatever K-Space solver is
 defined by the "kspace_style"_kspace_style.html command.  If the

src/bond.cpp

@@ -292,3 +292,14 @@ double Bond::memory_usage()
   bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }
+
+/* -----------------------------------------------------------------------
+   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");
+
+  init();
+}

src/bond.h

@@ -30,6 +30,8 @@ class Bond : protected Pointers {
   double virial[6];               // accumulated virial
   double *eatom,**vatom;          // accumulated per-atom energy/virial
 
+  int reinitflag;                // 1 if compatible with fix adapt and alike
+
   // KOKKOS host/device flag and data masks
 
   ExecutionSpace execution_space;
@@ -49,6 +51,8 @@ class Bond : protected Pointers {
   virtual void write_data(FILE *) {}
   virtual double single(int, double, int, int, double &) = 0;
   virtual double memory_usage();
+  virtual void *extract(char *, int &) {return NULL;}
+  virtual void reinit();
 
   void write_file(int, char**);
 

src/fix_adapt.cpp

@@ -16,6 +16,7 @@
 #include <stdlib.h>
 #include "fix_adapt.h"
 #include "atom.h"
+#include "bond.h"
 #include "update.h"
 #include "group.h"
 #include "modify.h"
@@ -35,7 +36,7 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{PAIR,KSPACE,ATOM};
+enum{PAIR,KSPACE,ATOM,BOND};
 enum{DIAMETER,CHARGE};
 
 /* ---------------------------------------------------------------------- */
@@ -68,6 +69,10 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
       if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command");
       nadapt++;
       iarg += 3;
+    } else if (strcmp(arg[iarg],"bond") == 0 ){
+      if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command");
+      nadapt++;
+      iarg += 5;
     } else break;
   }
 
@@ -103,6 +108,25 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
       } else error->all(FLERR,"Illegal fix adapt command");
       nadapt++;
       iarg += 6;
+    } else if (strcmp(arg[iarg],"bond") == 0 ){
+      if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command");
+      adapt[nadapt].which = BOND;
+      int n = strlen(arg[iarg+1]) + 1;
+      adapt[nadapt].bstyle = new char[n];
+      strcpy(adapt[nadapt].bstyle,arg[iarg+1]);
+      n = strlen(arg[iarg+2]) + 1;
+      adapt[nadapt].bparam = new char[n];
+      adapt[nadapt].bond = NULL;
+      strcpy(adapt[nadapt].bparam,arg[iarg+2]);
+      force->bounds(FLERR,arg[iarg+3],atom->ntypes,
+                    adapt[nadapt].ilo,adapt[nadapt].ihi);
+      if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) {
+        n = strlen(&arg[iarg+4][2]) + 1;
+        adapt[nadapt].var = new char[n];
+        strcpy(adapt[nadapt].var,&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;
@@ -160,6 +184,13 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
   for (int m = 0; m < nadapt; m++)
     if (adapt[m].which == PAIR)
       memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
+
+  // allocate bond style arrays:
+  n = atom->nbondtypes;
+  for (int m = 0; m < nadapt; ++m)
+    if (adapt[m].which == BOND)
+      // For now just use same storage and fake it to be one-dimensional:
+      memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -172,6 +203,10 @@ FixAdapt::~FixAdapt()
       delete [] adapt[m].pstyle;
       delete [] adapt[m].pparam;
       memory->destroy(adapt[m].array_orig);
+    } else if (adapt[m].which == BOND) {
+      delete [] adapt[m].bstyle;
+      delete [] adapt[m].bparam;
+      memory->destroy(adapt[m].vector_orig);
     }
   }
   delete [] adapt;
@@ -282,6 +317,7 @@ void FixAdapt::init()
   // setup and error checks
 
   anypair = 0;
+  anybond = 0;
 
   for (int m = 0; m < nadapt; m++) {
     Adapt *ad = &adapt[m];
@@ -350,7 +386,48 @@ void FixAdapt::init()
       }
 
       delete [] pstyle;
+    } else if (ad->which == BOND){
+      ad->bond = NULL;
+      anybond = 1;
+      // Use same routines from pair to strip any suffices:
+      
+      int n = strlen(ad->bstyle) + 1;
+      char *bstyle = new char[n];
+      strcpy(bstyle,ad->bstyle);
 
+      if (lmp->suffix_enable) {
+        int len = 2 + strlen(bstyle) + strlen(lmp->suffix);
+        char *bsuffix = new char[len];
+        strcpy(bsuffix,bstyle);
+        strcat(bsuffix,"/");
+        strcat(bsuffix,lmp->suffix);
+        ad->bond = force->bond_match(bsuffix);
+        delete [] bsuffix;
+      }
+      // If not set grab regular one instead:
+      if (ad->bond == NULL) ad->bond = force->bond_match(bstyle);
+      if (ad->bond == NULL )
+        error->all(FLERR,"Fix adapt bond style does not exist");
+
+      void *ptr = ad->bond->extract(ad->bparam,ad->bdim);
+      if( comm->me == 0 ){
+        fprintf(screen,"Attempting to extract %s\n", ad->bparam);
+        fprintf(screen,"Got back %p\n", ptr);
+      }
+      
+      if (ptr == NULL)
+        error->all(FLERR,"Fix adapt bond style param not supported");
+
+      // For bond styles you should use a vector
+
+      if (ad->bdim == 1) ad->vector = (double *) ptr;
+
+      if (strcmp(force->bond_style,"hybrid") == 0 ||
+          strcmp(force->bond_style,"hybrid_overlay") == 0)
+        error->all(FLERR,"Fix adapt does not support bond_style hybrid");
+
+      delete [] bstyle;
+        
     } else if (ad->which == KSPACE) {
       if (force->kspace == NULL)
         error->all(FLERR,"Fix adapt kspace style does not exist");
@@ -368,7 +445,7 @@ void FixAdapt::init()
     }
   }
 
-  // make copy of original pair array values
+  // make copy of original pair/bond array values
 
   for (int m = 0; m < nadapt; m++) {
     Adapt *ad = &adapt[m];
@@ -378,7 +455,12 @@ void FixAdapt::init()
           ad->array_orig[i][j] = ad->array[i][j];
     }else if (ad->which == PAIR && ad->pdim == 0){
       ad->scalar_orig = *ad->scalar;
+      
+    }else if (ad->which == BOND && ad->bdim == 1){
+      for (i = ad->ilo; i <= ad->ihi; ++i )
+        ad->vector_orig[i] = ad->vector[i];
     }
+    
   }
 
   // fixes that store initial per-atom values
@@ -470,6 +552,18 @@ void FixAdapt::change_settings()
               ad->array[i][j] = value;
       }
 
+    // set bond type array values:
+      
+    } else if (ad->which == BOND) {
+      if (ad->pdim == 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) {
@@ -532,6 +626,14 @@ void FixAdapt::change_settings()
       }
     }
   }
+  if (anybond) {
+    for (int m = 0; m < nadapt; ++m ) {
+      Adapt *ad = &adapt[m];
+      if (ad->which == BOND) {
+        ad->bond->reinit();
+      }
+    }
+  }
 
   // reset KSpace charges if charges have changed
 
@@ -554,6 +656,12 @@ void FixAdapt::restore_settings()
             ad->array[i][j] = ad->array_orig[i][j];
       }
 
+    } else if (ad->which == BOND) {
+      if (ad->pdim == 1) {
+        for (int i = ad->ilo; i <= ad->ihi; i++)
+          ad->vector[i] = ad->vector_orig[i];
+      }
+      
     } else if (ad->which == KSPACE) {
       *kspace_scale = 1.0;
 
@@ -588,6 +696,7 @@ void FixAdapt::restore_settings()
   }
 
   if (anypair) force->pair->reinit();
+  if (anybond) force->bond->reinit();
   if (chgflag && force->kspace) force->kspace->qsum_qsq();
 }
 

src/fix_adapt.h

@@ -43,7 +43,7 @@ class FixAdapt : public Fix {
 
  private:
   int nadapt,resetflag,scaleflag;
-  int anypair;
+  int anypair, anybond;
   int nlevels_respa;
   char *id_fix_diam,*id_fix_chg;
   class FixStore *fix_diam,*fix_chg;
@@ -52,12 +52,15 @@ class FixAdapt : public Fix {
     int which,ivar;
     char *var;
     char *pstyle,*pparam;
+    char *bstyle,*bparam;
     int ilo,ihi,jlo,jhi;
-    int pdim;
+    int pdim,bdim;
     double *scalar,scalar_orig;
+    double *vector,*vector_orig;
     double **array,**array_orig;
     int aparam;
     class Pair *pair;
+    class Bond *bond;
   };
 
   Adapt *adapt;