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

update angle.h and angle.cpp to compute per-atom centroid virial according to vflag=8

Browse Source
This commit is contained in:
Donatas Surblys
2019-10-04 18:05:25 +09:00
parent 0067a5b0fe
commit 239edc3442
2 changed files with 104 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

103
src/angle.cpp
View File

@ -34,9 +34,10 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
allocated = 0;
suffix_flag = Suffix::NONE;
maxeatom = maxvatom = 0;
maxeatom = maxvatom = maxcvatom = 0;
eatom = NULL;
vatom = NULL;
cvatom = NULL;
setflag = NULL;
execution_space = Host;
@ -54,6 +55,7 @@ Angle::~Angle()
memory->destroy(eatom);
memory->destroy(vatom);
memory->destroy(cvatom);
}
/* ----------------------------------------------------------------------
@ -85,9 +87,10 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
eflag_global = eflag % 2;
eflag_atom = eflag / 2;
vflag_either = vflag;
vflag_global = vflag % 4;
vflag_atom = vflag / 4;
vflag_atom = vflag & 4;
cvflag_atom = vflag & 8;
vflag_either = vflag_global || vflag_atom;
// reallocate per-atom arrays if necessary
@ -105,6 +108,13 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
memory->create(vatom,comm->nthreads*maxvatom,6,"angle:vatom");
}
}
if (cvflag_atom && atom->nmax > maxcvatom) {
maxcvatom = atom->nmax;
if (alloc) {
memory->destroy(cvatom);
memory->create(cvatom,comm->nthreads*maxcvatom,9,"angle:cvatom");
}
}
// zero accumulators
@ -127,6 +137,22 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
vatom[i][5] = 0.0;
}
}
if (cvflag_atom && alloc) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
for (i = 0; i < n; i++) {
cvatom[i][0] = 0.0;
cvatom[i][1] = 0.0;
cvatom[i][2] = 0.0;
cvatom[i][3] = 0.0;
cvatom[i][4] = 0.0;
cvatom[i][5] = 0.0;
cvatom[i][6] = 0.0;
cvatom[i][7] = 0.0;
cvatom[i][8] = 0.0;
cvatom[i][9] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
@ -230,6 +256,76 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
}
}
}
// per-atom centroid virial
if (cvflag_atom) {
// r0 = (r1+r2+r3)/3
// rij = ri-rj
// total virial = r10*f1 + r20*f2 + r30*f3
// del1: r12
// del2: r32
if (newton_bond || i < nlocal) {
double a1[3];
// a1 = r10 = (2*r12 - r32)/3
a1[0] = THIRD*(2*delx1-delx2);
a1[1] = THIRD*(2*dely1-dely2);
a1[2] = THIRD*(2*delz1-delz2);
cvatom[i][0] += a1[0]*f1[0];
cvatom[i][1] += a1[1]*f1[1];
cvatom[i][2] += a1[2]*f1[2];
cvatom[i][3] += a1[0]*f1[1];
cvatom[i][4] += a1[0]*f1[2];
cvatom[i][5] += a1[1]*f1[2];
cvatom[i][6] += a1[1]*f1[0];
cvatom[i][7] += a1[2]*f1[0];
cvatom[i][8] += a1[2]*f1[1];
}
if (newton_bond || j < nlocal) {
double a2[3];
double f2[3];
// a2 = r20 = ( -r12 - r32)/3
a2[0] = THIRD*(-delx1-delx2);
a2[1] = THIRD*(-dely1-dely2);
a2[2] = THIRD*(-delz1-delz2);
f2[0] = - f1[0] - f3[0];
f2[1] = - f1[1] - f3[1];
f2[2] = - f1[2] - f3[2];
cvatom[j][0] += a2[0]*f2[0];
cvatom[j][1] += a2[1]*f2[1];
cvatom[j][2] += a2[2]*f2[2];
cvatom[j][3] += a2[0]*f2[1];
cvatom[j][4] += a2[0]*f2[2];
cvatom[j][5] += a2[1]*f2[2];
cvatom[j][6] += a2[1]*f2[0];
cvatom[j][7] += a2[2]*f2[0];
cvatom[j][8] += a2[2]*f2[1];
}
if (newton_bond || k < nlocal) {
double a3[3];
// a3 = r30 = ( -r12 + 2*r32)/3
a3[0] = THIRD*(-delx1+2*delx2);
a3[1] = THIRD*(-dely1+2*dely2);
a3[2] = THIRD*(-delz1+2*delz2);
cvatom[k][0] += a3[0]*f3[0];
cvatom[k][1] += a3[1]*f3[1];
cvatom[k][2] += a3[2]*f3[2];
cvatom[k][3] += a3[0]*f3[1];
cvatom[k][4] += a3[0]*f3[2];
cvatom[k][5] += a3[1]*f3[2];
cvatom[k][6] += a3[1]*f3[0];
cvatom[k][7] += a3[2]*f3[0];
cvatom[k][8] += a3[2]*f3[1];
}
}
}
/* ---------------------------------------------------------------------- */
@ -238,5 +334,6 @@ double Angle::memory_usage()
{
double bytes = comm->nthreads*maxeatom * sizeof(double);
bytes += comm->nthreads*maxvatom*6 * sizeof(double);
bytes += comm->nthreads*maxcvatom*9 * sizeof(double);
return bytes;
}

7
src/angle.h
View File

@ -28,6 +28,7 @@ class Angle : protected Pointers {
double energy; // accumulated energies
double virial[6]; // accumulated virial
double *eatom,**vatom; // accumulated per-atom energy/virial
double **cvatom; // accumulated per-atom centroid virial
// KOKKOS host/device flag and data masks
@ -56,12 +57,12 @@ class Angle : protected Pointers {
int evflag;
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int maxeatom,maxvatom;
int vflag_either,vflag_global,vflag_atom,cvflag_atom;
int maxeatom,maxvatom,maxcvatom;
void ev_init(int eflag, int vflag, int alloc = 1) {
if (eflag||vflag) ev_setup(eflag, vflag, alloc);
else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = 0;
else evflag = eflag_either = eflag_global = eflag_atom = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0;
}
void ev_setup(int, int, int alloc = 1);
void ev_tally(int, int, int, int, int, double, double *, double *,