UB testing

src/AMOEBA/angle_amoeba.cpp

@@ -103,7 +103,9 @@ void AngleAmoeba::compute(int eflag, int vflag)
 
   ev_init(eflag,vflag);
 
+  // DEBUG
   int ubn = 0;
+  double eub = 0.0;
 
   for (n = 0; n < nanglelist; n++) {
     i1 = anglelist[n][0];
@@ -132,11 +134,12 @@ void AngleAmoeba::compute(int eflag, int vflag)
 
     if (uflag) {
       ubn++;
-      tinker_urey_bradley(i1,i3,type,eflag);
+      eub += tinker_urey_bradley(i1,i3,type,eflag);
     }
   }
 
-  printf("UB N %d\n",ubn);
+  // DEBUG
+  printf("UreyBradley n %d eng %g\n",ubn,eub);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -540,7 +543,7 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag)
 
 /* ---------------------------------------------------------------------- */
 
-void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
+double AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
 {
   double delx,dely,delz;
   double rsq,r,dr,rk;
@@ -582,6 +585,8 @@ void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
   }
 
   if (evflag) ev_tally2(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
+  
+  return ebond;
 }
 
 /* ---------------------------------------------------------------------- */

src/AMOEBA/angle_amoeba.h

@@ -46,7 +46,7 @@ class AngleAmoeba : public Angle {
   void tinker_angle(int, int, int, int, int);
   void tinker_anglep(int, int, int, int, int);
   void tinker_bondangle(int, int, int, int, int);
-  void tinker_urey_bradley(int, int, int, int);
+  double tinker_urey_bradley(int, int, int, int);
   void allocate();
 };
 

src/read_data.cpp

@@ -658,6 +658,13 @@ void ReadData::command(int narg, char **arg)
           error->all(FLERR,"Must define angle_style before BondAngle Coeffs");
         if (firstpass) anglecoeffs(2);
         else skip_lines(nangletypes);
+      } else if (strcmp(keyword,"UreyBradley Coeffs") == 0) {
+        if (atom->avec->angles_allow == 0)
+          error->all(FLERR,"Invalid data file section: UreyBradley Coeffs");
+        if (force->angle == nullptr)
+          error->all(FLERR,"Must define angle_style before UreyBradley Coeffs");
+        if (firstpass) anglecoeffs(3);
+        else skip_lines(nangletypes);
 
       } else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
         if (atom->avec->dihedrals_allow == 0)
@@ -1852,6 +1859,7 @@ void ReadData::anglecoeffs(int which)
     if (which == 0) parse_coeffs(buf,nullptr,0,1,aoffset);
     else if (which == 1) parse_coeffs(buf,"bb",0,1,aoffset);
     else if (which == 2) parse_coeffs(buf,"ba",0,1,aoffset);
+    else if (which == 3) parse_coeffs(buf,"ub",0,1,aoffset);
     if (ncoeffarg == 0) error->all(FLERR,"Unexpected empty line in AngleCoeffs section");
     force->angle->coeff(ncoeffarg,coeffarg);
     buf = next + 1;

tools/tinker/data.py

@@ -250,6 +250,11 @@ class data:
   def write(self,file):
     f = open(file,"w")
     print >>f,self.title
+ 
+    # write any keywords in standard list hkeywords
+    #   in the order they are in hkeywords
+    # then write any extra keywords at end of header section
+   
     for keyword in hkeywords:
       if self.headers.has_key(keyword):
         if keyword == "xlo xhi" or keyword == "ylo yhi" or \
@@ -261,12 +266,30 @@ class data:
 	  print >>f,triple[0],triple[1],triple[2],keyword
         else:
 	  print >>f,self.headers[keyword],keyword
+
+    for keyword in self.headers.keys():
+      if keyword not in hkeywords:
+        print >>f,self.headers[keyword],keyword
+
+    # write any sections in standard list skeywords
+    #   in the order they are in skeywords
+    # then write any extra sections at end of file
+
     for pair in skeywords:
       keyword = pair[0]
       if self.sections.has_key(keyword):
         print >>f,"\n%s\n" % keyword
         for line in self.sections[keyword]:
 	  print >>f,line,
+
+    skeyfirst = [pair[0] for pair in skeywords]
+    
+    for keyword in self.sections.keys():
+      if keyword not in skeyfirst:
+        print >>f,"\n%s\n" % keyword
+        for line in self.sections[keyword]:
+	  print >>f,line,
+
     f.close()
 
   # --------------------------------------------------------------------
@@ -348,35 +371,32 @@ class data:
     return self.headers["atom types"]
 
 # --------------------------------------------------------------------
-# data file keywords, both header and main sections
+# standard data file keywords, both header and main sections
 
 hkeywords = ["atoms","ellipsoids","lines","triangles","bodies",
-             "bonds","angles","dihedrals","impropers","pitorsions",
+             "bonds","angles","dihedrals","impropers",
 	     "atom types","bond types","angle types","dihedral types",
-	     "improper types","pitorsion types",
+	     "improper types",
              "xlo xhi","ylo yhi","zlo zhi","xy xz yz"]
 
 skeywords = [["Masses","atom types"],
              ["Atoms","atoms"],["Ellipsoids","ellipsoids"],
              ["Lines","lines"],["Triangles","triangles"],["Bodies","bodies"],
+             ["Velocities","atoms"],
              ["Bonds","bonds"],
 	     ["Angles","angles"],
              ["Dihedrals","dihedrals"],
 	     ["Impropers","impropers"],
-	     ["PiTorsions","pitorsions"],
-             ["Velocities","atoms"],
              ["Pair Coeffs","atom types"],
-	     ["Bond Coeffs","bond types"],["Angle Coeffs","angle types"],
+	     ["Bond Coeffs","bond types"],
+             ["Angle Coeffs","angle types"],
 	     ["Dihedral Coeffs","dihedral types"],
 	     ["Improper Coeffs","improper types"],
-	     ["PiTorsion Coeffs","pitorsion types"],
              ["BondBond Coeffs","angle types"],
              ["BondAngle Coeffs","angle types"],
-             ["UreyBradley Coeffs","angle types"],
              ["MiddleBondTorsion Coeffs","dihedral types"],
              ["EndBondTorsion Coeffs","dihedral types"],
              ["AngleTorsion Coeffs","dihedral types"],
              ["AngleAngleTorsion Coeffs","dihedral types"],
              ["BondBond13 Coeffs","dihedral types"],
-             ["AngleAngle Coeffs","improper types"],
+             ["AngleAngle Coeffs","improper types"]]
-             ["Molecules","atoms"]]