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

fix shake stats (again)

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-06-10 01:41:14 -04:00
parent 322bf1ef47
commit 1ee35bea61
2 changed files with 30 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
src/RIGID/fix_shake.cpp
View File

@ -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,

3
src/RIGID/fix_shake.h
View File

@ -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