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

flip signs for forces instead of gradients

Browse Source
This commit is contained in:
Steve Plimpton
2022-04-08 15:23:17 -06:00
parent 3039d10742
commit 2d2660487d
8 changed files with 105 additions and 106 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/AMOEBA/amoeba_charge_transfer.cpp
View File

@ -139,12 +139,12 @@ void PairAmoeba::charge_transfer()
// increment the total charge transfer energy and derivatives
f[i][0] -= frcx;
f[i][1] -= frcy;
f[i][2] -= frcz;
f[j][0] += frcx;
f[j][1] += frcy;
f[j][2] += frcz;
f[i][0] += frcx;
f[i][1] += frcy;
f[i][2] += frcz;
f[j][0] -= frcx;
f[j][1] -= frcy;
f[j][2] -= frcz;
// increment the internal virial tensor components

18
src/AMOEBA/amoeba_dispersion.cpp
View File

@ -203,12 +203,12 @@ void PairAmoeba::dispersion_real()
dedx = de * xr;
dedy = de * yr;
dedz = de * zr;
f[i][0] += dedx;
f[i][1] += dedy;
f[i][2] += dedz;
f[j][0] -= dedx;
f[j][1] -= dedy;
f[j][2] -= dedz;
f[i][0] -= dedx;
f[i][1] -= dedy;
f[i][2] -= dedz;
f[j][0] += dedx;
f[j][1] += dedy;
f[j][2] += dedz;
// increment the internal virial tensor components
@ -407,9 +407,9 @@ void PairAmoeba::dispersion_kspace()
}
fi = csix[iclass];
f[m][0] += fi * (recip[0][0]*de1 + recip[0][1]*de2 + recip[0][2]*de3);
f[m][1] += fi * (recip[1][0]*de1 + recip[1][1]*de2 + recip[1][2]*de3);
f[m][2] += fi * (recip[2][0]*de1 + recip[2][1]*de2 + recip[2][2]*de3);
f[m][0] -= fi * (recip[0][0]*de1 + recip[0][1]*de2 + recip[0][2]*de3);
f[m][1] -= fi * (recip[1][0]*de1 + recip[1][1]*de2 + recip[1][2]*de3);
f[m][2] -= fi * (recip[2][0]*de1 + recip[2][1]*de2 + recip[2][2]*de3);
}
// account for the energy and virial correction terms

36
src/AMOEBA/amoeba_hal.cpp
View File

@ -159,31 +159,31 @@ void PairAmoeba::hal()
"ghost comm is too short");
if (i == iv) {
f[i][0] += dedx;
f[i][1] += dedy;
f[i][2] += dedz;
f[i][0] -= dedx;
f[i][1] -= dedy;
f[i][2] -= dedz;
} else {
f[i][0] += dedx*redi;
f[i][1] += dedy*redi;
f[i][2] += dedz*redi;
f[iv][0] += dedx*rediv;
f[iv][1] += dedy*rediv;
f[iv][2] += dedz*rediv;
f[i][0] -= dedx*redi;
f[i][1] -= dedy*redi;
f[i][2] -= dedz*redi;
f[iv][0] -= dedx*rediv;
f[iv][1] -= dedy*rediv;
f[iv][2] -= dedz*rediv;
}
if (j == jv) {
f[j][0] -= dedx;
f[j][1] -= dedy;
f[j][2] -= dedz;
f[j][0] += dedx;
f[j][1] += dedy;
f[j][2] += dedz;
} else {
redj = kred[jclass];
redjv = 1.0 - redj;
f[j][0] -= dedx*redj;
f[j][1] -= dedy*redj;
f[j][2] -= dedz*redj;
f[jv][0] -= dedx*redjv;
f[jv][1] -= dedy*redjv;
f[jv][2] -= dedz*redjv;
f[j][0] += dedx*redj;
f[j][1] += dedy*redj;
f[j][2] += dedz*redj;
f[jv][0] += dedx*redjv;
f[jv][1] += dedy*redjv;
f[jv][2] += dedz*redjv;
}
// increment the internal virial tensor components

19
src/AMOEBA/amoeba_multipole.cpp
View File

@ -509,18 +509,18 @@ void PairAmoeba::multipole_real()
// increment force-based gradient and torque on first site
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] += ttmi[0];
tq[i][1] += ttmi[1];
tq[i][2] += ttmi[2];
// increment force-based gradient and torque on second site
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] += ttmk[0];
tq[j][1] += ttmk[1];
tq[j][2] += ttmk[2];
@ -807,12 +807,11 @@ void PairAmoeba::multipole_kspace()
h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3; // matvec?
h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
f[i][0] += h1;
f[i][1] += h2;
f[i][2] += h3;
f[i][0] -= h1;
f[i][1] -= h2;
f[i][2] -= h3;
}
empole += 0.5*e;
//printf("mpole_force %g %g %g \n", f[0][0], f[0][1], f[0][2]);
// augment the permanent multipole virial contributions

36
src/AMOEBA/amoeba_polar.cpp
View File

@ -1140,12 +1140,12 @@ void PairAmoeba::polar_real()
// increment force-based gradient on the interaction sites
f[i][0] -= frcx;
f[i][1] -= frcy;
f[i][2] -= frcz;
f[j][0] += frcx;
f[j][1] += frcy;
f[j][2] += frcz;
f[i][0] += frcx;
f[i][1] += frcy;
f[i][2] += frcz;
f[j][0] -= frcx;
f[j][1] -= frcy;
f[j][2] -= frcz;
// increment the virial due to pairwise Cartesian forces
@ -1518,9 +1518,9 @@ void PairAmoeba::polar_kspace()
h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;
h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
f[i][0] += h1;
f[i][1] += h2;
f[i][2] += h3;
f[i][0] -= h1;
f[i][1] -= h2;
f[i][2] -= h3;
}
// set the potential to be the induced dipole average
@ -1672,9 +1672,9 @@ void PairAmoeba::polar_kspace()
h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
f[i][0] += copm[k+m+1]*h1;
f[i][1] += copm[k+m+1]*h2;
f[i][2] += copm[k+m+1]*h3;
f[i][0] -= copm[k+m+1]*h1;
f[i][1] -= copm[k+m+1]*h2;
f[i][2] -= copm[k+m+1]*h3;
for (j = 1; j < 4; j++) {
cphid[j] = 0.0;
@ -1762,9 +1762,9 @@ void PairAmoeba::polar_kspace()
h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3;
h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
f[i][0] += h1;
f[i][1] += h2;
f[i][2] += h3;
f[i][0] -= h1;
f[i][1] -= h2;
f[i][2] -= h3;
for (j = 1; j < 4; j++) {
cphid[j] = 0.0;
@ -1838,9 +1838,9 @@ void PairAmoeba::polar_kspace()
h1 = recip[0][0]*f1 + recip[0][1]*f2 + recip[0][2]*f3; // matvec
h2 = recip[1][0]*f1 + recip[1][1]*f2 + recip[1][2]*f3;
h3 = recip[2][0]*f1 + recip[2][1]*f2 + recip[2][2]*f3;
f[i][0] += h1;
f[i][1] += h2;
f[i][2] += h3;
f[i][0] -= h1;
f[i][1] -= h2;
f[i][2] -= h3;
for (j = 1; j < 4; j++) {
cphid[j] = 0.0;

24
src/AMOEBA/angle_amoeba.cpp
View File

@ -393,27 +393,27 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
// 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];
f[i1][0] -= f1[0];
f[i1][1] -= f1[1];
f[i1][2] -= f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] += f2[0];
f[i2][1] += f2[1];
f[i2][2] += f2[2];
f[i2][0] -= f2[0];
f[i2][1] -= f2[1];
f[i2][2] -= f2[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
f[i3][0] -= f3[0];
f[i3][1] -= f3[1];
f[i3][2] -= f3[2];
}
if (newton_bond || i4 < nlocal) {
f[i4][0] += f4[0];
f[i4][1] += f4[1];
f[i4][2] += f4[2];
f[i4][0] -= f4[0];
f[i4][1] -= f4[1];
f[i4][2] -= f4[2];
}
if (evflag) ev_tally4(i1,i2,i3,14,nlocal,newton_bond,eangle,f1,f2,f3,f4);

30
src/AMOEBA/fix_amoeba_bitorsion.cpp
View File

@ -605,37 +605,37 @@ void FixAmoebaBiTorsion::post_force(int vflag)
if (ia < nlocal) {
ebitorsion += engfraction;
f[ia][0] += dedxia;
f[ia][1] += dedyia;
f[ia][2] += dedzia;
f[ia][0] -= dedxia;
f[ia][1] -= dedyia;
f[ia][2] -= dedzia;
}
if (ib < nlocal) {
ebitorsion += engfraction;
f[ib][0] += dedxib + dedxib2;
f[ib][1] += dedyib + dedyib2;
f[ib][2] += dedzib + dedzib2;
f[ib][0] -= dedxib + dedxib2;
f[ib][1] -= dedyib + dedyib2;
f[ib][2] -= dedzib + dedzib2;
}
if (ic < nlocal) {
ebitorsion += engfraction;
f[ic][0] += dedxic + dedxic2;
f[ic][1] += dedyic + dedyic2;
f[ic][2] += dedzic + dedzic2;
f[ic][0] -= dedxic + dedxic2;
f[ic][1] -= dedyic + dedyic2;
f[ic][2] -= dedzic + dedzic2;
}
if (id < nlocal) {
ebitorsion += engfraction;
f[id][0] += dedxid + dedxid2;
f[id][1] += dedyid + dedyid2;
f[id][2] += dedzid + dedzid2;
f[id][0] -= dedxid + dedxid2;
f[id][1] -= dedyid + dedyid2;
f[id][2] -= dedzid + dedzid2;
}
if (ie < nlocal) {
ebitorsion += engfraction;
f[ie][0] += dedxie2;
f[ie][1] += dedyie2;
f[ie][2] += dedzie2;
f[ie][0] -= dedxie2;
f[ie][1] -= dedyie2;
f[ie][2] -= dedzie2;
}
// increment the internal virial tensor components

36
src/AMOEBA/fix_amoeba_pitorsion.cpp
View File

@ -499,44 +499,44 @@ void FixAmoebaPiTorsion::post_force(int vflag)
if (ia < nlocal) {
epitorsion += engfraction;
f[ia][0] += dedxia;
f[ia][1] += dedyia;
f[ia][2] += dedzia;
f[ia][0] -= dedxia;
f[ia][1] -= dedyia;
f[ia][2] -= dedzia;
}
if (ib < nlocal) {
epitorsion += engfraction;
f[ib][0] += dedxib;
f[ib][1] += dedyib;
f[ib][2] += dedzib;
f[ib][0] -= dedxib;
f[ib][1] -= dedyib;
f[ib][2] -= dedzib;
}
if (ic < nlocal) {
epitorsion += engfraction;
f[ic][0] += dedxic;
f[ic][1] += dedyic;
f[ic][2] += dedzic;
f[ic][0] -= dedxic;
f[ic][1] -= dedyic;
f[ic][2] -= dedzic;
}
if (id < nlocal) {
epitorsion += engfraction;
f[id][0] += dedxid;
f[id][1] += dedyid;
f[id][2] += dedzid;
f[id][0] -= dedxid;
f[id][1] -= dedyid;
f[id][2] -= dedzid;
}
if (ie < nlocal) {
epitorsion += engfraction;
f[ie][0] += dedxie;
f[ie][1] += dedyie;
f[ie][2] += dedzie;
f[ie][0] -= dedxie;
f[ie][1] -= dedyie;
f[ie][2] -= dedzie;
}
if (ig < nlocal) {
epitorsion += engfraction;
f[ig][0] += dedxig;
f[ig][1] += dedyig;
f[ig][2] += dedzig;
f[ig][0] -= dedxig;
f[ig][1] -= dedyig;
f[ig][2] -= dedzig;
}
// virial tensor components