diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 99836483dd..c284d99011 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -52,10 +52,10 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr), vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), - b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), - b_max(nullptr), b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), - a_count(nullptr), a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), - a_ave_all(nullptr), a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr) + b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr), + b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr), + a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr), + a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr) { energy_global_flag = energy_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1; @@ -204,8 +204,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : int nb = atom->nbondtypes + 1; b_count = new int[nb]; b_count_all = new int[nb]; - b_atom = new int[nb]; - b_atom_all = new int[nb]; b_ave = new double[nb]; b_ave_all = new double[nb]; b_max = new double[nb]; @@ -298,8 +296,6 @@ FixShake::~FixShake() if (output_every) { delete[] b_count; delete[] b_count_all; - delete[] b_atom; - delete[] b_atom_all; delete[] b_ave; delete[] b_ave_all; delete[] b_max; @@ -2547,7 +2543,6 @@ void FixShake::bond_force(tagint id1, tagint id2, double length) void FixShake::stats() { - int i,j,m,n,iatom,jatom,katom; double delx,dely,delz; double r,r1,r2,r3,angle; @@ -2556,13 +2551,12 @@ void FixShake::stats() int nb = atom->nbondtypes + 1; int na = atom->nangletypes + 1; - for (i = 0; i < nb; i++) { + for (int i = 0; i < nb; i++) { b_count[i] = 0; b_ave[i] = b_max[i] = 0.0; b_min[i] = BIG; - b_atom[i] = -1; } - for (i = 0; i < na; i++) { + for (int i = 0; i < na; i++) { a_count[i] = 0; a_ave[i] = a_max[i] = 0.0; a_min[i] = BIG; @@ -2574,25 +2568,26 @@ void FixShake::stats() double **x = atom->x; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) { - if (shake_flag[i] == 0) continue; + for (int ii = 0; ii < nlist; ++ii) { + int i = list[ii]; + int n = shake_flag[i]; + if (n == 0) continue; // bond stats - n = shake_flag[i]; if (n == 1) n = 3; - iatom = atom->map(shake_atom[i][0]); - for (j = 1; j < n; j++) { - jatom = atom->map(shake_atom[i][j]); + int iatom = atom->map(shake_atom[i][0]); + for (int j = 1; j < n; j++) { + int jatom = atom->map(shake_atom[i][j]); + if (jatom >= nlocal) continue; delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; delz = x[iatom][2] - x[jatom][2]; domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); - m = shake_type[i][j-1]; + r = sqrt(delx*delx + dely*dely + delz*delz); + int m = shake_type[i][j-1]; b_count[m]++; - b_atom[m] = n; b_ave[m] += r; b_max[m] = MAX(b_max[m],r); b_min[m] = MIN(b_min[m],r); @@ -2601,9 +2596,13 @@ void FixShake::stats() // angle stats if (shake_flag[i] == 1) { - iatom = atom->map(shake_atom[i][0]); - jatom = atom->map(shake_atom[i][1]); - katom = atom->map(shake_atom[i][2]); + int iatom = atom->map(shake_atom[i][0]); + int jatom = atom->map(shake_atom[i][1]); + int katom = atom->map(shake_atom[i][2]); + int n = 0; + if (iatom < nlocal) ++n; + if (jatom < nlocal) ++n; + if (katom < nlocal) ++n; delx = x[iatom][0] - x[jatom][0]; dely = x[iatom][1] - x[jatom][1]; @@ -2625,9 +2624,9 @@ void FixShake::stats() angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2)); angle *= 180.0/MY_PI; - m = shake_type[i][2]; - a_count[m]++; - a_ave[m] += angle; + int m = shake_type[i][2]; + a_count[m] += n; + a_ave[m] += n*angle; a_max[m] = MAX(a_max[m],angle); a_min[m] = MIN(a_min[m],angle); } @@ -2636,7 +2635,6 @@ void FixShake::stats() // sum across all procs MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world); - MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world); MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world); MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world); @@ -2652,13 +2650,13 @@ void FixShake::stats() const int width = log10((MAX(MAX(1,nb),na)))+2; auto mesg = fmt::format("{} stats (type/ave/delta/count) on step {}\n", utils::uppercase(style), update->ntimestep); - for (i = 1; i < nb; i++) { + for (int i = 1; i < nb; i++) { const auto bcnt = b_count_all[i]; if (bcnt) mesg += fmt::format("Bond: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, - b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]); + b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt); } - for (i = 1; i < na; i++) { + for (int i = 1; i < na; i++) { const auto acnt = a_count_all[i]; if (acnt) mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index d6f9a8f5d6..820d68ebe7 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -113,8 +113,7 @@ class FixShake : public Fix { int nlist, maxlist; // size and max-size of list // stat quantities - int *b_count, *b_count_all, *b_atom, - *b_atom_all; // counts for each bond type, atoms in bond cluster + int *b_count, *b_count_all; // counts for each bond type, atoms in bond cluster double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays int *a_count, *a_count_all; // ditto for angle types