diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp
index 6faf2edbdd..57079cdd83 100644
--- a/src/AMOEBA/angle_amoeba.cpp
+++ b/src/AMOEBA/angle_amoeba.cpp
@@ -41,6 +41,11 @@ AngleAmoeba::AngleAmoeba(LAMMPS *lmp) : Angle(lmp)
   k4 = nullptr;
   k5 = nullptr;
   k6 = nullptr;
+
+  ba_k1 = nullptr;
+  ba_k2 = nullptr;
+  ba_r1 = nullptr;
+  ba_r2 = nullptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -51,6 +56,8 @@ AngleAmoeba::~AngleAmoeba()
 
   if (allocated) {
     memory->destroy(setflag);
+    memory->destroy(setflag_a);
+    memory->destroy(setflag_ba);
 
     memory->destroy(pflag);
     memory->destroy(theta0);
@@ -59,6 +66,11 @@ AngleAmoeba::~AngleAmoeba()
     memory->destroy(k4);
     memory->destroy(k5);
     memory->destroy(k6);
+
+    memory->destroy(ba_k1);
+    memory->destroy(ba_k2);
+    memory->destroy(ba_r1);
+    memory->destroy(ba_r2);
   }
 }
 
@@ -72,6 +84,7 @@ void AngleAmoeba::compute(int eflag, int vflag)
   double dtheta,dtheta2,dtheta3,dtheta4,dtheta5,dtheta6,de_angle;
   double dr1,dr2,tk1,tk2,aa1,aa2,aa11,aa12,aa21,aa22;
   double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22,b1,b2;
+  double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22;
 
   eangle = 0.0;
   ev_init(eflag,vflag);
@@ -159,6 +172,50 @@ void AngleAmoeba::compute(int eflag, int vflag)
     if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + 
                  k4[type]*dtheta4 + k5[type]*dtheta5 + k6[type]*dtheta6;
 
+    // force & energy for bond-angle term
+    // bond-stretch cross term in Tinker
+
+    dr1 = r1 - ba_r1[type];
+    dr2 = r2 - ba_r2[type];
+
+    aa1 = s * dr1 * ba_k1[type];
+    aa2 = s * dr2 * ba_k2[type];
+
+    aa11 = aa1 * c / rsq1;
+    aa12 = -aa1 / (r1 * r2);
+    aa21 = aa2 * c / rsq1;
+    aa22 = -aa2 / (r1 * r2);
+
+    vx11 = (aa11 * delx1) + (aa12 * delx2);
+    vx12 = (aa21 * delx1) + (aa22 * delx2);
+    vy11 = (aa11 * dely1) + (aa12 * dely2);
+    vy12 = (aa21 * dely1) + (aa22 * dely2);
+    vz11 = (aa11 * delz1) + (aa12 * delz2);
+    vz12 = (aa21 * delz1) + (aa22 * delz2);
+
+    aa11 = aa1 * c / rsq2;
+    aa21 = aa2 * c / rsq2;
+
+    vx21 = (aa11 * delx2) + (aa12 * delx1);
+    vx22 = (aa21 * delx2) + (aa22 * delx1);
+    vy21 = (aa11 * dely2) + (aa12 * dely1);
+    vy22 = (aa21 * dely2) + (aa22 * dely1);
+    vz21 = (aa11 * delz2) + (aa12 * delz1);
+    vz22 = (aa21 * delz2) + (aa22 * delz1);
+
+    b1 = ba_k1[type] * dtheta / r1;
+    b2 = ba_k2[type] * dtheta / r2;
+
+    f1[0] -= vx11 + b1*delx1 + vx12;
+    f1[1] -= vy11 + b1*dely1 + vy12;
+    f1[2] -= vz11 + b1*delz1 + vz12;
+
+    f3[0] -= vx21 + b2*delx2 + vx22;
+    f3[1] -= vy21 + b2*dely2 + vy22;
+    f3[2] -= vz21 + b2*delz2 + vz22;
+
+    if (eflag) eangle += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
+
     // apply force to each of 3 atoms
 
     if (newton_bond || i1 < nlocal) {
@@ -382,8 +439,13 @@ void AngleAmoeba::allocate()
   memory->create(k5,n+1,"angle:k5");
   memory->create(k6,n+1,"angle:k6");
 
+  memory->create(ba_k1,n+1,"angle:ba_k1");
+  memory->create(ba_k2,n+1,"angle:ba_k2");
+  memory->create(ba_r1,n+1,"angle:ba_r1");
+  memory->create(ba_r2,n+1,"angle:ba_r2");
+
   memory->create(setflag,n+1,"angle:setflag");
-  for (int i = 1; i <= n; i++) setflag[i] = 0;
+  for (int i = 1; i <= n; i++) setflag[i] = setflag_a[i] = setflag_ba[i] = 0;
 }
 
 /* ----------------------------------------------------------------------
@@ -392,6 +454,7 @@ void AngleAmoeba::allocate()
 
 void AngleAmoeba::coeff(int narg, char **arg)
 {
+  if (narg < 2) error->all(FLERR,"Incorrect args for angle coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi;
@@ -399,32 +462,52 @@ void AngleAmoeba::coeff(int narg, char **arg)
 
   int count = 0;
   
-  if (narg != 8) error->all(FLERR,"Incorrect args for angle coefficients");
+  if (strcmp(arg[1],"ba") == 0) {
+    if (narg != 6) error->all(FLERR,"Incorrect args for angle coefficients");
 
-  int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
-  double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
-  double k2_one = utils::numeric(FLERR,arg[3],false,lmp);
-  double k3_one = utils::numeric(FLERR,arg[4],false,lmp);
-  double k4_one = utils::numeric(FLERR,arg[5],false,lmp);
-  double k5_one = utils::numeric(FLERR,arg[6],false,lmp);
-  double k6_one = utils::numeric(FLERR,arg[7],false,lmp);
+    double ba_k1_one = utils::numeric(FLERR,arg[2],false,lmp);
+    double ba_k2_one = utils::numeric(FLERR,arg[3],false,lmp);
+    double ba_r1_one = utils::numeric(FLERR,arg[4],false,lmp);
+    double ba_r2_one = utils::numeric(FLERR,arg[5],false,lmp);
 
-  // convert theta0 from degrees to radians
+    for (int i = ilo; i <= ihi; i++) {
+      ba_k1[i] = ba_k1_one;
+      ba_k2[i] = ba_k2_one;
+      ba_r1[i] = ba_r1_one;
+      ba_r2[i] = ba_r2_one;
+      setflag_ba[i] = 1;
+      count++;
+    }
+
+  } else { 
+    if (narg != 8) error->all(FLERR,"Incorrect args for angle coefficients");
+
+    int pflag_one = utils::inumeric(FLERR,arg[1],false,lmp);
+    double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
+    double k2_one = utils::numeric(FLERR,arg[3],false,lmp);
+    double k3_one = utils::numeric(FLERR,arg[4],false,lmp);
+    double k4_one = utils::numeric(FLERR,arg[5],false,lmp);
+    double k5_one = utils::numeric(FLERR,arg[6],false,lmp);
+    double k6_one = utils::numeric(FLERR,arg[7],false,lmp);
+
+    // convert theta0 from degrees to radians
   
-  for (int i = ilo; i <= ihi; i++) {
-    pflag[i] = pflag_one;
-    theta0[i] = theta0_one/180.0 * MY_PI;
-    k2[i] = k2_one;
-    k3[i] = k3_one;
-    k4[i] = k4_one;
-    k5[i] = k5_one;
-    k6[i] = k6_one;
-    count++;
+    for (int i = ilo; i <= ihi; i++) {
+      pflag[i] = pflag_one;
+      theta0[i] = theta0_one/180.0 * MY_PI;
+      k2[i] = k2_one;
+      k3[i] = k3_one;
+      k4[i] = k4_one;
+      k5[i] = k5_one;
+      k6[i] = k6_one;
+      count++;
+    }
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
 
-  for (int i = ilo; i <= ihi; i++) setflag[i] = 1;
+  for (int i = ilo; i <= ihi; i++)
+    if (setflag_a[i] == 1 && setflag_ba[i] == 1) setflag[i] = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -447,6 +530,11 @@ void AngleAmoeba::write_restart(FILE *fp)
   fwrite(&k4[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&k5[1],sizeof(double),atom->nangletypes,fp);
   fwrite(&k6[1],sizeof(double),atom->nangletypes,fp);
+
+  fwrite(&ba_k1[1],sizeof(double),atom->nangletypes,fp);
+  fwrite(&ba_k2[1],sizeof(double),atom->nangletypes,fp);
+  fwrite(&ba_r1[1],sizeof(double),atom->nangletypes,fp);
+  fwrite(&ba_r2[1],sizeof(double),atom->nangletypes,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -465,6 +553,11 @@ void AngleAmoeba::read_restart(FILE *fp)
     utils::sfread(FLERR,&k4[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
     utils::sfread(FLERR,&k5[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
     utils::sfread(FLERR,&k6[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
+
+    utils::sfread(FLERR,&ba_k1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
+    utils::sfread(FLERR,&ba_k2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
+    utils::sfread(FLERR,&ba_r1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
+    utils::sfread(FLERR,&ba_r2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
   }
 
   MPI_Bcast(&pflag[1],atom->nangletypes,MPI_INT,0,world);
@@ -475,6 +568,11 @@ void AngleAmoeba::read_restart(FILE *fp)
   MPI_Bcast(&k5[1],atom->nangletypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&k6[1],atom->nangletypes,MPI_DOUBLE,0,world);
 
+  MPI_Bcast(&ba_k1[1],atom->nangletypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&ba_k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&ba_r1[1],atom->nangletypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&ba_r2[1],atom->nangletypes,MPI_DOUBLE,0,world);
+
   for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
 }
 
@@ -487,6 +585,10 @@ void AngleAmoeba::write_data(FILE *fp)
   for (int i = 1; i <= atom->nangletypes; i++)
     fprintf(fp,"%d %d %g %g %g %g %g %g\n",
             i,pflag[i],theta0[i]/MY_PI*180.0,k2[i],k3[i],k4[i],k5[i],k6[i]);
+
+  fprintf(fp,"\nBondAngle Coeffs\n\n");
+  for (int i = 1; i <= atom->nangletypes; i++)
+    fprintf(fp,"%d %g %g %g %g\n",i,ba_k1[i],ba_k2[i],ba_r1[i],ba_r2[i]);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -526,5 +628,9 @@ double AngleAmoeba::single(int type, int i1, int i2, int i3)
   double energy = k2[type]*dtheta2 + k3[type]*dtheta3 + k4[type]*dtheta4
            + k5[type]*dtheta5 + k6[type]*dtheta6;
 
+  double dr1 = r1 - ba_r1[type];
+  double dr2 = r2 - ba_r2[type];
+  energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
+
   return energy;
 }
diff --git a/src/AMOEBA/angle_amoeba.h b/src/AMOEBA/angle_amoeba.h
index 66a0c2c6b7..0fa50880e9 100644
--- a/src/AMOEBA/angle_amoeba.h
+++ b/src/AMOEBA/angle_amoeba.h
@@ -39,6 +39,8 @@ class AngleAmoeba : public Angle {
  protected:
   int *pflag;
   double *theta0, *k2, *k3, *k4, *k5, *k6;
+  double *ba_k1, *ba_k2, *ba_r1, *ba_r2;
+  int *setflag_a, *setflag_ba;
 
   void anglep(int, int, int, int, int);
   void allocate();
diff --git a/src/AMOEBA/improper_amoeba.cpp b/src/AMOEBA/improper_amoeba.cpp
index d32b23ae2a..9741939241 100644
--- a/src/AMOEBA/improper_amoeba.cpp
+++ b/src/AMOEBA/improper_amoeba.cpp
@@ -17,9 +17,10 @@
 #include <cmath>
 #include "atom.h"
 #include "comm.h"
+#include "update.h"
 #include "neighbor.h"
 #include "force.h"
-#include "update.h"
+#include "pair.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
@@ -56,7 +57,13 @@ void ImproperAmoeba::compute(int eflag, int vflag)
   double xab,yab,zab,xcb,ycb,zcb,xdb,ydb,zdb,xad,yad,zad,xcd,ycd,zcd;
   double rad2,rcd2,rdb2,dot,cc,ee;
   double sine,angle;
-  double eimproper,f1[3],f2[3],f3[3],f4[3];
+  double dt,dt2,dt3,dt4,e;
+  double deddt,sign,dedcos,term;
+  double dccdxia,dccdyia,dccdzia,dccdxic,dccdyic,dccdzic;
+  double dccdxid,dccdyid,dccdzid;
+  double deedxia,deedyia,deedzia,deedxic,deedyic,deedzic;
+  double deedxid,deedyid,deedzid;
+  double eimproper,fa[3],fb[3],fc[3],fd[3];
 
   eimproper = 0.0;
   ev_init(eflag,vflag);
@@ -129,48 +136,91 @@ void ImproperAmoeba::compute(int eflag, int vflag)
 
     sine = fabs(ee) / sqrt(cc*rdb2);
     sine = MIN(1.0,sine);
-    angle = asin(sine) * MY_PI/180.0;
+    angle = 180.0/MY_PI * asin(sine);   // angle = degrees
 
-    //dt = angle;
-    //dt2 = dt * dt;
-    //dt3 = dt2 * dt;
-    //dt4 = dt2 * dt2;
-    //e = opbunit * force * dt2 * (1.0d0+copb*dt+qopb*dt2+popb*dt3+sopb*dt4);
-    //deddt = opbunit * force * dt * radian * 
-    // (2.0d0 + 3.0d0*copb*dt + 4.0d0*qopb*dt2 + 5.0d0*popb*dt3 + 6.0d0*sopb*dt4);
-    //dedcos = -deddt * sign(1.0d0,ee) / sqrt(cc*rdb2-ee*ee);
+    dt = angle;
+    dt2 = dt * dt;
+    dt3 = dt2 * dt;
+    dt4 = dt2 * dt2;
+    double prefactor = 1.0 / (180.0/MY_PI) / (180.0/MY_PI);
+    e = prefactor * k[type] * dt2 * (1.0 + opbend_cubic*dt + opbend_quartic*dt2 + 
+                         opbend_pentic*dt3 + opbend_sextic*dt4);
+    eimproper += e;
+
+    deddt = k[type] * dt * 
+      (2.0 + 3.0*opbend_cubic*dt + 4.0*opbend_quartic*dt2 + 
+       5.0*opbend_pentic*dt3 + 6.0*opbend_sextic*dt4);
+    sign = (ee >= 0.0) ? 1.0 : -1.0;
+    dedcos = -deddt * sign / sqrt(cc*rdb2 - ee*ee);
+
+    // chain rule terms for first derivative components
+
+    term = ee / cc;
+    dccdxia = (xad*rcd2-xcd*dot) * term;
+    dccdyia = (yad*rcd2-ycd*dot) * term;
+    dccdzia = (zad*rcd2-zcd*dot) * term;
+    dccdxic = (xcd*rad2-xad*dot) * term;
+    dccdyic = (ycd*rad2-yad*dot) * term;
+    dccdzic = (zcd*rad2-zad*dot) * term;
+    dccdxid = -dccdxia - dccdxic;
+    dccdyid = -dccdyia - dccdyic;
+    dccdzid = -dccdzia - dccdzic;
+
+    term = ee / rdb2;
+    deedxia = ydb*zcb - zdb*ycb;
+    deedyia = zdb*xcb - xdb*zcb;
+    deedzia = xdb*ycb - ydb*xcb;
+    deedxic = yab*zdb - zab*ydb;
+    deedyic = zab*xdb - xab*zdb;
+    deedzic = xab*ydb - yab*xdb;
+    deedxid = ycb*zab - zcb*yab + xdb*term;
+    deedyid = zcb*xab - xcb*zab + ydb*term;
+    deedzid = xcb*yab - ycb*xab + zdb*term;
+
+    //  compute first derivative components for this angle
+
+    fa[0] = dedcos * (dccdxia+deedxia);
+    fa[1] = dedcos * (dccdyia+deedyia);
+    fa[2] = dedcos * (dccdzia+deedzia);
+    fc[0] = dedcos * (dccdxic+deedxic);
+    fc[1] = dedcos * (dccdyic+deedyic);
+    fc[2] = dedcos * (dccdzic+deedzic);
+    fd[0] = dedcos * (dccdxid+deedxid);
+    fd[1] = dedcos * (dccdyid+deedyid);
+    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];
 
-    /*
     // apply force to each of 4 atoms
 
-    if (newton_bond || i1 < nlocal) {
-      f[i1][0] += f1[0];
-      f[i1][1] += f1[1];
-      f[i1][2] += f1[2];
+    if (newton_bond || id < nlocal) {
+      f[id][0] += fd[0];
+      f[id][1] += fd[1];
+      f[id][2] += fd[2];
     }
 
-    if (newton_bond || i2 < nlocal) {
-      f[i2][0] += f2[0];
-      f[i2][1] += f2[1];
-      f[i2][2] += f2[2];
+    if (newton_bond || ib < nlocal) {
+      f[ib][0] += fb[0];
+      f[ib][1] += fb[1];
+      f[ib][2] += fb[2];
     }
 
-    if (newton_bond || i3 < nlocal) {
-      f[i3][0] += f3[0];
-      f[i3][1] += f3[1];
-      f[i3][2] += f3[2];
+    if (newton_bond || ia < nlocal) {
+      f[ia][0] += fa[0];
+      f[ia][1] += fa[1];
+      f[ia][2] += fa[2];
     }
 
-    if (newton_bond || i4 < nlocal) {
-      f[i4][0] += f4[0];
-      f[i4][1] += f4[1];
-      f[i4][2] += f4[2];
+    if (newton_bond || ic < nlocal) {
+      f[ic][0] += fc[0];
+      f[ic][1] += fc[1];
+      f[ic][2] += fc[2];
     }
 
     if (evflag)
-      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
-               vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
-    */
+      ev_tally(id,ib,ia,ic,nlocal,newton_bond,e,fd,fa,fc,
+               xdb,ydb,zdb,xab,yab,zab,xic-xia,yic-yia,zic-zia);
   }
 }
 
@@ -213,6 +263,27 @@ void ImproperAmoeba::coeff(int narg, char **arg)
   if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
 }
 
+/* ----------------------------------------------------------------------
+   set opbend higher-order term weights from PairAmoeba
+------------------------------------------------------------------------- */
+
+void ImproperAmoeba::init_style()
+{
+  Pair *pair = force->pair_match("amoeba",1,0);
+  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);
+}
+
 /* ----------------------------------------------------------------------
    proc 0 writes out coeffs to restart file
 ------------------------------------------------------------------------- */
diff --git a/src/AMOEBA/improper_amoeba.h b/src/AMOEBA/improper_amoeba.h
index 439a6c39fb..b25a60ac47 100644
--- a/src/AMOEBA/improper_amoeba.h
+++ b/src/AMOEBA/improper_amoeba.h
@@ -30,11 +30,13 @@ class ImproperAmoeba : public Improper {
   ~ImproperAmoeba();
   void compute(int, int);
   void coeff(int, char **);
+  void init_style();
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_data(FILE *);
 
  protected:
+  double opbend_cubic,opbend_quartic,opbend_pentic,opbend_sextic;
   double *k;
 
   virtual void allocate();
diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp
index f9e098e884..474a764f22 100644
--- a/src/AMOEBA/pair_amoeba.cpp
+++ b/src/AMOEBA/pair_amoeba.cpp
@@ -2107,6 +2107,18 @@ void PairAmoeba::zero_energy_force_virial()
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+void *PairAmoeba::extract(const char *str, int &dim)
+{
+  dim = 0;
+  if (strcmp(str,"opbend_cubic") == 0) return (void *) &opbend_cubic;
+  if (strcmp(str,"opbend_quartic") == 0) return (void *) &opbend_quartic;
+  if (strcmp(str,"opbend_pentic") == 0) return (void *) &opbend_pentic;
+  if (strcmp(str,"opbend_sextic") == 0) return (void *) &opbend_sextic;
+  return nullptr;
+}
+
 /* ----------------------------------------------------------------------
   grow local vectors and arrays if necessary
   keep them atom->nmax in length
diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h
index b28a00fb84..607c261695 100644
--- a/src/AMOEBA/pair_amoeba.h
+++ b/src/AMOEBA/pair_amoeba.h
@@ -57,6 +57,8 @@ class PairAmoeba : public Pair {
   void pack_reverse_grid(int, void *, int, int *);
   void unpack_reverse_grid(int, void *, int, int *);
 
+  void *extract(const char *, int &);
+
  protected:
   int me,nprocs;
   int nmax;        // allocation for owned+ghost