ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

more force flips in torque2force

Browse Source
This commit is contained in:
Steve Plimpton
2022-04-08 15:41:08 -06:00
parent 2d2660487d
commit 2111797ed8
1 changed files with 18 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
src/AMOEBA/amoeba_utils.cpp
View File

@ -1111,8 +1111,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;
}
@ -1122,9 +1122,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;
}
@ -1135,9 +1135,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;
}
@ -1149,10 +1149,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;
@ -1191,8 +1191,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];
@ -1216,8 +1216,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];
@ -1241,8 +1241,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;
}
}