sign flips for virial terms

This commit is contained in:
Steve Plimpton
2022-04-28 17:12:01 -06:00
parent f7cdfdd884
commit 221142a36d
9 changed files with 105 additions and 97 deletions

View File

@ -155,12 +155,12 @@ void PairAmoeba::charge_transfer()
vyz = zr * frcy;
vzz = zr * frcz;
virqxfer[0] += vxx;
virqxfer[1] += vyy;
virqxfer[2] += vzz;
virqxfer[3] += vxy;
virqxfer[4] += vxz;
virqxfer[5] += vyz;
virqxfer[0] -= vxx;
virqxfer[1] -= vyy;
virqxfer[2] -= vzz;
virqxfer[3] -= vxy;
virqxfer[4] -= vxz;
virqxfer[5] -= vyz;
// energy = e
// virial = 6-vec vir

View File

@ -219,12 +219,12 @@ void PairAmoeba::dispersion_real()
vzy = zr * dedy;
vzz = zr * dedz;
virdisp[0] += vxx;
virdisp[1] += vyy;
virdisp[2] += vzz;
virdisp[3] += vyx;
virdisp[4] += vzx;
virdisp[5] += vzy;
virdisp[0] -= vxx;
virdisp[1] -= vyy;
virdisp[2] -= vzz;
virdisp[3] -= vyx;
virdisp[4] -= vzx;
virdisp[5] -= vzy;
// energy = e
// virial = 6-vec vir
@ -418,8 +418,8 @@ void PairAmoeba::dispersion_kspace()
if (me == 0) {
edisp -= term;
virdisp[0] += term;
virdisp[1] += term;
virdisp[2] += term;
virdisp[0] -= term;
virdisp[1] -= term;
virdisp[2] -= term;
}
}

View File

@ -195,12 +195,12 @@ void PairAmoeba::hal()
vzy = zr * dedy;
vzz = zr * dedz;
virhal[0] += vxx;
virhal[1] += vyy;
virhal[2] += vzz;
virhal[3] += vyx;
virhal[4] += vzx;
virhal[5] += vzy;
virhal[0] -= vxx;
virhal[1] -= vyy;
virhal[2] -= vzz;
virhal[3] -= vyx;
virhal[4] -= vzx;
virhal[5] -= vzy;
// energy = e
// virial = 6-vec vir

View File

@ -533,12 +533,12 @@ void PairAmoeba::multipole_real()
vyz = -0.5 * (zr*frcy+yr*frcz);
vzz = -zr * frcz;
virmpole[0] += vxx;
virmpole[1] += vyy;
virmpole[2] += vzz;
virmpole[3] += vxy;
virmpole[4] += vxz;
virmpole[5] += vyz;
virmpole[0] -= vxx;
virmpole[1] -= vyy;
virmpole[2] -= vzz;
virmpole[3] -= vxy;
virmpole[4] -= vxz;
virmpole[5] -= vyz;
// energy = e
// virial = 6-vec vir
@ -583,12 +583,12 @@ void PairAmoeba::multipole_real()
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
vzz = zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
virmpole[0] += vxx;
virmpole[1] += vyy;
virmpole[2] += vzz;
virmpole[3] += vxy;
virmpole[4] += vxz;
virmpole[5] += vyz;
virmpole[0] -= vxx;
virmpole[1] -= vyy;
virmpole[2] -= vzz;
virmpole[3] -= vxy;
virmpole[4] -= vxz;
virmpole[5] -= vyz;
}
}
@ -781,7 +781,6 @@ void PairAmoeba::multipole_kspace()
for (k = 0; k < 20; k++)
fphi[i][k] *= felec;
}
//printf("fphi elec %e %e %e %e \n",fphi[0][0],fphi[0][1],fphi[0][2],fphi[0][3]);
// convert field from fractional to Cartesian
@ -811,7 +810,6 @@ void PairAmoeba::multipole_kspace()
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
@ -878,12 +876,12 @@ void PairAmoeba::multipole_kspace()
// increment total internal virial tensor components
virmpole[0] += vxx;
virmpole[1] += vyy;
virmpole[2] += vzz;
virmpole[3] += vxy;
virmpole[4] += vxz;
virmpole[5] += vyz;
virmpole[0] -= vxx;
virmpole[1] -= vyy;
virmpole[2] -= vzz;
virmpole[3] -= vxy;
virmpole[4] -= vxz;
virmpole[5] -= vyz;
// free local memory

View File

@ -129,12 +129,12 @@ void PairAmoeba::polar()
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
virpolar[0] += vxx;
virpolar[1] += vyy;
virpolar[2] += vzz;
virpolar[3] += vxy;
virpolar[4] += vxz;
virpolar[5] += vyz;
virpolar[0] -= vxx;
virpolar[1] -= vyy;
virpolar[2] -= vzz;
virpolar[3] -= vxy;
virpolar[4] -= vxz;
virpolar[5] -= vyz;
}
// clean up
@ -1155,12 +1155,12 @@ void PairAmoeba::polar_real()
vyz = 0.5 * (zr*frcy+yr*frcz);
vzz = zr * frcz;
virpolar[0] += vxx;
virpolar[1] += vyy;
virpolar[2] += vzz;
virpolar[3] += vxy;
virpolar[4] += vxz;
virpolar[5] += vyz;
virpolar[0] -= vxx;
virpolar[1] -= vyy;
virpolar[2] -= vzz;
virpolar[3] -= vxy;
virpolar[4] -= vxz;
virpolar[5] -= vyz;
// energy = e
// virial = 6-vec vir
@ -1224,12 +1224,12 @@ void PairAmoeba::polar_real()
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
virpolar[0] += vxx;
virpolar[1] += vyy;
virpolar[2] += vzz;
virpolar[3] += vxy;
virpolar[4] += vxz;
virpolar[5] += vyz;
virpolar[0] -= vxx;
virpolar[1] -= vyy;
virpolar[2] -= vzz;
virpolar[3] -= vxy;
virpolar[4] -= vxz;
virpolar[5] -= vyz;
}
}
@ -2125,12 +2125,12 @@ void PairAmoeba::polar_kspace()
qgrid[k3][k2][k1][1]*qgrip[k3][k2][k1][1];
eterm = 0.5 * felec * expterm * struc2;
vterm = (2.0/hsq) * (1.0-term) * eterm;
virpolar[0] += h1*h1*vterm - eterm;
virpolar[1] += h2*h2*vterm - eterm;
virpolar[2] += h3*h3*vterm - eterm;
virpolar[3] += h1*h2*vterm;
virpolar[4] += h1*h3*vterm;
virpolar[5] += h2*h3*vterm;
virpolar[0] -= h1*h1*vterm - eterm;
virpolar[1] -= h2*h2*vterm - eterm;
virpolar[2] -= h3*h3*vterm - eterm;
virpolar[3] -= h1*h2*vterm;
virpolar[4] -= h1*h3*vterm;
virpolar[5] -= h2*h3*vterm;
}
}
@ -2189,12 +2189,12 @@ void PairAmoeba::polar_kspace()
qgrid[k3][k2][k1][1]*qgrip[k3][k2][k1][1];
eterm = 0.5 * felec * expterm * struc2;
vterm = (2.0/hsq) * (1.0-term) * eterm;
virpolar[0] += h1*h1*vterm - eterm;
virpolar[1] += h2*h2*vterm - eterm;
virpolar[2] += h3*h3*vterm - eterm;
virpolar[3] += h1*h2*vterm;
virpolar[4] += h1*h3*vterm;
virpolar[5] += h2*h3*vterm;
virpolar[0] -= h1*h1*vterm - eterm;
virpolar[1] -= h2*h2*vterm - eterm;
virpolar[2] -= h3*h3*vterm - eterm;
virpolar[3] -= h1*h2*vterm;
virpolar[4] -= h1*h3*vterm;
virpolar[5] -= h2*h3*vterm;
}
}
}
@ -2203,12 +2203,12 @@ void PairAmoeba::polar_kspace()
// increment the total internal virial tensor components
virpolar[0] += vxx;
virpolar[1] += vyy;
virpolar[2] += vzz;
virpolar[3] += vxy;
virpolar[4] += vxz;
virpolar[5] += vyz;
virpolar[0] -= vxx;
virpolar[1] -= vyy;
virpolar[2] -= vzz;
virpolar[3] -= vxy;
virpolar[4] -= vxz;
virpolar[5] -= vyz;
// deallocation of local arrays, some from induce

View File

@ -361,12 +361,12 @@ void PairAmoeba::repulsion()
vyz = -0.5 * (zr*frcy+yr*frcz);
vzz = -zr * frcz;
virrepulse[0] += vxx;
virrepulse[1] += vyy;
virrepulse[2] += vzz;
virrepulse[3] += vxy;
virrepulse[4] += vxz;
virrepulse[5] += vyz;
virrepulse[0] -= vxx;
virrepulse[1] -= vyy;
virrepulse[2] -= vzz;
virrepulse[3] -= vxy;
virrepulse[4] -= vxz;
virrepulse[5] -= vyz;
// energy = e
// virial = 6-vec vir
@ -414,12 +414,12 @@ void PairAmoeba::repulsion()
vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
virrepulse[0] += vxx;
virrepulse[1] += vyy;
virrepulse[2] += vzz;
virrepulse[3] += vxy;
virrepulse[4] += vxz;
virrepulse[5] += vyz;
virrepulse[0] -= vxx;
virrepulse[1] -= vyy;
virrepulse[2] -= vzz;
virrepulse[3] -= vxy;
virrepulse[4] -= vxz;
virrepulse[5] -= vyz;
// virial = 6-vec vir
// NOTE: add tally function

View File

@ -416,7 +416,13 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag)
f[i4][2] -= f4[2];
}
if (evflag) ev_tally4(i1,i2,i3,14,nlocal,newton_bond,eangle,f1,f2,f3,f4);
if (evflag) {
f1[0] = -f1[0]; f1[1] = -f1[1]; f1[2] = -f1[2];
f2[0] = -f2[0]; f2[1] = -f2[1]; f2[2] = -f2[2];
f3[0] = -f3[0]; f3[1] = -f3[1]; f3[2] = -f3[2];
f4[0] = -f4[0]; f4[1] = -f4[1]; f4[2] = -f4[2];
ev_tally4(i1,i2,i3,14,nlocal,newton_bond,eangle,f1,f2,f3,f4);
}
}
/* ---------------------------------------------------------------------- */

View File

@ -555,12 +555,12 @@ void FixAmoebaPiTorsion::post_force(int vflag)
vxterm = dedxid + dedxia + dedxib;
vyterm = dedyid + dedyia + dedyib;
vzterm = dedzid + dedzia + dedzib;
v[0] = xdc*vxterm + xcp*dedxip - xqd*dedxiq;
v[1] = ydc*vyterm + ycp*dedyip - yqd*dedyiq;
v[2] = zdc*vzterm + zcp*dedzip - zqd*dedziq;
v[3] = ydc*vxterm + ycp*dedxip - yqd*dedxiq;
v[4] = zdc*vxterm + zcp*dedxip - zqd*dedxiq;
v[5] = zdc*vyterm + zcp*dedyip - zqd*dedyiq;
v[0] = -xdc*vxterm - xcp*dedxip + xqd*dedxiq;
v[1] = -ydc*vyterm - ycp*dedyip + yqd*dedyiq;
v[2] = -zdc*vzterm - zcp*dedzip + zqd*dedziq;
v[3] = -ydc*vxterm - ycp*dedxip + yqd*dedxiq;
v[4] = -zdc*vxterm - zcp*dedxip + zqd*dedxiq;
v[5] = -zdc*vyterm - zcp*dedyip + zqd*dedyiq;
ev_tally(nlist,list,6.0,e,v);
}

View File

@ -227,9 +227,13 @@ void ImproperAmoeba::compute(int eflag, int vflag)
f[ic][2] -= fc[2];
}
if (evflag)
if (evflag) {
fd[0] = -fd[0]; fd[1] = -fd[1]; fd[2] = -fd[2];
fa[0] = -fa[0]; fa[1] = -fa[1]; fa[2] = -fa[2];
fc[0] = -fc[0]; fc[1] = -fc[1]; fc[2] = -fc[2];
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);
}
}
}