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

Make compute stress/mop compatible with triclinic boxes

Browse Source
This commit is contained in:
Evangelos Voyiatzis
2024-05-09 15:12:28 +02:00
committed by GitHub
parent 46a441451d
commit 7b007d82a0
1 changed files with 10 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/EXTRA-COMPUTE/compute_stress_mop.cpp
View File

@ -89,12 +89,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
error->all(FLERR, "Plane for compute stress/mop is out of bounds");
}
if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) {
pos1 = pos + domain->prd[dir];
} else {
pos1 = pos - domain->prd[dir];
}
// parse values until one isn't recognized
which = new int[3 * (narg - 5)];
@ -151,10 +145,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(
if (domain->dimension == 2 && dir == Z)
error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system");
// orthogonal simulation box
if (domain->triclinic != 0)
error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box");
// Initialize some variables
values_local = values_global = vector = nullptr;
@ -397,6 +387,9 @@ void ComputeStressMop::compute_pairs()
jlist = firstneigh[i];
jnum = numneigh[i];
xi[dir] -= pos;
domain->minimum_image(xi[0], xi[1], xi[2]);
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
@ -409,6 +402,10 @@ void ComputeStressMop::compute_pairs()
xj[0] = atom->x[j][0];
xj[1] = atom->x[j][1];
xj[2] = atom->x[j][2];
xj[dir] -= pos;
domain->minimum_image(xj[0], xj[1], xj[2]);
delx = xi[0] - xj[0];
dely = xi[1] - xj[1];
delz = xi[2] - xj[2];
@ -419,20 +416,19 @@ void ComputeStressMop::compute_pairs()
if (newton_pair || j < nlocal) {
//check if ij pair is across plane, add contribution to pressure
if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) {
if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) {
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p;
values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p;
values_local[m + 2] += fpair * (xi[2] - xj[2]) / area * nktv2p;
} else if (((xi[dir] < pos) && (xj[dir] > pos)) ||
((xi[dir] < pos1) && (xj[dir] > pos1))) {
} else if ((xi[dir] < 0.0) && (xj[dir] > 0.0)) {
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
values_local[m] -= fpair * (xi[0] - xj[0]) / area * nktv2p;
values_local[m + 1] -= fpair * (xi[1] - xj[1]) / area * nktv2p;
values_local[m + 2] -= fpair * (xi[2] - xj[2]) / area * nktv2p;
}
} else {
if (((xi[dir] > pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) {
if ((xi[dir] > 0.0) && (xj[dir] < 0.0)) {
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
values_local[m] += fpair * (xi[0] - xj[0]) / area * nktv2p;
values_local[m + 1] += fpair * (xi[1] - xj[1]) / area * nktv2p;