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 *);