From 2c7badfa43f3b83485be9bd47e84f1234744ba72 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Wed, 13 Apr 2022 17:05:44 -0600
Subject: [PATCH] debug for improper amoeba

---
 examples/amoeba/in.ubiquitin   |  2 +-
 src/AMOEBA/angle_amoeba.cpp    | 12 ++++++------
 src/AMOEBA/improper_amoeba.cpp |  8 +++++++-
 3 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/examples/amoeba/in.ubiquitin b/examples/amoeba/in.ubiquitin
index 8870fcdf63..8abcac750f 100644
--- a/examples/amoeba/in.ubiquitin
+++ b/examples/amoeba/in.ubiquitin
@@ -29,7 +29,7 @@ read_data       data.ubiquitin fix amtype NULL "Tinker Types" &
                 fix pit pitorsions PiTorsions &
                 fix bit bitorsions BiTorsions
 
-pair_style      amoeba include angle
+pair_style      amoeba include improper
 pair_coeff      * * amoeba_ubiquitin.prm amoeba_ubiquitin.key
 
 special_bonds   lj/coul 0.5 0.5 0.5 one/five yes
diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp
index 6f4f72b432..fb59ecb39d 100644
--- a/src/AMOEBA/angle_amoeba.cpp
+++ b/src/AMOEBA/angle_amoeba.cpp
@@ -116,17 +116,17 @@ void AngleAmoeba::compute(int eflag, int vflag)
     if (enable_angle) {
       tflag = pflag[type];
 
-      //if (tflag && nspecial[i2][0] == 3)
-      //  tinker_anglep(i1,i2,i3,type,eflag);
-      //else
-      //  tinker_angle(i1,i2,i3,type,eflag);
+      if (tflag && nspecial[i2][0] == 3)
+        tinker_anglep(i1,i2,i3,type,eflag);
+      else
+        tinker_angle(i1,i2,i3,type,eflag);
 
       // bondangle = bond-stretch cross term in Tinker
-
+      
       if (ba_k1[type] != 0.0)
         tinker_bondangle(i1,i2,i3,type,eflag);
     }
-
+    
     // Urey-Bradley H-H bond term within water molecules
 
     if (enable_urey) {
diff --git a/src/AMOEBA/improper_amoeba.cpp b/src/AMOEBA/improper_amoeba.cpp
index d821c9411a..64c2f77030 100644
--- a/src/AMOEBA/improper_amoeba.cpp
+++ b/src/AMOEBA/improper_amoeba.cpp
@@ -197,7 +197,13 @@ void ImproperAmoeba::compute(int eflag, int vflag)
     fd[2] = dedcos * (dccdzid+deedzid);
     fb[0] = -fa[0] - fc[0] - fd[0];
     fb[1] = -fa[1] - fc[1] - fd[1];
-    fb[2] = -fa[1] - fc[2] - fd[2];
+    fb[2] = -fa[2] - fc[2] - fd[2];
+
+    printf("IMP DBAC %d %d %d %d: angle %g eng %g dedcos %g "
+           "fD %g %g %g fA %g %g %g fC %g %g %g\n",
+           atom->tag[id],atom->tag[ib],atom->tag[ia],atom->tag[ic],
+           angle,e,dedcos,fd[0],fd[1],fd[2],fa[0],fa[1],fa[2],
+           fc[0],fc[1],fc[2]);
 
     // apply force to each of 4 atoms