From 1638e59f30dfd37dcd5a17d6ce153b973d2b4d61 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 10 Jun 2015 13:50:47 -0400 Subject: [PATCH] various bugfixes to reference /tally computes - 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 --- src/compute_force_tally.cpp | 59 +++++++++++----------- src/compute_force_tally.h | 1 - src/compute_pe_tally.cpp | 57 ++++++++++----------- src/compute_pe_tally.h | 1 - src/compute_stress_tally.cpp | 96 ++++++++++++++++++------------------ src/compute_stress_tally.h | 1 - 6 files changed, 103 insertions(+), 112 deletions(-) diff --git a/src/compute_force_tally.cpp b/src/compute_force_tally.cpp index 47e9870bb9..31c7fd62d7 100644 --- a/src/compute_force_tally.cpp +++ b/src/compute_force_tally.cpp @@ -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; + } } /* ---------------------------------------------------------------------- diff --git a/src/compute_force_tally.h b/src/compute_force_tally.h index b1a1808055..fce6b0fdd5 100644 --- a/src/compute_force_tally.h +++ b/src/compute_force_tally.h @@ -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 *); diff --git a/src/compute_pe_tally.cpp b/src/compute_pe_tally.cpp index edb69219c4..3af6ba323d 100644 --- a/src/compute_pe_tally.cpp +++ b/src/compute_pe_tally.cpp @@ -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 diff --git a/src/compute_pe_tally.h b/src/compute_pe_tally.h index 920ce0e5c4..d502570dbc 100644 --- a/src/compute_pe_tally.h +++ b/src/compute_pe_tally.h @@ -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 *); diff --git a/src/compute_stress_tally.cpp b/src/compute_stress_tally.cpp index b8182f19e8..5658676861 100644 --- a/src/compute_stress_tally.cpp +++ b/src/compute_stress_tally.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/compute_stress_tally.h b/src/compute_stress_tally.h index 23dd2f816c..4ce2466bd3 100644 --- a/src/compute_stress_tally.h +++ b/src/compute_stress_tally.h @@ -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 *);