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

various bugfixes to reference /tally computes

Browse Source 
- remove vector compute to avoid ambiguity when using atom style variables. fold reduction into compute_scalar.
- count pair contributions only if one of the two is one group and the other atom in the other group. groups can overlap (all/all counts all pairs)
- after using reverse communication to collect contributions from ghost atoms, we need to zero out the data from ghosts, or they'll be counted twice with compute reduce.

thanks to loris ercole for demonstrating the issues
This commit is contained in:
Axel Kohlmeyer
2015-06-10 13:50:47 -04:00
parent dee993230d
commit 1638e59f30
6 changed files with 103 additions and 112 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
src/compute_force_tally.cpp
View File

@ -36,19 +36,17 @@ ComputeForceTally::ComputeForceTally(LAMMPS *lmp, int narg, char **arg) :
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 3;
extscalar = 0;
extvector = 0;
size_vector = 3;
extscalar = 1;
peflag = 1; // we need Pair::ev_tally() to be run
nmax = -1;
fatom = NULL;
vector = new double[size_vector];
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
@ -107,19 +105,23 @@ void ComputeForceTally::pair_tally_callback(int i, int j, int nlocal, int newton
fatom[i][j] = 0.0;
}
for (int i=0; i < size_vector; ++i)
for (int i=0; i < size_peratom_cols; ++i)
vector[i] = ftotal[i] = 0.0;
}
if ((mask[i] & groupbit) && (newton || i < nlocal)) {
ftotal[0] += fpair*dx; fatom[i][0] += fpair*dx;
ftotal[1] += fpair*dy; fatom[i][1] += fpair*dy;
ftotal[2] += fpair*dz; fatom[i][2] += fpair*dz;
}
if ((mask[j] & groupbit) && (newton || j < nlocal)) {
ftotal[0] -= fpair*dx; fatom[i][0] -= fpair*dx;
ftotal[1] -= fpair*dy; fatom[i][1] -= fpair*dy;
ftotal[2] -= fpair*dz; fatom[i][2] -= fpair*dz;
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ) {
if (newton || i < nlocal) {
ftotal[0] += fpair*dx; fatom[i][0] += fpair*dx;
ftotal[1] += fpair*dy; fatom[i][1] += fpair*dy;
ftotal[2] += fpair*dz; fatom[i][2] += fpair*dz;
}
if (newton || j < nlocal) {
ftotal[0] -= fpair*dx; fatom[i][0] -= fpair*dx;
ftotal[1] -= fpair*dy; fatom[i][1] -= fpair*dy;
ftotal[2] -= fpair*dz; fatom[i][2] -= fpair*dz;
}
}
}
@ -162,7 +164,9 @@ double ComputeForceTally::compute_scalar()
if (update->eflag_global != invoked_scalar)
error->all(FLERR,"Energy was not tallied on needed timestep");
compute_vector();
// sum accumulated forces across procs
MPI_Allreduce(ftotal,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
scalar = sqrt(vector[0]*vector[0]+vector[1]*vector[1]+vector[2]*vector[2]);
return scalar;
@ -170,27 +174,22 @@ double ComputeForceTally::compute_scalar()
/* ---------------------------------------------------------------------- */
void ComputeForceTally::compute_vector()
{
invoked_vector = update->ntimestep;
if (update->eflag_global != invoked_vector)
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated force across procs
MPI_Allreduce(ftotal,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
}
/* ---------------------------------------------------------------------- */
void ComputeForceTally::compute_peratom()
{
invoked_peratom = update->ntimestep;
if (update->eflag_global != invoked_peratom)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (force->newton_pair)
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
for (int j = 0; j < size_peratom_cols; ++j)
fatom[i][j] = 0.0;
}
}
/* ----------------------------------------------------------------------

1
src/compute_force_tally.h
View File

@ -33,7 +33,6 @@ class ComputeForceTally : public Compute {
void init();
double compute_scalar();
void compute_vector();
void compute_peratom();
int pack_reverse_comm(int, int, double *);

57
src/compute_pe_tally.cpp
View File

@ -36,19 +36,17 @@ ComputePETally::ComputePETally(LAMMPS *lmp, int narg, char **arg) :
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 2;
extscalar = 0;
extvector = 0;
size_vector = 2;
extscalar = 1;
peflag = 1; // we need Pair::ev_tally() to be run
nmax = -1;
eatom = NULL;
vector = new double[size_vector];
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
@ -108,14 +106,18 @@ void ComputePETally::pair_tally_callback(int i, int j, int nlocal, int newton,
vector[0] = etotal[0] = vector[1] = etotal[1] = 0.0;
}
evdwl *= 0.5; ecoul *= 0.5;
if ((mask[i] & groupbit) && (newton || i < nlocal)) {
etotal[0] += evdwl; eatom[i][0] += evdwl;
etotal[1] += ecoul; eatom[i][1] += ecoul;
}
if ((mask[j] & groupbit) && (newton || j < nlocal)) {
etotal[0] += evdwl; eatom[j][0] += evdwl;
etotal[1] += ecoul; eatom[j][1] += ecoul;
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ){
evdwl *= 0.5; ecoul *= 0.5;
if (newton || i < nlocal) {
etotal[0] += evdwl; eatom[i][0] += evdwl;
etotal[1] += ecoul; eatom[i][1] += ecoul;
}
if (newton || j < nlocal) {
etotal[0] += evdwl; eatom[j][0] += evdwl;
etotal[1] += ecoul; eatom[j][1] += ecoul;
}
}
}
@ -156,22 +158,12 @@ double ComputePETally::compute_scalar()
if (update->eflag_global != invoked_scalar)
error->all(FLERR,"Energy was not tallied on needed timestep");
compute_vector();
scalar = vector[0]+vector[1];
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputePETally::compute_vector()
{
invoked_vector = update->ntimestep;
if (update->eflag_global != invoked_vector)
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated energies across procs
MPI_Allreduce(etotal,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(etotal,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
scalar = vector[0]+vector[1];
return scalar;
}
/* ---------------------------------------------------------------------- */
@ -182,10 +174,15 @@ void ComputePETally::compute_peratom()
if (update->eflag_global != invoked_peratom)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (force->newton_pair)
comm->reverse_comm_compute(this);
}
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
eatom[i][0] = eatom[i][1] = 0.0;
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array

1
src/compute_pe_tally.h
View File

@ -33,7 +33,6 @@ class ComputePETally : public Compute {
void init();
double compute_scalar();
void compute_vector();
void compute_peratom();
int pack_reverse_comm(int, int, double *);

96
src/compute_stress_tally.cpp
View File

@ -36,19 +36,17 @@ ComputeStressTally::ComputeStressTally(LAMMPS *lmp, int narg, char **arg) :
groupbit2 = group->bitmask[igroup2];
scalar_flag = 1;
vector_flag = 1;
vector_flag = 0;
peratom_flag = 1;
timeflag = 1;
comm_reverse = size_peratom_cols = 6;
extscalar = 0;
extvector = 0;
size_vector = 6;
peflag = 1; // we need Pair::ev_tally() to be run
nmax = -1;
stress = NULL;
vector = new double[size_vector];
vector = new double[size_peratom_cols];
}
/* ---------------------------------------------------------------------- */
@ -107,34 +105,37 @@ void ComputeStressTally::pair_tally_callback(int i, int j, int nlocal, int newto
stress[i][j] = 0.0;
}
for (int i=0; i < size_vector; ++i)
for (int i=0; i < size_peratom_cols; ++i)
vector[i] = virial[i] = 0.0;
}
fpair *= 0.5;
const double v0 = dx*dx*fpair;
const double v1 = dy*dy*fpair;
const double v2 = dz*dz*fpair;
const double v3 = dx*dy*fpair;
const double v4 = dx*dz*fpair;
const double v5 = dy*dz*fpair;
if ( ((mask[i] & groupbit) && (mask[j] & groupbit2))
|| ((mask[i] & groupbit2) && (mask[j] & groupbit)) ) {
if ((mask[i] & groupbit) && (newton || i < nlocal)) {
virial[0] += v0; stress[i][0] += v0;
virial[1] += v1; stress[i][1] += v1;
virial[2] += v2; stress[i][2] += v2;
virial[3] += v3; stress[i][3] += v3;
virial[4] += v4; stress[i][4] += v4;
virial[5] += v5; stress[i][5] += v5;
}
if ((mask[j] & groupbit) && (newton || j < nlocal)) {
virial[0] += v0; stress[j][0] += v0;
virial[1] += v1; stress[j][1] += v1;
virial[2] += v2; stress[j][2] += v2;
virial[3] += v3; stress[j][3] += v3;
virial[4] += v4; stress[j][4] += v4;
virial[5] += v5; stress[j][5] += v5;
fpair *= 0.5;
const double v0 = dx*dx*fpair;
const double v1 = dy*dy*fpair;
const double v2 = dz*dz*fpair;
const double v3 = dx*dy*fpair;
const double v4 = dx*dz*fpair;
const double v5 = dy*dz*fpair;
if (newton || i < nlocal) {
virial[0] += v0; stress[i][0] += v0;
virial[1] += v1; stress[i][1] += v1;
virial[2] += v2; stress[i][2] += v2;
virial[3] += v3; stress[i][3] += v3;
virial[4] += v4; stress[i][4] += v4;
virial[5] += v5; stress[i][5] += v5;
}
if (newton || j < nlocal) {
virial[0] += v0; stress[j][0] += v0;
virial[1] += v1; stress[j][1] += v1;
virial[2] += v2; stress[j][2] += v2;
virial[3] += v3; stress[j][3] += v3;
virial[4] += v4; stress[j][4] += v4;
virial[5] += v5; stress[j][5] += v5;
}
}
}
@ -183,39 +184,37 @@ double ComputeStressTally::compute_scalar()
if (update->eflag_global != invoked_scalar)
error->all(FLERR,"Energy was not tallied on needed timestep");
compute_vector();
scalar = (vector[0]+vector[1]+vector[2])/3.0;
// sum accumulated forces across procs
MPI_Allreduce(virial,vector,size_peratom_cols,MPI_DOUBLE,MPI_SUM,world);
if (domain->dimension == 3)
scalar = (vector[0]+vector[1]+vector[2])/3.0;
else
scalar = (vector[0]+vector[1])/2.0;
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::compute_vector()
{
invoked_vector = update->ntimestep;
if (update->eflag_global != invoked_vector)
error->all(FLERR,"Energy was not tallied on needed timestep");
// sum accumulated virial across procs
MPI_Allreduce(virial,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
const double nktv2p = -force->nktv2p;
for (int i=0; i < size_vector; ++i)
vector[i] *= nktv2p;
}
/* ---------------------------------------------------------------------- */
void ComputeStressTally::compute_peratom()
{
invoked_peratom = update->ntimestep;
if (update->eflag_global != invoked_peratom)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (force->newton_pair)
// collect contributions from ghost atoms
if (force->newton_pair) {
comm->reverse_comm_compute(this);
const int nall = atom->nlocal + atom->nghost;
for (int i = atom->nlocal; i < nall; ++i)
for (int j = 0; j < size_peratom_cols; ++j)
stress[i][j] = 0.0;
}
// convert to stress*volume units = -pressure*volume
const double nktv2p = -force->nktv2p;
@ -229,7 +228,6 @@ void ComputeStressTally::compute_peratom()
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */

1
src/compute_stress_tally.h
View File

@ -33,7 +33,6 @@ class ComputeStressTally : public Compute {
void init();
double compute_scalar();
void compute_vector();
void compute_peratom();
int pack_reverse_comm(int, int, double *);