From a2f62ae2dbb609fc7d768421113cb697c1d75ac9 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Tue, 19 Oct 2021 16:10:51 -0600
Subject: [PATCH] angles issue with angle vs anglep

---
 src/AMOEBA/angle_amoeba.cpp | 19 ++++++-------------
 tools/tinker2lmp.py         | 18 ------------------
 2 files changed, 6 insertions(+), 31 deletions(-)

diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp
index 4476615a50..6faf2edbdd 100644
--- a/src/AMOEBA/angle_amoeba.cpp
+++ b/src/AMOEBA/angle_amoeba.cpp
@@ -79,6 +79,8 @@ void AngleAmoeba::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   int **anglelist = neighbor->anglelist;
+  int **nspecial = atom->nspecial;
+
   int nanglelist = neighbor->nanglelist;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
@@ -90,11 +92,11 @@ void AngleAmoeba::compute(int eflag, int vflag)
     type = anglelist[n][3];
 
     // tflag = 0 for "angle", 1 for "anglep" in Tinker PRM file
-    // atom 2 must have 3 bond partners to invoke "anglep" option
+    // atom 2 must have 3 bond partners to invoke anglep() variant
 
     tflag = pflag[type];
 
-    if (tflag && atom->num_bond[i2] == 3) {
+    if (tflag && nspecial[i2][0] == 3) {
       anglep(i1,i2,i3,type,eflag);
       continue;
     }
@@ -203,7 +205,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
 
   double **x = atom->x;
   double **f = atom->f;
-  tagint **bond_atom = atom->bond_atom;
+  tagint **special = atom->special;
   int nlocal = atom->nlocal;
   int newton_bond = force->newton_bond;
 
@@ -213,7 +215,7 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
   i3tag = atom->tag[i3];
 
   for (int ibond = 0; ibond < 3; ibond++) {
-    i4tag = bond_atom[i2][ibond];
+    i4tag = special[i2][ibond];
     if (i4tag != i1tag && i4tag != i3tag) break;
   }
 
@@ -290,15 +292,6 @@ void AngleAmoeba::anglep(int i1, int i2, int i3, int type, int eflag)
   if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
                k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
 
-  /*
-  printf("ANGLEP: %d %d %d %d: eng %g\n",
-         atom->tag[i1],
-         atom->tag[i2],
-         atom->tag[i3],
-         atom->tag[i4],
-         eangle);
-  */
-
   // chain rule terms for first derivative components
 
   terma = -deddt / (rap2*rm);
diff --git a/tools/tinker2lmp.py b/tools/tinker2lmp.py
index 9b4aec3a5f..4d272c3fa9 100644
--- a/tools/tinker2lmp.py
+++ b/tools/tinker2lmp.py
@@ -580,9 +580,6 @@ atype = []
 aparams = []
 baparams = []
 
-#DEBUG
-opcount = 0
-
 for atom1,atom2,atom3 in alist:
   type1 = type[atom1-1]
   type2 = type[atom2-1]
@@ -609,9 +606,6 @@ for atom1,atom2,atom3 in alist:
     #     option 2 if one of 2 additional bonds is to an H atom
     #     option 3 if both of 2 additional bonds is to an H atom
 
-    # DEBUG
-    debugflag = 0
-    
     if len(params[3]) == 1:
       which = 1
       
@@ -638,10 +632,6 @@ for atom1,atom2,atom3 in alist:
       print "  Nbonds and hydrogen count:",nbonds,hcount
       print "  which:",which,m
 
-      # DEBUG
-      debugflag = 1
-      opcount += 1
-      
     elif len(params[3]) == 3:
       nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
       
@@ -668,24 +658,16 @@ for atom1,atom2,atom3 in alist:
       print "  Nbonds and hydrogen count:",nbonds,hcount
       print "  which:",which,m
 
-      # DEBUG
-      debugflag = 1
-      opcount += 1
-
   else:
     print "Angle not found",atom1,atom2,atom3,c1,c2,c3
     sys.exit()
     
   if not flags[m]:
     pflag,v1,v2,v3,v4,v5,v6 = params[3][which-1]
-    # DEBUG single line
-    if debugflag: pflag = 2
     aparams.append((pflag,v1,v2,v3,v4,v5,v6))
     flags[m] = len(aparams)
   atype.append(flags[m])
 
-print "OPTION angles",opcount
-
   # NOTE: baparams may need to be flipped if match is 3,2,1 instead of 1,2,3
   
   # NOTE: mismatch between angle and bondangle params may not be handled right