add stretch-bend cross term

examples/amoeba/in.ubiquitin

@@ -8,6 +8,7 @@ atom_style amoeba
 bond_style      class2
 angle_style     amoeba
 dihedral_style  fourier
+improper_style  amoeba
 
 # per-atom properties required by AMOEBA or HIPPO
 

src/AMOEBA/angle_amoeba.cpp

@@ -169,8 +169,9 @@ void AngleAmoeba::compute(int eflag, int vflag)
     f3[1] = a22*dely2 + a12*dely1;
     f3[2] = a22*delz2 + a12*delz1;
 
-    if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
+    //if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
-                 k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
+    //             k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
+    eangle = 0.0;
 
     // force & energy for bond-angle term
     // bond-stretch cross term in Tinker
@@ -216,6 +217,9 @@ void AngleAmoeba::compute(int eflag, int vflag)
 
     if (eflag) eangle += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
 
+    printf("Stretch-Bend: %ld %d %d: eng %g\n",
+           atom->tag[i1],atom->tag[i2],atom->tag[i3],eangle);
+
     // apply force to each of 3 atoms
 
     if (newton_bond || i1 < nlocal) {
@@ -445,6 +449,9 @@ void AngleAmoeba::allocate()
   memory->create(ba_r2,n+1,"angle:ba_r2");
 
   memory->create(setflag,n+1,"angle:setflag");
+  memory->create(setflag_a,n+1,"angle:setflag");
+  memory->create(setflag_ba,n+1,"angle:setflag");
+
   for (int i = 1; i <= n; i++) setflag[i] = setflag_a[i] = setflag_ba[i] = 0;
 }
 
@@ -500,6 +507,7 @@ void AngleAmoeba::coeff(int narg, char **arg)
       k4[i] = k4_one;
       k5[i] = k5_one;
       k6[i] = k6_one;
+      setflag_a[i] = 1;
       count++;
     }
   }

src/AMOEBA/improper_amoeba.cpp

@@ -273,15 +273,11 @@ void ImproperAmoeba::init_style()
   if (!pair) pair = force->pair_match("hippo",1,0);
   if (!pair) error->all(FLERR,"Improper amoeba could not find pair amoega");
   
-  
   int dim;
   opbend_cubic = *(double *) pair->extract("opbend_cubic",dim);
   opbend_quartic = *(double *) pair->extract("opbend_quartic",dim);
   opbend_pentic = *(double *) pair->extract("opbend_pentic",dim);
   opbend_sextic = *(double *) pair->extract("opbend_sextic",dim);
-
-  printf("OPBEND %g %g %g %g\n",
-         opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic);
 }
 
 /* ----------------------------------------------------------------------
@@ -301,9 +297,8 @@ void ImproperAmoeba::read_restart(FILE *fp)
 {
   allocate();
 
-  if (comm->me == 0) {
+  if (comm->me == 0)
     utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
-  }
   MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
 
   for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;

tools/tinker2lmp.py

@@ -36,7 +36,7 @@ def error(txt=""):
     print "  -nopbc"
     print "  -pbc xhi yhi zhi"
   else: print "ERROR:",txt
-  sys.exit()
+  #sys.exit()
 
 # read and store values from a Tinker xyz file
 
@@ -251,13 +251,13 @@ class PRMfile:
   
   def bondangles(self):
     params = []
-    iline = self.find_section("Stretch Bend Parameters")
+    iline = self.find_section("Stretch-Bend Parameters")
     if iline < 0: return params
     iline += 3
 
     bdict = {}
-    for m,params in enumerate(self.bondparams):
+    for m,bparams in enumerate(self.bondparams):
-      bdict[(params[0],params[1])] = params[2]
+      bdict[(bparams[0],bparams[1])] = bparams[2]
     
     while iline < self.nlines:
       words = self.lines[iline].split()
@@ -271,14 +271,27 @@ class PRMfile:
           value2 = float(words[5])
           lmp1 = value1
           lmp2 = value2
-          lmp3 = lmp4 = 0.0
+          
-          if (class2,class1) in bdict: lmp3 = bdict[(class2,class1)]
+          if (class1,class2) in bdict:
-          if (class1,class2) in bdict: lmp3 = bdict[(class1,class2)]
+            lmp3 = bdict[(class1,class2)]
-          if (class2,class3) in bdict: lmp4 = bdict[(class2,class3)]
+          elif (class2,class1) in bdict:
-          if (class3,class2) in bdict: lmp4 = bdict[(class3,class2)]
+            lmp3 = bdict[(class2,class1)]
-          if lmp3 == 0.0 or lmp4 == 0.0:
+          else:
-            print "Bond in BondAngle term not found",class1,class2,class3
+            error("1st bond in BondAngle term not found: %d %d %d" % \
-            sys.exit()
+                  (class1,class2,class3))
+            # NOTE: just for debugging
+            lmp3 = 0.0
+
+          if (class2,class3) in bdict:
+            lmp4 = bdict[(class2,class3)]
+          elif (class3,class2) in bdict:
+            lmp4 = bdict[(class3,class2)]
+          else:
+            error("2nd bond in BondAngle term not found: %d %d %d" % \
+                  (class1,class2,class3))
+            # NOTE: just for debugging
+            lmp4 = 0.0
+            
           params.append((class1,class2,class3,lmp1,lmp2,lmp3,lmp4))
       iline += 1
     return params
@@ -301,12 +314,11 @@ class PRMfile:
           class4 = int(words[4])
           
           if len(words) <= 5:
-            print "Torsion has no params",class1,class2,class3,class4
+            error("Torsion has no params: %d %d %d %d" % \
-            sys.exit()
+                  (class1,class2,class3,class4))
           if (len(words)-5) % 3: 
-            print "Torsion does not have triplets of params", \
+            error("Torsion does not have triplets of params: %d %d %d %d" % \
-              class1,class2,class3,class4
+                  (class1,class2,class3,class4))
-            sys.exit()
 
           mfourier = (len(words)-5) / 3
           oneparams = [class1,class2,class3,class4,mfourier]
@@ -603,9 +615,7 @@ for atom1,atom2 in blist:
   
   if (c1,c2) in bdict: m,params = bdict[(c1,c2)]
   elif (c2,c1) in bdict: m,params = bdict[(c2,c1)]
-  else:
+  else: error("Bond not found: %d %d: %d %d" % (atom1,atom2,c1,c2))
-    print "Bond not found",atom1,atom2,c1,c2
-    sys.exit()
     
   if not flags[m]:
     v1,v2,v3,v4 = params[2:]
@@ -615,7 +625,6 @@ for atom1,atom2 in blist:
 
 # generate atype = LAMMPS type of each angle
 # generate aparams = LAMMPS params for each angle type
-# generate baparams = LAMMPS bond-angle params for each angle type
 # flags[i] = which LAMMPS angle type (1-N) the Tinker FF file angle I is
 #        0 = none
 # Tinker FF file angle entries can have 1, 2, or 3 options
@@ -636,12 +645,11 @@ for m,params in enumerate(prm.angleparams):
   noptions += n
   
 flags = noptions*[0]
-#baflags = len(baprm)*[0]
 atype = []
 aparams = []
-baparams = []
 
-for atom1,atom2,atom3 in alist:
+for i,one in enumerate(alist):
+  atom1,atom2,atom3 = one
   type1 = type[atom1-1]
   type2 = type[atom2-1]
   type3 = type[atom3-1]
@@ -651,7 +659,18 @@ for atom1,atom2,atom3 in alist:
 
   if (c1,c2,c3) in adict or (c3,c2,c1) in adict:
     if (c1,c2,c3) in adict: m,params = adict[(c1,c2,c3)]
-    if (c3,c2,c1) in adict: m,params = adict[(c3,c2,c1)]
+
+    # IMPORTANT subtlety
+    # flip order of 3 atoms in alist if the angle
+    #   matches Angle Bending section of PRM file in reverse order
+    # necessary b/c BondAngle coeffs will be generated with r1,r2 params
+    #   from Bond Stretching section of PRM file
+    # since in general r1 != r2, the LAMMPS AngleAmoeba class requires
+    #   the 3 atoms in the angle be in the order that matches r1 and r2
+    
+    if (c3,c2,c1) in adict:
+      m,params = adict[(c3,c2,c1)]
+      alist[i] = (atom3,atom2,atom1)
       
     # params is a sequence of 1 or 2 or 3 options
     # which = which of 1,2,3 options this atom triplet matches
@@ -679,7 +698,7 @@ for atom1,atom2,atom3 in alist:
         print "  angle atom classes:",c1,c2,c3
         print "  Tinker FF file param options:",len(params[3])
         print "  Nbonds and hydrogen count:",nbonds,hcount
-        #sys.exit()      // NOTE: allow this for now
+        #sys.exit()      NOTE: allow this for now
 
       if hcount == 0: which = 1
       elif hcount == 1:
@@ -702,7 +721,7 @@ for atom1,atom2,atom3 in alist:
         print "  angle atom classes:",c1,c2,c3
         print "  Tinker FF file param options:",len(params[3])
         print "  Nbonds and hydrogen count:",nbonds,hcount
-        #sys.exit()     // NOTE: allow this for now
+        #sys.exit()     NOTE: allow this for now
         
       if hcount == 0: which = 1
       elif hcount == 1:
@@ -713,8 +732,7 @@ for atom1,atom2,atom3 in alist:
         m += 2
 
   else:
-    print "Angle not found",atom1,atom2,atom3,c1,c2,c3
+    error("Angle not found: %d %d %d: %d %d %d" % (atom1,atom2,atom3,c1,c2,c3))
-    sys.exit()
     
   if not flags[m]:
     pflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
@@ -722,24 +740,36 @@ for atom1,atom2,atom3 in alist:
     flags[m] = len(aparams)
   atype.append(flags[m])
 
-  # NOTE: baparams may need to be flipped if match is 3,2,1 instead of 1,2,3
+# augment the aparams with bond-angle cross terms from bondangleparams
+# generate baparams = LAMMPS bond-angle params for each angle type
+# sbdict = dictionary for angle tuples in bongangleparams
 
-  # NOTE: mismatch between angle and bondangle params may not be handled right
+sbdict = {}
-  # should be a new LAMMPS type if either angle or bondangle params do not match?
+for v1,v2,v3,v4,v5,v6,v7 in prm.bondangleparams:
-  #for m,params in enumerate(baprm):
+  if (v1,v2,v3) in sbdict: continue
-  #  c1,c2,c3,v1,v2,v3,v4 = params
+  sbdict[(v1,v2,v3)] = (v4,v5,v6,v7)
-  #  if (c1 == class1 and c2 == class2 and c3 == class3) or \
+
-  #     (c1 == class3 and c2 == class2 and c3 == class1):
+baparams = []
-  #    found += 1
+
-  #    if baflags[m]:
+for itype in range(len(aparams)):
-  #      continue
+  iangle = atype.index(itype+1) 
-  #      #atype.append(baflags[m])
+  atom1,atom2,atom3 = alist[iangle]
-  #    else:
+  type1 = type[atom1-1]
-  #      baparams.append((v1,v2,v3,v4))
+  type2 = type[atom2-1]
-  #      baflags[m] = len(baparams)
+  type3 = type[atom3-1]
-  #      #atype.append(baflags[m])
+  c1 = classes[type1-1]
-  #    break
+  c2 = classes[type2-1]
-  #  if found != 1: print "Not found",atom1,atom2,atom3,class1,class2,class3
+  c3 = classes[type3-1]
+
+  if (c1,c2,c3) in sbdict:
+    n1,n2,r1,r2 = sbdict[(c1,c2,c3)]
+  elif (c3,c2,c1) in sbdict:
+    n1,n2,r1,r2 = sbdict[(c3,c2,c1)]
+  else:
+    print "Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3)
+    n1,n2,r1,r2 = 4*[0.0]
+    
+  baparams.append((n1,n2,r1,r2))
 
 # generate dtype = LAMMPS type of each dihedral
 # generate dparams = LAMMPS params for each dihedral type
@@ -774,8 +804,8 @@ for atom1,atom2,atom3,atom4 in dlist:
   if (c1,c2,c3,c4) in ddict: m,params = ddict[(c1,c2,c3,c4)]
   elif (c4,c3,c2,c1) in ddict: m,params = ddict[(c4,c3,c2,c1)]
   else:
-    print "Dihedral not found",atom1,atom2,atom3,atom4,c1,c2,c3,c4
+    error("Dihedral not found: %d %d %d %d: %d %d %d %d" % \
-    sys.exit()
+          (atom1,atom2,atom3,atom4,c1,c2,c3,c4))
 
   if not flags[m]:
     oneparams = params[4:]
@@ -948,17 +978,12 @@ if nangles:
     lines.append(line+'\n')
   d.sections["Angle Coeffs"] = lines
 
-  #lines = []
+  lines = []
-  #for i,one in enumerate(aparams):
+  for i,one in enumerate(baparams):
-  #  line = "%d %g %g %g" % (i+1,0.0,0.0,0.0)
+    strone = [str(single) for single in one]
-  #  lines.append(line+'\n')
+    line = "%d %s" % (i+1,' '.join(strone))
-  #d.sections["BondBond Coeffs"] = lines
+    lines.append(line+'\n')
-
+  d.sections["BondAngle Coeffs"] = lines
-  #lines = []
-  #for i,one in enumerate(aparams):
-  #  line = "%d %g %g %g %g" % (i+1,0.0,0.0,0.0,0.0)
-  #  lines.append(line+'\n')
-  # d.sections["BondAngle Coeffs"] = lines
 
   lines = [] 
   for i,one in enumerate(alist):