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

energy/virial flags with CENTROID enum

Browse Source
This commit is contained in:
Steve Plimpton
2020-11-18 12:17:14 -07:00
parent 35035189e6
commit 0876684780
74 changed files with 362 additions and 262 deletions
Download Patch File Download Diff File Expand all files Collapse all files

57
src/compute_centroid_stress_atom.cpp
View File

@ -124,10 +124,35 @@ void ComputeCentroidStressAtom::init()
else biasflag = NOBIAS;
} else biasflag = NOBIAS;
// check if pair styles support centroid atom stress
// check if force components and fixes support centroid atom stress
// all bond styles support it as CENTROID_SAME
if (pairflag && force->pair)
if (force->pair->centroidstressflag & CENTROID_NOTAVAIL)
if (force->pair->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "Pair style does not support compute centroid/stress/atom");
if (angleflag && force->angle)
if (force->angle->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "Angle style does not support compute centroid/stress/atom");
if (dihedralflag && force->dihedral)
if (force->dihedral->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "Dihedral style does not support compute centroid/stress/atom");
if (improperflag && force->improper)
if (force->improper->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "Improper style does not support compute centroid/stress/atom");
if (kspaceflag && force->kspace)
if (force->kspace->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "KSpace style does not support compute centroid/stress/atom");
if (fixflag) {
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (modify->fix[ifix]->virial_flag &&
modify->fix[ifix]->centroidstressflag == CENTROID_NOTAVAIL)
error->all(FLERR, "Fix style does not support compute centroid/stress/atom");
}
}
/* ---------------------------------------------------------------------- */
@ -173,12 +198,12 @@ void ComputeCentroidStressAtom::compute_peratom()
for (j = 0; j < 9; j++)
stress[i][j] = 0.0;
// add in per-atom contributions from each force
// per-atom virial and per-atom centroid virial are the same for two-body
// many-body pair styles not yet implemented
// add in per-atom contributions from all force components and fixes
// pair styles are either CENTROID_SAME or CENTROID_AVAIL or CENTROID_NOTAVAIL
if (pairflag && force->pair && force->pair->compute_flag) {
if (force->pair->centroidstressflag & CENTROID_AVAIL) {
if (force->pair->centroidstressflag == CENTROID_AVAIL) {
double **cvatom = force->pair->cvatom;
for (i = 0; i < npair; i++)
for (j = 0; j < 9; j++)
@ -195,6 +220,9 @@ void ComputeCentroidStressAtom::compute_peratom()
}
// per-atom virial and per-atom centroid virial are the same for bonds
// bond styles are all CENTROID_SAME
// all other styles are CENTROID_AVAIL or CENTROID_NOTAVAIL
if (bondflag && force->bond) {
double **vatom = force->bond->vatom;
for (i = 0; i < nbond; i++) {
@ -228,12 +256,9 @@ void ComputeCentroidStressAtom::compute_peratom()
if (kspaceflag && force->kspace && force->kspace->compute_flag) {
double **vatom = force->kspace->vatom;
for (i = 0; i < nkspace; i++) {
for (j = 0; j < 6; j++)
stress[i][j] += vatom[i][j];
for (j = 6; j < 9; j++)
stress[i][j] += vatom[i][j-3];
}
for (i = 0; i < nkspace; i++)
for (j = 0; j < 9; j++)
stress[i][j] += cvatom[i][j];
}
// add in per-atom contributions from relevant fixes
@ -241,7 +266,8 @@ void ComputeCentroidStressAtom::compute_peratom()
// possible during setup phase if fix has not initialized its vatom yet
// e.g. fix ave/spatial defined before fix shake,
// and fix ave/spatial uses a per-atom stress from this compute as input
// fix styles are CENTROID_SAME or CENTROID_NOTAVAIL
if (fixflag) {
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (modify->fix[ifix]->virial_flag) {
@ -258,7 +284,8 @@ void ComputeCentroidStressAtom::compute_peratom()
// communicate ghost virials between neighbor procs
if (force->newton || (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag))
if (force->newton ||
(force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag))
comm->reverse_comm_compute(this);
// zero virial of atoms not in group