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

Fix mistake in virial calculation for improper_fourier and improper_umbrella

Browse Source
This commit is contained in:
Steven Vandenbrande
2019-01-15 10:06:42 +01:00
parent a5ce656c32
commit 6702f65fbb
4 changed files with 91 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/USER-MISC/improper_fourier.cpp
View File

@ -206,17 +206,17 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
dahy = ary-c*hry;
dahz = arz-c*hrz;
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
f4[0] = dahx*rhr;
f4[1] = dahy*rhr;
f4[2] = dahz*rhr;
f4[0] = dahx*rhr*a;
f4[1] = dahy*rhr*a;
f4[2] = dahz*rhr*a;
f1[0] = -(f2[0] + f3[0] + f4[0]);
f1[1] = -(f2[1] + f3[1] + f4[1]);
@ -225,32 +225,32 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
// apply force to each of 4 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0]*a;
f[i1][1] += f1[1]*a;
f[i1][2] += f1[2]*a;
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] += f3[0]*a;
f[i2][1] += f3[1]*a;
f[i2][2] += f3[2]*a;
f[i2][0] += f3[0];
f[i2][1] += f3[1];
f[i2][2] += f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f2[0]*a;
f[i3][1] += f2[1]*a;
f[i3][2] += f2[2]*a;
f[i3][0] += f2[0];
f[i3][1] += f2[1];
f[i3][2] += f2[2];
}
if (newton_bond || i4 < nlocal) {
f[i4][0] += f4[0]*a;
f[i4][1] += f4[1]*a;
f[i4][2] += f4[2]*a;
f[i4][0] += f4[0];
f[i4][1] += f4[1];
f[i4][2] += f4[2];
}
if (evflag)
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
-vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z);
}
/* ---------------------------------------------------------------------- */