diff --git a/src/AMOEBA/amoeba_repulsion.cpp b/src/AMOEBA/amoeba_repulsion.cpp index 2823a82de2..0342eeb513 100644 --- a/src/AMOEBA/amoeba_repulsion.cpp +++ b/src/AMOEBA/amoeba_repulsion.cpp @@ -334,22 +334,20 @@ void PairAmoeba::repulsion() erepulse += e; - //fprintf(fp,"REPUL %d %d %15.12g %15.12g\n",atom->tag[i],atom->tag[j],r,e); - // increment force-based gradient and torque on atom I - f[i][0] += frcx; - f[i][1] += frcy; - f[i][2] += frcz; + f[i][0] -= frcx; + f[i][1] -= frcy; + f[i][2] -= frcz; tq[i][0] += ttri[0]; tq[i][1] += ttri[1]; tq[i][2] += ttri[2]; // increment force-based gradient and torque on atom J - f[j][0] -= frcx; - f[j][1] -= frcy; - f[j][2] -= frcz; + f[j][0] += frcx; + f[j][1] += frcy; + f[j][2] += frcz; tq[j][0] += ttrk[0]; tq[j][1] += ttrk[1]; tq[j][2] += ttrk[2]; diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp index e5831f9348..422b905949 100644 --- a/src/AMOEBA/amoeba_utils.cpp +++ b/src/AMOEBA/amoeba_utils.cpp @@ -1016,8 +1016,8 @@ void PairAmoeba::torque2force(int i, double *trq, if (axetyp == ZONLY) { for (j = 0; j < 3; j++) { du = uv[j]*dphidv/(usiz*uvsin) + uw[j]*dphidw/usiz; - f[ia][j] -= du; - f[ib][j] += du; + f[ia][j] += du; + f[ib][j] -= du; frcz[j] += du; } @@ -1027,9 +1027,9 @@ void PairAmoeba::torque2force(int i, double *trq, for (j = 0; j < 3; j++) { du = uv[j]*dphidv/(usiz*uvsin) + uw[j]*dphidw/usiz; dv = -uv[j]*dphidu/(vsiz*uvsin); - f[ia][j] -= du; - f[ic][j] -= dv; - f[ib][j] += du + dv; + f[ia][j] += du; + f[ic][j] += dv; + f[ib][j] -= du + dv; frcz[j] += du; frcx[j] += dv; } @@ -1040,9 +1040,9 @@ void PairAmoeba::torque2force(int i, double *trq, for (j = 0; j < 3; j++) { du = uv[j]*dphidv/(usiz*uvsin) + 0.5*uw[j]*dphidw/usiz; dv = -uv[j]*dphidu/(vsiz*uvsin) + 0.5*vw[j]*dphidw/vsiz; - f[ia][j] -= du; - f[ic][j] -= dv; - f[ib][j] += du + dv; + f[ia][j] += du; + f[ic][j] += dv; + f[ib][j] -= du + dv; frcz[j] += du; frcx[j] += dv; } @@ -1054,10 +1054,10 @@ void PairAmoeba::torque2force(int i, double *trq, du = ur[j]*dphidr/(usiz*ursin) + us[j]*dphids/usiz; dv = (vssin*s[j]-vscos*t1[j])*dphidu / (vsiz*(ut1sin+ut2sin)); dw = (wssin*s[j]-wscos*t2[j])*dphidu / (wsiz*(ut1sin+ut2sin)); - f[ia][j] -= du; - f[ic][j] -= dv; - f[id][j] -= dw; - f[ib][j] += du + dv + dw; + f[ia][j] += du; + f[ic][j] += dv; + f[id][j] += dw; + f[ib][j] -= du + dv + dw; frcz[j] += du; frcx[j] += dv; frcy[j] += dw; @@ -1096,8 +1096,8 @@ void PairAmoeba::torque2force(int i, double *trq, eps[2] = del[0]*w[1] - del[1]*w[0]; for (j = 0; j < 3; j++) { dw = del[j]*dphidr/(wsiz*rwsin) + eps[j]*dphiddel*wpcos/(wsiz*psiz); - f[id][j] -= dw; - f[ib][j] += dw; + f[id][j] += dw; + f[ib][j] -= dw; frcy[j] += dw; } r[0] = v[0] + w[0]; @@ -1121,8 +1121,8 @@ void PairAmoeba::torque2force(int i, double *trq, eps[2] = del[0]*u[1] - del[1]*u[0]; for (j = 0; j < 3; j++) { du = del[j]*dphidr/(usiz*rusin) + eps[j]*dphiddel*upcos/(usiz*psiz); - f[ia][j] -= du; - f[ib][j] += du; + f[ia][j] += du; + f[ib][j] -= du; frcz[j] += du; } r[0] = u[0] + w[0]; @@ -1146,8 +1146,8 @@ void PairAmoeba::torque2force(int i, double *trq, eps[2] = del[0]*v[1] - del[1]*v[0]; for (j = 0; j < 3; j++) { dv = del[j]*dphidr/(vsiz*rvsin) + eps[j]*dphiddel*vpcos/(vsiz*psiz); - f[ic][j] -= dv; - f[ib][j] += dv; + f[ic][j] += dv; + f[ib][j] -= dv; frcx[j] += dv; } }