git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14809 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/bond.cpp

@@ -16,6 +16,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "neighbor.h"
 #include "suffix.h"
 #include "atom_masks.h"
 #include "memory.h"
@@ -23,6 +24,8 @@
 
 using namespace LAMMPS_NS;
 
+enum{NONE,LINEAR,SPLINE};
+
 /* -----------------------------------------------------------------------
    set bond contribution to Vdwl energy to 0.0
    a particular bond style can override this
@@ -212,6 +215,74 @@ void Bond::ev_tally(int i, int j, int nlocal, int newton_bond,
   }
 }
 
+/* ----------------------------------------------------------------------
+   write a table of bond potential energy/force vs distance to a file
+------------------------------------------------------------------------- */
+
+void Bond::write_file(int narg, char **arg)
+{
+  if (narg != 6 && narg !=8) error->all(FLERR,"Illegal bond_write command");
+
+  // parse optional arguments
+
+  int itype = 0;
+  int jtype = 0;
+  if (narg == 8) {
+    itype = force->inumeric(FLERR,arg[6]);
+    jtype = force->inumeric(FLERR,arg[7]);
+    if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes)
+    error->all(FLERR,"Invalid atom types in bond_write command");
+  }
+
+  int btype = force->inumeric(FLERR,arg[0]);
+  int n = force->inumeric(FLERR,arg[1]);
+  double inner = force->numeric(FLERR,arg[2]);
+  double outer = force->numeric(FLERR,arg[3]);
+  if (inner <= 0.0 || inner >= outer)
+    error->all(FLERR,"Invalid rlo/rhi values in bond_write command");
+
+
+  double r0 = equilibrium_distance(btype);
+
+  // open file in append mode
+  // print header in format used by bond_style table
+
+  int me;
+  MPI_Comm_rank(world,&me);
+  FILE *fp;
+  if (me == 0) {
+    fp = fopen(arg[4],"a");
+    if (fp == NULL) error->one(FLERR,"Cannot open bond_write file");
+  }
+
+  // initialize potentials before evaluating bond potential
+  // insures all bond coeffs are set and force constants
+  // also initialize neighbor so that neighbor requests are processed
+  // NOTE: might be safest to just do lmp->init()
+
+  force->init();
+  neighbor->init();
+
+  if (me == 0) {
+    double r,e,f;
+
+    // evaluate energy and force at each of N distances
+    // note that Bond::single() takes r**2 and returns f/r.
+
+    fprintf(fp,"# Bond potential %s for bond type %d: i,r,energy,force\n",
+            force->bond_style,btype);
+    fprintf(fp,"\n%s\nN %d EQ %.15g\n\n",arg[5],n,r0);
+
+    const double dr = (outer-inner) / static_cast<double>(n-1);
+    for (int i = 0; i < n; i++) {
+      r = inner + dr * static_cast<double>(i);
+      e = single(btype,r*r,itype,jtype,f);
+      fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f*r);
+    }
+    fclose(fp);
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 double Bond::memory_usage()

src/bond.h

@@ -54,6 +54,8 @@ class Bond : protected Pointers {
   virtual unsigned int data_mask() {return datamask;}
   virtual unsigned int data_mask_ext() {return datamask_ext;}
 
+  void write_file(int, char**);
+
  protected:
   int suffix_flag;             // suffix compatibility flag
 

src/input.cpp

@@ -644,6 +644,7 @@ int Input::execute_command()
   else if (!strcmp(command,"atom_style")) atom_style();
   else if (!strcmp(command,"bond_coeff")) bond_coeff();
   else if (!strcmp(command,"bond_style")) bond_style();
+  else if (!strcmp(command,"bond_write")) bond_write();
   else if (!strcmp(command,"boundary")) boundary();
   else if (!strcmp(command,"box")) box();
   else if (!strcmp(command,"comm_modify")) comm_modify();
@@ -1281,6 +1282,17 @@ void Input::bond_style()
 
 /* ---------------------------------------------------------------------- */
 
+void Input::bond_write()
+{
+  if (atom->avec->bonds_allow == 0)
+    error->all(FLERR,"Bond_write command when no bonds allowed");
+  if (force->bond == NULL)
+    error->all(FLERR,"Bond_write command before bond_style is defined");
+  else force->bond->write_file(narg,arg);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Input::boundary()
 {
   if (domain->box_exist)

src/input.h

@@ -86,6 +86,7 @@ class Input : protected Pointers {
   void atom_style();
   void bond_coeff();
   void bond_style();
+  void bond_write();
   void boundary();
   void box();
   void comm_modify();

src/pair.cpp

@@ -1574,9 +1574,9 @@ void Pair::write_file(int narg, char **arg)
     fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n",
             force->pair_style,itype,jtype);
     if (style == RLINEAR)
-      fprintf(fp,"\n%s\nN %d R %g %g\n\n",arg[7],n,inner,outer);
+      fprintf(fp,"\n%s\nN %d R %.15g %.15g\n\n",arg[7],n,inner,outer);
     if (style == RSQ)
-      fprintf(fp,"\n%s\nN %d RSQ %g %g\n\n",arg[7],n,inner,outer);
+      fprintf(fp,"\n%s\nN %d RSQ %.15g %.15g\n\n",arg[7],n,inner,outer);
   }
 
   // initialize potentials before evaluating pair potential
@@ -1618,7 +1618,7 @@ void Pair::write_file(int narg, char **arg)
     init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits);
     int ntable = 1 << n;
     if (me == 0)
-      fprintf(fp,"\n%s\nN %d BITMAP %g %g\n\n",arg[7],ntable,inner,outer);
+      fprintf(fp,"\n%s\nN %d BITMAP %.15g %.15g\n\n",arg[7],ntable,inner,outer);
     n = ntable;
   }
 
@@ -1647,7 +1647,7 @@ void Pair::write_file(int narg, char **arg)
       e = single(0,1,itype,jtype,rsq,1.0,1.0,f);
       f *= r;
     } else e = f = 0.0;
-    if (me == 0) fprintf(fp,"%d %g %g %g\n",i+1,r,e,f);
+    if (me == 0) fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f);
   }
 
   // restore original vecs that were swapped in for

src/pair_table.cpp

@@ -395,13 +395,16 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
   union_int_float_t rsq_lookup;
 
   int rerror = 0;
+  int cerror = 0;
 
   fgets(line,MAXLINE,fp);
   for (int i = 0; i < tb->ninput; i++) {
-    fgets(line,MAXLINE,fp);
-    sscanf(line,"%d %lg %lg %lg",&itmp,&rfile,&tb->efile[i],&tb->ffile[i]);
-    rnew = rfile;
+    if (NULL == fgets(line,MAXLINE,fp))
+      error->one(FLERR,"Premature end of file in pair table");
+    if (4 != sscanf(line,"%d %lg %lg %lg",
+                    &itmp,&rfile,&tb->efile[i],&tb->ffile[i]))  ++cerror;
 
+    rnew = rfile;
     if (tb->rflag == RLINEAR)
       rnew = tb->rlo + (tb->rhi - tb->rlo)*i/(tb->ninput-1);
     else if (tb->rflag == RSQ) {
@@ -419,6 +422,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
     }
 
     if (tb->rflag && fabs(rnew-rfile)/rfile > EPSILONR) rerror++;
+
     tb->rfile[i] = rnew;
   }
 
@@ -452,8 +456,8 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
 
   if (ferror) {
     char str[128];
-    sprintf(str,"%d force values in table are inconsistent with -dE/dr; "
-            "should only be mistakenly flagged at inflection points",ferror);
+    sprintf(str,"%d of %d force values in table are inconsistent with -dE/dr.\n"
+            "  Should only be flagged at inflection points",ferror,tb->ninput);
     error->warning(FLERR,str);
   }
 
@@ -461,8 +465,17 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
 
   if (rerror) {
     char str[128];
-    sprintf(str,"%d distance values in table differ signifcantly "
-            "from re-computed values",rerror);
+    sprintf(str,"%d of %d distance values in table with relative error\n"
+            "  over %g to re-computed values",rerror,tb->ninput,EPSILONR);
+    error->warning(FLERR,str);
+  }
+
+  // warn if data was read incompletely, e.g. columns were missing
+
+  if (cerror) {
+    char str[128];
+    sprintf(str,"%d of %d lines in table were incomplete\n"
+            "  or could not be parsed completely",cerror,tb->ninput);
     error->warning(FLERR,str);
   }
 }