From 578a9ab1611ff700ce59ef38b142c92fdfac8a7d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 13 Apr 2022 14:52:01 -0600 Subject: [PATCH] bug fix for bondangle term --- src/AMOEBA/angle_amoeba.cpp | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp index 6ada36ea2a..6f4f72b432 100644 --- a/src/AMOEBA/angle_amoeba.cpp +++ b/src/AMOEBA/angle_amoeba.cpp @@ -500,21 +500,17 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag) 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; + 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; + f3[0] = -(vx21 + b2*delx2 + vx22); + f3[1] = -(vy21 + b2*dely2 + vy22); + f3[2] = -(vz21 + b2*delz2 + vz22); eangle = 0.0; if (eflag) eangle = ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta; - printf("BA: ijk %d %d %d eng %g fi %g %g %g fk %g %g %g\n", - atom->tag[i1],atom->tag[i2],atom->tag[i3], - eangle,f1[0],f1[1],f1[2],f3[0],f3[1],f3[2]); - // apply force to each of 3 atoms if (newton_bond || i1 < nlocal) {