From 092720b69a1a453ad5b4ffd770596479e5915e69 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 1 Mar 2022 06:39:19 -0500 Subject: [PATCH 01/20] correctly compute bond count for shake statistics --- src/RIGID/fix_shake.cpp | 72 ++++++++++++++++++++--------------------- src/RIGID/fix_shake.h | 8 ++--- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 49116cb8cc..735a36ce55 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -51,12 +51,11 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : type_flag(nullptr), mass_list(nullptr), bond_distance(nullptr), angle_distance(nullptr), 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_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) + shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr), + b_count(nullptr), b_count_all(nullptr), b_atom(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) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); @@ -172,8 +171,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : if (imol == -1) error->all(FLERR,"Molecule template ID for fix shake does not exist"); if (atom->molecules[imol]->nset > 1 && comm->me == 0) - error->warning(FLERR,"Molecule template for " - "fix shake has multiple molecules"); + error->warning(FLERR,"Molecule template for fix shake has multiple molecules"); onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; iarg += 2; @@ -199,6 +197,7 @@ 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_ave = new double[nb]; b_ave_all = new double[nb]; b_max = new double[nb]; @@ -228,8 +227,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : find_clusters(); if (comm->me == 0) - utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n", - platform::walltime()-time1); + utils::logmesg(lmp," find clusters CPU = {:.3f} seconds\n",platform::walltime()-time1); // initialize list of SHAKE clusters to constrain @@ -281,32 +279,33 @@ FixShake::~FixShake() memory->destroy(vtmp); - delete [] bond_flag; - delete [] angle_flag; - delete [] type_flag; - delete [] mass_list; + delete[] bond_flag; + delete[] angle_flag; + delete[] type_flag; + delete[] mass_list; - delete [] bond_distance; - delete [] angle_distance; + delete[] bond_distance; + delete[] angle_distance; if (output_every) { - delete [] b_count; - delete [] b_count_all; - delete [] b_ave; - delete [] b_ave_all; - delete [] b_max; - delete [] b_max_all; - delete [] b_min; - delete [] b_min_all; + delete[] b_count; + delete[] b_count_all; + delete[] b_atom; + delete[] b_ave; + delete[] b_ave_all; + delete[] b_max; + delete[] b_max_all; + delete[] b_min; + delete[] b_min_all; - delete [] a_count; - delete [] a_count_all; - delete [] a_ave; - delete [] a_ave_all; - delete [] a_max; - delete [] a_max_all; - delete [] a_min; - delete [] a_min_all; + delete[] a_count; + delete[] a_count_all; + delete[] a_ave; + delete[] a_ave_all; + delete[] a_max; + delete[] a_max_all; + delete[] a_min; + delete[] a_min_all; } memory->destroy(list); @@ -2504,6 +2503,7 @@ void FixShake::stats() 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); @@ -2563,16 +2563,16 @@ void FixShake::stats() auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n", update->ntimestep); for (i = 1; i < nb; i++) { - const auto bcnt = b_count_all[i]/2; + 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/2.0,b_max_all[i]-b_min_all[i],bcnt); + b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom[i]); } for (i = 1; i < na; i++) { - const auto acnt = a_count_all[i]/3; + const auto acnt = a_count_all[i]; if (acnt) mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width, - a_ave_all[i]/acnt/3.0,a_max_all[i]-a_min_all[i],acnt); + a_ave_all[i]/acnt,a_max_all[i]-a_min_all[i],acnt/3); } utils::logmesg(lmp,mesg); } diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 5a956bd798..8a076afaba 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -108,10 +108,10 @@ class FixShake : public Fix { int nlist, maxlist; // size and max-size of list // stat quantities - int *b_count, *b_count_all; // counts for each bond type - 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 + int *b_count, *b_count_all, *b_atom; // counts for each bond type + 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 double *a_ave, *a_max, *a_min; double *a_ave_all, *a_max_all, *a_min_all; From 2a0e66164e4ec09491be35e92f8d441159e447bf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 1 Mar 2022 21:02:06 -0500 Subject: [PATCH 02/20] enable and apply clang-format --- src/CLASS2/pair_lj_class2.cpp | 391 +++++++------- src/CLASS2/pair_lj_class2_coul_cut.cpp | 297 ++++++----- src/CLASS2/pair_lj_class2_coul_long.cpp | 647 ++++++++++++------------ 3 files changed, 658 insertions(+), 677 deletions(-) diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp index 7762e261cd..066869e851 100644 --- a/src/CLASS2/pair_lj_class2.cpp +++ b/src/CLASS2/pair_lj_class2.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -65,13 +64,13 @@ PairLJClass2::~PairLJClass2() void PairLJClass2::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -104,34 +103,32 @@ void PairLJClass2::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -144,10 +141,10 @@ void PairLJClass2::compute(int eflag, int vflag) void PairLJClass2::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -165,8 +162,8 @@ void PairLJClass2::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -187,28 +184,28 @@ void PairLJClass2::compute_inner() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -219,10 +216,10 @@ void PairLJClass2::compute_inner() void PairLJClass2::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -243,10 +240,10 @@ void PairLJClass2::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -267,32 +264,32 @@ void PairLJClass2::compute_middle() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; jtype = type[j]; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -303,13 +300,13 @@ void PairLJClass2::compute_middle() void PairLJClass2::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcelj, factor_lj, rsw; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -327,8 +324,8 @@ void PairLJClass2::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -349,55 +346,54 @@ void PairLJClass2::compute_outer(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { if (rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } if (eflag) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; } if (vflag) { if (rsq <= cut_in_off_sq) { - r2inv = 1.0/rsq; - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fpair = factor_lj*forcelj*r2inv; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fpair = factor_lj * forcelj * r2inv; } else if (rsq < cut_in_on_sq) - fpair = factor_lj*forcelj*r2inv; + fpair = factor_lj * forcelj * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); } } } @@ -409,23 +405,21 @@ void PairLJClass2::compute_outer(int eflag, int vflag) void PairLJClass2::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(cut, np1, np1, "pair:cut"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -434,14 +428,14 @@ void PairLJClass2::allocate() void PairLJClass2::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; @@ -454,22 +448,22 @@ void PairLJClass2::settings(int narg, char **arg) void PairLJClass2::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); double cut_one = cut_global; - if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp); + if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut[i][j] = cut_one; @@ -478,7 +472,7 @@ void PairLJClass2::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -492,12 +486,12 @@ void PairLJClass2::init_style() int irequest; int respa = 0; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; } - irequest = neighbor->request(this,instance_me); + irequest = neighbor->request(this, instance_me); if (respa >= 1) { neighbor->requests[irequest]->respaouter = 1; @@ -507,10 +501,11 @@ void PairLJClass2::init_style() // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; } /* ---------------------------------------------------------------------- @@ -523,23 +518,22 @@ double PairLJClass2::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut[i][j] > 0.0)) { double ratio = sigma[i][j] / cut[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; @@ -550,7 +544,7 @@ double PairLJClass2::init_one(int i, int j) // check interior rRESPA cutoff if (cut_respa && cut[i][j] < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -559,22 +553,21 @@ double PairLJClass2::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut[i][j] * cut[i][j] * cut[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut[i][j]; @@ -588,14 +581,14 @@ void PairLJClass2::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut[i][j], sizeof(double), 1, fp); } } } @@ -609,21 +602,21 @@ void PairLJClass2::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -634,10 +627,10 @@ void PairLJClass2::read_restart(FILE *fp) void PairLJClass2::write_restart_settings(FILE *fp) { - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -648,26 +641,24 @@ void PairLJClass2::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } - /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairLJClass2::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -678,27 +669,25 @@ void PairLJClass2::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut[i][j]); } /* ---------------------------------------------------------------------- */ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, - double /*factor_coul*/, double factor_lj, - double &fforce) + double /*factor_coul*/, double factor_lj, double &fforce) { - double r2inv,rinv,r3inv,r6inv,forcelj,philj; + double r2inv, rinv, r3inv, r6inv, forcelj, philj; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - fforce = factor_lj*forcelj*r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + fforce = factor_lj * forcelj * r2inv; - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - return factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + return factor_lj * philj; } /* ---------------------------------------------------------------------- */ @@ -706,7 +695,7 @@ double PairLJClass2::single(int /*i*/, int /*j*/, int itype, int jtype, double r void *PairLJClass2::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index c4ac5ae7c8..844c979994 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,17 +13,17 @@ #include "pair_lj_class2_coul_cut.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -65,14 +64,14 @@ PairLJClass2CoulCut::~PairLJClass2CoulCut() void PairLJClass2CoulCut::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj; - double factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj; + double factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -110,47 +109,49 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - else forcecoul = 0.0; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + fpair = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq[itype][jtype]) - ecoul = factor_coul * qqrd2e * qtmp*q[j]*sqrt(r2inv); - else ecoul = 0.0; + ecoul = factor_coul * qqrd2e * qtmp * q[j] * sqrt(r2inv); + else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -165,26 +166,25 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) void PairLJClass2CoulCut::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(cut_coul,n+1,n+1,"pair:cut_coul"); - memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(cut_coul, np1, np1, "pair:cut_coul"); + memory->create(cut_coulsq, np1, np1, "pair:cut_coulsq"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -193,16 +193,18 @@ void PairLJClass2CoulCut::allocate() void PairLJClass2CoulCut::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul_global = cut_lj_global; - else cut_coul_global = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul_global = cut_lj_global; + else + cut_coul_global = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) { @@ -218,26 +220,25 @@ void PairLJClass2CoulCut::settings(int narg, char **arg) void PairLJClass2CoulCut::coeff(int narg, char **arg) { - if (narg < 4 || narg > 6) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 6) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; - if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp); - if (narg == 6) cut_coul_one = utils::numeric(FLERR,arg[5],false,lmp); + if (narg >= 5) cut_coul_one = cut_lj_one = utils::numeric(FLERR, arg[4], false, lmp); + if (narg == 6) cut_coul_one = utils::numeric(FLERR, arg[5], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; @@ -247,7 +248,7 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -256,10 +257,9 @@ void PairLJClass2CoulCut::coeff(int narg, char **arg) void PairLJClass2CoulCut::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style lj/class2/coul/cut requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/cut requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); } /* ---------------------------------------------------------------------- @@ -272,28 +272,27 @@ double PairLJClass2CoulCut::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); - cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]); + cut_coul[i][j] = mix_distance(cut_coul[i][i], cut_coul[j][j]); } - double cut = MAX(cut_lj[i][j],cut_coul[i][j]); + double cut = MAX(cut_lj[i][j], cut_coul[i][j]); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; cut_coulsq[j][i] = cut_coulsq[i][j]; @@ -310,22 +309,21 @@ double PairLJClass2CoulCut::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut; @@ -339,15 +337,15 @@ void PairLJClass2CoulCut::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); - fwrite(&cut_coul[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); + fwrite(&cut_coul[i][j], sizeof(double), 1, fp); } } } @@ -361,23 +359,23 @@ void PairLJClass2CoulCut::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -388,11 +386,11 @@ void PairLJClass2CoulCut::read_restart(FILE *fp) void PairLJClass2CoulCut::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -402,17 +400,17 @@ void PairLJClass2CoulCut::write_restart_settings(FILE *fp) void PairLJClass2CoulCut::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- @@ -421,8 +419,7 @@ void PairLJClass2CoulCut::read_restart_settings(FILE *fp) void PairLJClass2CoulCut::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -433,39 +430,38 @@ void PairLJClass2CoulCut::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { - double r2inv,rinv,r3inv,r6inv,forcecoul,forcelj,phicoul,philj; + double r2inv, rinv, r3inv, r6inv, forcecoul, forcelj, phicoul, philj; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - else forcecoul = 0.0; + forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; - fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; + fforce = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; if (rsq < cut_coulsq[itype][jtype]) { - phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - eng += factor_coul*phicoul; + phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + eng += factor_coul * phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + eng += factor_lj * philj; } return eng; @@ -476,9 +472,8 @@ double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, void *PairLJClass2CoulCut::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } - diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index 63265509c9..79773d6408 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,32 +13,32 @@ #include "pair_lj_class2_coul_long.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "kspace.h" -#include "update.h" -#include "respa.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "respa.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace MathConst; -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 +static constexpr double EWALD_F = 1.12837917; +static constexpr double EWALD_P = 0.3275911; +static constexpr double A1 = 0.254829592; +static constexpr double A2 = -0.284496736; +static constexpr double A3 = 1.421413741; +static constexpr double A4 = -1.453152027; +static constexpr double A5 = 1.061405429; /* ---------------------------------------------------------------------- */ @@ -79,16 +78,16 @@ PairLJClass2CoulLong::~PairLJClass2CoulLong() void PairLJClass2CoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itable,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj; - double grij,expm2,prefactor,t,erfc; - double factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itable, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double rsq, r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj; + double grij, expm2, prefactor, t, erfc; + double factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -126,75 +125,77 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - ecoul = prefactor*erfc; + ecoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; } - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -206,11 +207,11 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) void PairLJClass2CoulLong::compute_inner() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -231,8 +232,8 @@ void PairLJClass2CoulLong::compute_inner() double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -255,34 +256,35 @@ void PairLJClass2CoulLong::compute_inner() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -293,11 +295,11 @@ void PairLJClass2CoulLong::compute_inner() void PairLJClass2CoulLong::compute_middle() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair; - double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, fpair; + double rsq, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; double **f = atom->f; @@ -321,10 +323,10 @@ void PairLJClass2CoulLong::compute_middle() double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; + double cut_out_on_sq = cut_out_on * cut_out_on; + double cut_out_off_sq = cut_out_off * cut_out_off; // loop over neighbors of my atoms @@ -347,38 +349,39 @@ void PairLJClass2CoulLong::compute_middle() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { - r2inv = 1.0/rsq; - forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul; + r2inv = 1.0 / rsq; + forcecoul = qqrd2e * qtmp * q[j] * sqrt(r2inv); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * forcecoul; jtype = type[j]; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + fpair *= rsw * rsw * (3.0 - 2.0 * rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); + rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff; + fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0); } - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } } } @@ -389,17 +392,17 @@ void PairLJClass2CoulLong::compute_middle() void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double fraction,table; - double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; - double grij,expm2,prefactor,t,erfc; + int i, j, ii, jj, inum, jnum, itype, jtype, itable; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double fraction, table; + double r, rinv, r2inv, r3inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj; + double grij, expm2, prefactor, t, erfc; double rsw; - int *ilist,*jlist,*numneigh,**firstneigh; + int *ilist, *jlist, *numneigh, **firstneigh; double rsq; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -420,8 +423,8 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; + double cut_in_off_sq = cut_in_off * cut_in_off; + double cut_in_on_sq = cut_in_on * cut_in_on; // loop over neighbors of my atoms @@ -444,32 +447,30 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0); + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = qqrd2e * qtmp * q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2 - 1.0); if (rsq > cut_in_off_sq) { if (rsq < cut_in_on_sq) { - rsw = (r - cut_in_off)/cut_in_diff; - forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw); + rsw = (r - cut_in_off) / cut_in_diff; + forcecoul += prefactor * rsw * rsw * (3.0 - 2.0 * rsw); if (factor_coul < 1.0) - forcecoul -= - (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw); + forcecoul -= (1.0 - factor_coul) * prefactor * rsw * rsw * (3.0 - 2.0 * rsw); } else { forcecoul += prefactor; - if (factor_coul < 1.0) - forcecoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } } } else { @@ -478,96 +479,99 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = qtmp*q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - forcelj *= rsw*rsw*(3.0 - 2.0*rsw); + rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff; + forcelj *= rsw * rsw * (3.0 - 2.0 * rsw); } - } else forcelj = 0.0; + } else + forcelj = 0.0; fpair = (forcecoul + forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - ecoul = prefactor*erfc; - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + ecoul = prefactor * erfc; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; } else { - table = etable[itable] + fraction*detable[itable]; - ecoul = qtmp*q[j] * table; + table = etable[itable] + fraction * detable[itable]; + ecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - ecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + ecoul -= (1.0 - factor_coul) * prefactor; } } - } else ecoul = 0.0; + } else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + evdwl = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } if (vflag) { if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { - table = vtable[itable] + fraction*dvtable[itable]; - forcecoul = qtmp*q[j] * table; + table = vtable[itable] + fraction * dvtable[itable]; + forcecoul = qtmp * q[j] * table; if (factor_coul < 1.0) { - table = ptable[itable] + fraction*dptable[itable]; - prefactor = qtmp*q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ptable[itable] + fraction * dptable[itable]; + prefactor = qtmp * q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq <= cut_in_off_sq) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } else if (rsq <= cut_in_on_sq) { - rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); } - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + fpair = (forcecoul + factor_lj * forcelj) * r2inv; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -580,24 +584,22 @@ void PairLJClass2CoulLong::compute_outer(int eflag, int vflag) void PairLJClass2CoulLong::allocate() { allocated = 1; - int n = atom->ntypes; + const int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + memory->create(setflag, np1, np1, "pair:setflag"); + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(epsilon, np1, np1, "pair:epsilon"); + memory->create(sigma, np1, np1, "pair:sigma"); + memory->create(lj1, np1, np1, "pair:lj1"); + memory->create(lj2, np1, np1, "pair:lj2"); + memory->create(lj3, np1, np1, "pair:lj3"); + memory->create(lj4, np1, np1, "pair:lj4"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -606,16 +608,18 @@ void PairLJClass2CoulLong::allocate() void PairLJClass2CoulLong::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul = cut_lj_global; - else cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul = cut_lj_global; + else + cut_coul = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; @@ -628,24 +632,23 @@ void PairLJClass2CoulLong::settings(int narg, char **arg) void PairLJClass2CoulLong::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); - double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); + double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[3], false, lmp); - double cut_lj_one = cut_lj_global; - if (narg == 5) cut_lj_one = utils::numeric(FLERR,arg[4],false,lmp); + double cut_lj_one = cut_lj_global; + if (narg == 5) cut_lj_one = utils::numeric(FLERR, arg[4], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; cut_lj[i][j] = cut_lj_one; @@ -654,7 +657,7 @@ void PairLJClass2CoulLong::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -663,25 +666,23 @@ void PairLJClass2CoulLong::coeff(int narg, char **arg) void PairLJClass2CoulLong::init_style() { - if (!atom->q_flag) - error->all(FLERR, - "Pair style lj/class2/coul/long requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style lj/class2/coul/long requires atom attribute q"); // request regular or rRESPA neighbor list int irequest; int respa = 0; - if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { + if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) { if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; } - irequest = neighbor->request(this,instance_me); + irequest = neighbor->request(this, instance_me); if (respa >= 1) { neighbor->requests[irequest]->respaouter = 1; - neighbor->requests[irequest]->respainner = 1; + neighbor->requests[irequest]->respainner = 1; } if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; @@ -689,19 +690,19 @@ void PairLJClass2CoulLong::init_style() // set rRESPA cutoffs - if (utils::strmatch(update->integrate_style,"^respa") && + if (utils::strmatch(update->integrate_style, "^respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = nullptr; + else + cut_respa = nullptr; // insure use of KSpace long-range solver, set g_ewald - if (force->kspace == nullptr) - error->all(FLERR,"Pair style requires a KSpace style"); + if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables - if (ncoultablebits) init_tables(cut_coul,cut_respa); + if (ncoultablebits) init_tables(cut_coul, cut_respa); } /* ---------------------------------------------------------------------- @@ -714,26 +715,25 @@ double PairLJClass2CoulLong::init_one(int i, int j) // mix distance via user-defined rule if (setflag[i][j] == 0) { - epsilon[i][j] = 2.0 * sqrt(epsilon[i][i]*epsilon[j][j]) * - pow(sigma[i][i],3.0) * pow(sigma[j][j],3.0) / - (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0)); - sigma[i][j] = - pow((0.5 * (pow(sigma[i][i],6.0) + pow(sigma[j][j],6.0))),1.0/6.0); - cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + epsilon[i][j] = 2.0 * sqrt(epsilon[i][i] * epsilon[j][j]) * pow(sigma[i][i], 3.0) * + pow(sigma[j][j], 3.0) / (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0)); + sigma[i][j] = pow((0.5 * (pow(sigma[i][i], 6.0) + pow(sigma[j][j], 6.0))), 1.0 / 6.0); + cut_lj[i][j] = mix_distance(cut_lj[i][i], cut_lj[j][j]); } - double cut = MAX(cut_lj[i][j],cut_coul); + double cut = MAX(cut_lj[i][j], cut_coul); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j],9.0); - lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + lj3[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + lj4[i][j] = 3.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); if (offset_flag && (cut_lj[i][j] > 0.0)) { double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = epsilon[i][j] * (2.0*pow(ratio,9.0) - 3.0*pow(ratio,6.0)); - } else offset[i][j] = 0.0; + offset[i][j] = epsilon[i][j] * (2.0 * pow(ratio, 9.0) - 3.0 * pow(ratio, 6.0)); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; lj1[j][i] = lj1[i][j]; @@ -744,8 +744,8 @@ double PairLJClass2CoulLong::init_one(int i, int j) // check interior rRESPA cutoff - if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + if (cut_respa && MIN(cut_lj[i][j], cut_coul) < cut_respa[3]) + error->all(FLERR, "Pair cutoff < Respa interior cutoff"); // compute I,J contribution to long-range tail correction // count total # of atoms of type I and J via Allreduce @@ -754,22 +754,21 @@ double PairLJClass2CoulLong::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); - double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j]; - double sig6 = sig3*sig3; - double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j]; - double rc6 = rc3*rc3; - etail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 3.0*rc3) / (3.0*rc6); - ptail_ij = 2.0*MY_PI*all[0]*all[1]*epsilon[i][j] * - sig6 * (sig3 - 2.0*rc3) / rc6; + double sig3 = sigma[i][j] * sigma[i][j] * sigma[i][j]; + double sig6 = sig3 * sig3; + double rc3 = cut_lj[i][j] * cut_lj[i][j] * cut_lj[i][j]; + double rc6 = rc3 * rc3; + etail_ij = + 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 3.0 * rc3) / (3.0 * rc6); + ptail_ij = 2.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 * (sig3 - 2.0 * rc3) / rc6; } return cut; @@ -783,14 +782,14 @@ void PairLJClass2CoulLong::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); + fwrite(&epsilon[i][j], sizeof(double), 1, fp); + fwrite(&sigma[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); } } } @@ -804,21 +803,21 @@ void PairLJClass2CoulLong::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -829,13 +828,13 @@ void PairLJClass2CoulLong::read_restart(FILE *fp) void PairLJClass2CoulLong::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); - fwrite(&ncoultablebits,sizeof(int),1,fp); - fwrite(&tabinner,sizeof(double),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); + fwrite(&ncoultablebits, sizeof(int), 1, fp); + fwrite(&tabinner, sizeof(double), 1, fp); } /* ---------------------------------------------------------------------- @@ -845,21 +844,21 @@ void PairLJClass2CoulLong::write_restart_settings(FILE *fp) void PairLJClass2CoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); - MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); - MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world); + MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- @@ -868,8 +867,7 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp) void PairLJClass2CoulLong::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]); } /* ---------------------------------------------------------------------- @@ -880,69 +878,68 @@ void PairLJClass2CoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut_lj[i][j]); + fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut_lj[i][j]); } /* ---------------------------------------------------------------------- */ -double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { - double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor; - double fraction,table,forcecoul,forcelj,phicoul,philj; + double r2inv, r, rinv, r3inv, r6inv, grij, expm2, t, erfc, prefactor; + double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + expm2 = exp(-grij * grij); + t = 1.0 / (1.0 + EWALD_P * grij); + erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { union_int_float_t rsq_lookup; rsq_lookup.f = rsq; itable = rsq_lookup.i & ncoulmask; itable >>= ncoulshiftbits; fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; - table = ftable[itable] + fraction*dftable[itable]; - forcecoul = atom->q[i]*atom->q[j] * table; + table = ftable[itable] + fraction * dftable[itable]; + forcecoul = atom->q[i] * atom->q[j] * table; if (factor_coul < 1.0) { - table = ctable[itable] + fraction*dctable[itable]; - prefactor = atom->q[i]*atom->q[j] * table; - forcecoul -= (1.0-factor_coul)*prefactor; + table = ctable[itable] + fraction * dctable[itable]; + prefactor = atom->q[i] * atom->q[j] * table; + forcecoul -= (1.0 - factor_coul) * prefactor; } } - } else forcecoul = 0.0; + } else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { rinv = sqrt(r2inv); - r3inv = r2inv*rinv; - r6inv = r3inv*r3inv; - forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]); - } else forcelj = 0.0; - fforce = (forcecoul + factor_lj*forcelj) * r2inv; + r3inv = r2inv * rinv; + r6inv = r3inv * r3inv; + forcelj = r6inv * (lj1[itype][jtype] * r3inv - lj2[itype][jtype]); + } else + forcelj = 0.0; + fforce = (forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor*erfc; + phicoul = prefactor * erfc; else { - table = etable[itable] + fraction*detable[itable]; - phicoul = atom->q[i]*atom->q[j] * table; + table = etable[itable] + fraction * detable[itable]; + phicoul = atom->q[i] * atom->q[j] * table; } - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - - offset[itype][jtype]; - eng += factor_lj*philj; + philj = r6inv * (lj3[itype][jtype] * r3inv - lj4[itype][jtype]) - offset[itype][jtype]; + eng += factor_lj * philj; } return eng; @@ -953,9 +950,9 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype, void *PairLJClass2CoulLong::extract(const char *str, int &dim) { dim = 0; - if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul; dim = 2; - if (strcmp(str,"epsilon") == 0) return (void *) epsilon; - if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + if (strcmp(str, "sigma") == 0) return (void *) sigma; return nullptr; } From 7fb684b802b66983b483ec5b9ffa4c0d024e56ae Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 1 Mar 2022 21:18:51 -0500 Subject: [PATCH 03/20] for 2d systems, rigid bodies always have a moment of inertia and no DOFs need to be subtracted --- src/RIGID/fix_rigid_nh_small.cpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index 105b24f83b..8b773a30ea 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -1192,13 +1192,7 @@ void FixRigidNHSmall::compute_dof() for (int k = 0; k < dimension; k++) if (fabs(b->inertia[k]) < EPSILON) nf_r--; } - } else if (dimension == 2) { - nf_r = nlocal_body; - for (int ibody = 0; ibody < nlocal_body; ibody++) { - Body *b = &body[ibody]; - if (fabs(b->inertia[2]) < EPSILON) nf_r--; - } - } + } else if (dimension == 2) nf_r = nlocal_body; double nf[2], nfall[2]; nf[0] = nf_t; From b4da01c23a6218fa62c22e7bff8c923bac462164 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 2 Mar 2022 22:19:28 -0500 Subject: [PATCH 04/20] reduce compiler warnings and need for private variables in OpenMP pragmas --- src/OPENMP/fix_rigid_nh_omp.cpp | 10 ++++------ src/OPENMP/pair_comb_omp.cpp | 9 ++++----- src/OPENMP/pair_reaxff_omp.cpp | 19 +++++++++---------- 3 files changed, 17 insertions(+), 21 deletions(-) diff --git a/src/OPENMP/fix_rigid_nh_omp.cpp b/src/OPENMP/fix_rigid_nh_omp.cpp index d95585e0ac..4e73eaa89b 100644 --- a/src/OPENMP/fix_rigid_nh_omp.cpp +++ b/src/OPENMP/fix_rigid_nh_omp.cpp @@ -246,12 +246,11 @@ void FixRigidNHOMP::compute_forces_and_torques() if (rstyle == SINGLE) { // we have just one rigid body. use OpenMP reduction to get sum[] double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; - int i; #if defined(_OPENMP) -#pragma omp parallel for LMP_DEFAULT_NONE private(i) reduction(+:s0,s1,s2,s3,s4,s5) +#pragma omp parallel for LMP_DEFAULT_NONE reduction(+:s0,s1,s2,s3,s4,s5) #endif - for (i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if (ibody < 0) continue; @@ -285,12 +284,11 @@ void FixRigidNHOMP::compute_forces_and_torques() for (int ib = 0; ib < nbody; ++ib) { double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; - int i; #if defined(_OPENMP) -#pragma omp parallel for LMP_DEFAULT_NONE private(i) LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5) +#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(ib) reduction(+:s0,s1,s2,s3,s4,s5) #endif - for (i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if (ibody != ib) continue; diff --git a/src/OPENMP/pair_comb_omp.cpp b/src/OPENMP/pair_comb_omp.cpp index 116db7194d..b3306d678c 100644 --- a/src/OPENMP/pair_comb_omp.cpp +++ b/src/OPENMP/pair_comb_omp.cpp @@ -384,7 +384,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) double PairCombOMP::yasu_char(double *qf_fix, int &igroup) { - int ii; double potal,fac11,fac11e; const double * const * const x = atom->x; @@ -401,7 +400,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) const int groupbit = group->bitmask[igroup]; qf = qf_fix; - for (ii = 0; ii < inum; ii++) { + for (int ii = 0; ii < inum; ii++) { const int i = ilist[ii]; if (mask[i] & groupbit) qf[i] = 0.0; @@ -417,9 +416,9 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) // loop over full neighbor list of my atoms #if defined(_OPENMP) -#pragma omp parallel for private(ii) LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e) +#pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(potal,fac11e) #endif - for (ii = 0; ii < inum; ii ++) { + for (int ii = 0; ii < inum; ii ++) { double fqi,fqij,fqji,fqjj,delr1[3]; double sr1,sr2,sr3; int mr1,mr2,mr3; @@ -533,7 +532,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) // sum charge force on each node and return it double eneg = 0.0; - for (ii = 0; ii < inum; ii++) { + for (int ii = 0; ii < inum; ii++) { const int i = ilist[ii]; if (mask[i] & groupbit) eneg += qf[i]; diff --git a/src/OPENMP/pair_reaxff_omp.cpp b/src/OPENMP/pair_reaxff_omp.cpp index 1f320f76fc..a329809c4e 100644 --- a/src/OPENMP/pair_reaxff_omp.cpp +++ b/src/OPENMP/pair_reaxff_omp.cpp @@ -402,8 +402,7 @@ int PairReaxFFOMP::estimate_reax_lists() int PairReaxFFOMP::write_reax_lists() { - int itr_i, itr_j, i, j, num_mynbrs; - int *jlist; + int num_mynbrs; double d_sqr, dist, cutoff_sqr; rvec dvec; @@ -425,19 +424,19 @@ int PairReaxFFOMP::write_reax_lists() num_nbrs = 0; - for (itr_i = 0; itr_i < numall; ++itr_i) { - i = ilist[itr_i]; + for (int itr_i = 0; itr_i < numall; ++itr_i) { + int i = ilist[itr_i]; num_nbrs_offset[i] = num_nbrs; num_nbrs += numneigh[i]; } #if defined(_OPENMP) #pragma omp parallel for schedule(dynamic,50) default(shared) \ - private(itr_i, itr_j, i, j, jlist, cutoff_sqr, num_mynbrs, d_sqr, dvec, dist) + private(cutoff_sqr, num_mynbrs, d_sqr, dvec, dist) #endif - for (itr_i = 0; itr_i < numall; ++itr_i) { - i = ilist[itr_i]; - jlist = firstneigh[i]; + for (int itr_i = 0; itr_i < numall; ++itr_i) { + int i = ilist[itr_i]; + auto jlist = firstneigh[i]; Set_Start_Index(i, num_nbrs_offset[i], far_nbrs); if (i < inum) @@ -447,8 +446,8 @@ int PairReaxFFOMP::write_reax_lists() num_mynbrs = 0; - for (itr_j = 0; itr_j < numneigh[i]; ++itr_j) { - j = jlist[itr_j]; + for (int itr_j = 0; itr_j < numneigh[i]; ++itr_j) { + int j = jlist[itr_j]; j &= NEIGHMASK; get_distance(x[j], x[i], &d_sqr, &dvec); From 04dbf307d06fba71268cd675e77d663ca59cbf63 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 3 Mar 2022 06:43:52 -0500 Subject: [PATCH 05/20] remove obsolete TODO item --- doc/src/Developer_org.rst | 6 ------ 1 file changed, 6 deletions(-) diff --git a/doc/src/Developer_org.rst b/doc/src/Developer_org.rst index ce36d9a590..63ff30abfe 100644 --- a/doc/src/Developer_org.rst +++ b/doc/src/Developer_org.rst @@ -252,12 +252,6 @@ follows: - The Timer class logs timing information, output at the end of a run. -.. TODO section on "Spatial decomposition and parallel operations" -.. diagram of 3d processor grid, brick vs. tiled. local vs. ghost -.. atoms, 6-way communication with pack/unpack functions, -.. PBC as part of the communication, forward and reverse communication -.. rendezvous communication, ring communication. - .. TODO section on "Fixes, Computes, and Variables" .. how and when data is computed and provided and how it is .. referenced. flags in Fix/Compute/Variable classes tell From 98f83f9b3aaec2c2e23592d5962fb38b4fd6e099 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 4 Mar 2022 06:25:16 -0500 Subject: [PATCH 06/20] fix typos and make output section more readable --- doc/src/compute_tally.rst | 60 ++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 26 deletions(-) diff --git a/doc/src/compute_tally.rst b/doc/src/compute_tally.rst index 313b54645c..45ebb42843 100644 --- a/doc/src/compute_tally.rst +++ b/doc/src/compute_tally.rst @@ -31,7 +31,7 @@ Syntax compute ID group-ID style group2-ID * ID, group-ID are documented in :doc:`compute ` command -* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or * or *pe/tally* or *pe/mol/tally* or *stress/tally* +* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or *pe/tally* or *pe/mol/tally* or *stress/tally* * group2-ID = group ID of second (or same) group Examples @@ -61,7 +61,7 @@ mechanism. Compute *pe/mol/tally* is one such style, that can - through using this mechanism - separately tally intermolecular and intramolecular energies. Something that would otherwise be impossible without integrating this as a core functionality into -the based classes of LAMMPS. +the base classes of LAMMPS. ---------- @@ -148,30 +148,38 @@ pairwise property computations. Output info """"""""""" -Compute *pe/tally* calculates a global scalar (the energy) and a per -atom scalar (the contributions of the single atom to the global -scalar). Compute *pe/mol/tally* calculates a global 4-element vector -containing (in this order): *evdwl* and *ecoul* for intramolecular pairs -and *evdwl* and *ecoul* for intermolecular pairs. Since molecules are -identified by their molecule IDs, the partitioning does not have to be -related to molecules, but the energies are tallied into the respective -slots depending on whether the molecule IDs of a pair are the same or -different. Compute *force/tally* calculates a global scalar (the force -magnitude) and a per atom 3-element vector (force contribution from -each atom). Compute *stress/tally* calculates a global scalar -(average of the diagonal elements of the stress tensor) and a per atom -vector (the 6 elements of stress tensor contributions from the -individual atom). As in :doc:`compute heat/flux `, -compute *heat/flux/tally* calculates a global vector of length 6, -where the first 3 components are the :math:`x`, :math:`y`, :math:`z` -components of the full heat flow vector, -and the next 3 components are the corresponding components -of just the convective portion of the flow, i.e. the -first term in the equation for :math:`\mathbf{Q}`. -Compute *heat/flux/virial/tally* calculates a global scalar (heat flow) -and a per atom 3-element vector -(contribution to the force acting over atoms in the first group -from individual atoms in both groups). +- Compute *pe/tally* calculates a global scalar (the energy) and a per + atom scalar (the contributions of the single atom to the global + scalar). + +- Compute *pe/mol/tally* calculates a global 4-element vector containing + (in this order): *evdwl* and *ecoul* for intramolecular pairs and + *evdwl* and *ecoul* for intermolecular pairs. Since molecules are + identified by their molecule IDs, the partitioning does not have to be + related to molecules, but the energies are tallied into the respective + slots depending on whether the molecule IDs of a pair are the same or + different. + +- Compute *force/tally* calculates a global scalar (the force magnitude) + and a per atom 3-element vector (force contribution from each atom). + +- Compute *stress/tally* calculates a global scalar + (average of the diagonal elements of the stress tensor) and a per atom + vector (the 6 elements of stress tensor contributions from the + individual atom). + +- As in :doc:`compute heat/flux `, + compute *heat/flux/tally* calculates a global vector of length 6, + where the first 3 components are the :math:`x`, :math:`y`, :math:`z` + components of the full heat flow vector, + and the next 3 components are the corresponding components + of just the convective portion of the flow, i.e. the + first term in the equation for :math:`\mathbf{Q}`. + +- Compute *heat/flux/virial/tally* calculates a global scalar (heat flow) + and a per atom 3-element vector + (contribution to the force acting over atoms in the first group + from individual atoms in both groups). Both the scalar and vector values calculated by this compute are "extensive". From f9c02b97350de06aa8a1cb17079f755c354efe8b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Mar 2022 10:18:00 -0500 Subject: [PATCH 07/20] clarify. improve typesetting --- doc/src/compute_momentum.rst | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/doc/src/compute_momentum.rst b/doc/src/compute_momentum.rst index c86901933b..40d2f527f5 100644 --- a/doc/src/compute_momentum.rst +++ b/doc/src/compute_momentum.rst @@ -23,11 +23,10 @@ Examples Description """"""""""" -Define a computation that calculates the translational momentum -of a group of particles. - -The momentum of each particles is computed as m v, where m and v are -the mass and velocity of the particle. +Define a computation that calculates the translational momentum *p* +of a group of particles. It is computed as the sum :math:`\vec{p} = \sum_i m_i \cdot \vec{v}_i` +over all particles in the compute group, where *m* and *v* are +the mass and velocity vector of the particle, respectively. Output info """"""""""" From a9eae51d27cafab26b0712fd296b97be64133bf3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Mar 2022 14:31:24 -0500 Subject: [PATCH 08/20] print name of non-existent compute/fix/variable in error message --- src/thermo.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index 27d74c58b6..f95b665e4e 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -905,7 +905,8 @@ void Thermo::parse_fields(char *str) if (argi.get_type() == ArgInfo::COMPUTE) { auto icompute = modify->get_compute_by_id(argi.get_name()); - if (!icompute) error->all(FLERR,"Could not find thermo custom compute ID"); + if (!icompute) + error->all(FLERR,"Could not find thermo custom compute ID: {}", argi.get_name()); if (argindex1[nfield] == 0 && icompute->scalar_flag == 0) error->all(FLERR,"Thermo compute does not compute scalar"); if (argindex1[nfield] > 0 && argindex2[nfield] == 0) { @@ -935,7 +936,7 @@ void Thermo::parse_fields(char *str) } else if (argi.get_type() == ArgInfo::FIX) { auto ifix = modify->get_fix_by_id(argi.get_name()); - if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID"); + if (!ifix) error->all(FLERR,"Could not find thermo custom fix ID: {}", argi.get_name()); if (argindex1[nfield] == 0 && ifix->scalar_flag == 0) error->all(FLERR,"Thermo fix does not compute scalar"); if (argindex1[nfield] > 0 && argindex2[nfield] == 0) { @@ -961,7 +962,7 @@ void Thermo::parse_fields(char *str) } else if (argi.get_type() == ArgInfo::VARIABLE) { int n = input->variable->find(argi.get_name()); if (n < 0) - error->all(FLERR,"Could not find thermo custom variable name"); + error->all(FLERR,"Could not find thermo custom variable name: {}", argi.get_name()); if (argindex1[nfield] == 0 && input->variable->equalstyle(n) == 0) error->all(FLERR,"Thermo custom variable is not equal-style variable"); if (argindex1[nfield] && input->variable->vectorstyle(n) == 0) From 8945050423bc76cd18891ee541e5e8eccbbb6151 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Mar 2022 16:47:31 -0500 Subject: [PATCH 09/20] cosmetic --- src/MANYBODY/pair_sw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index bac65a2464..239b125f93 100644 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -292,7 +292,7 @@ void PairSW::read_file(char *file) if (comm->me == 0) { PotentialFileReader reader(lmp, file, "sw", unit_convert_flag); - char * line; + char *line; // transparently convert units for supported conversions From 83f08ff8eff2cbbdcaa3fa6b6fe0a27c64267d51 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 5 Mar 2022 17:03:19 -0500 Subject: [PATCH 10/20] small update and correct broken link --- doc/src/Developer_notes.rst | 23 +++++++++++------------ doc/src/Developer_utils.rst | 2 +- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/doc/src/Developer_notes.rst b/doc/src/Developer_notes.rst index cd3c51fed6..e58dec94ac 100644 --- a/doc/src/Developer_notes.rst +++ b/doc/src/Developer_notes.rst @@ -25,18 +25,17 @@ following benefits: - transparent support for translating unsupported UTF-8 characters to their ASCII equivalents (the text-to-value conversion functions **only** accept ASCII characters) -In most cases (e.g. potential files) the same data is needed on all -MPI ranks. Then it is best to do the reading and parsing only on MPI -rank 0, and communicate the data later with one or more -``MPI_Bcast()`` calls. For reading generic text and potential -parameter files the custom classes :cpp:class:`TextFileReader -` and :cpp:class:`PotentialFileReader -` are available. They allow reading -the file as individual lines for which they can return a tokenizer -class (see below) for parsing the line. Or they can return blocks of -numbers as a vector directly. The documentation on `File reader -classes `_ contains an example for a typical -case. +In most cases (e.g. potential files) the same data is needed on all MPI +ranks. Then it is best to do the reading and parsing only on MPI rank +0, and communicate the data later with one or more ``MPI_Bcast()`` +calls. For reading generic text and potential parameter files the +custom classes :cpp:class:`TextFileReader ` +and :cpp:class:`PotentialFileReader ` +are available. They allow reading the file as individual lines for which +they can return a tokenizer class (see below) for parsing the line. Or +they can return blocks of numbers as a vector directly. The +documentation on :ref:`File reader classes ` +contains an example for a typical case. When reading per-atom data, the data on each line of the file usually needs to include an atom ID so it can be associated with a particular diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index a9df85c899..720eececcc 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -396,7 +396,7 @@ A typical code segment would look like this: ---------- -.. file-reader-classes: +.. _file-reader-classes: File reader classes ------------------- From a2bf40df0a81b3d1429e1ffa1a9fb8a2520fd56a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 6 Mar 2022 11:45:26 -0500 Subject: [PATCH 11/20] correct mini-ToC links and add missing entry --- doc/src/Build_settings.rst | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/doc/src/Build_settings.rst b/doc/src/Build_settings.rst index b68313aaed..6053bdbf8a 100644 --- a/doc/src/Build_settings.rst +++ b/doc/src/Build_settings.rst @@ -4,15 +4,15 @@ Optional build settings LAMMPS can be built with several optional settings. Each sub-section explain how to do this for building both with CMake and make. -* :ref:`C++11 standard compliance ` when building all of LAMMPS -* :ref:`FFT library ` for use with the :doc:`kspace_style pppm ` command -* :ref:`Size of LAMMPS integer types ` -* :ref:`Read or write compressed files ` -* :ref:`Output of JPG and PNG files ` via the :doc:`dump image ` command -* :ref:`Output of movie files ` via the :doc:`dump_movie ` command -* :ref:`Memory allocation alignment ` -* :ref:`Workaround for long long integers ` -* :ref:`Error handling exceptions ` when using LAMMPS as a library +* `C++11 standard compliance`_ when building all of LAMMPS +* `FFT library`_ for use with the :doc:`kspace_style pppm ` command +* `Size of LAMMPS integer types and size limits`_ +* `Read or write compressed files`_ +* `Output of JPG, PNG, and move files` via the :doc:`dump image ` or :doc:`dump movie ` commands +* `Memory allocation alignment`_ +* `Workaround for long long integers`_ +* `Exception handling when using LAMMPS as a library`_ to capture errors +* `Trigger selected floating-point exceptions`_ ---------- From 6edfb495275e738cc224e168d0f5c7a96b0bfb97 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Mar 2022 14:37:01 -0500 Subject: [PATCH 12/20] update .gitignore for recent additions --- src/.gitignore | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/.gitignore b/src/.gitignore index a9eef4b371..8a18f5b29f 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1,6 +1,7 @@ /Makefile.package /Makefile.package.settings /MAKE/MINE +/ADIOS/Makefile.lammps /Make.py.last /lmp_* @@ -267,6 +268,8 @@ /fix_drag.h /fix_numdiff.cpp /fix_numdiff.h +/fix_numdiff_virial.cpp +/fix_numdiff_virial.h /fix_spring_rg.cpp /fix_spring_rg.h /fix_temp_csld.cpp @@ -930,6 +933,8 @@ /msm_cg.h /neb.cpp /neb.h +/netcdf_units.cpp +/netcdf_units.h /pair_adp.cpp /pair_adp.h /pair_agni.cpp @@ -1065,6 +1070,8 @@ /pair_hdnnp.h /pair_ilp_graphene_hbn.cpp /pair_ilp_graphene_hbn.h +/pair_ilp_tmd.cpp +/pair_ilp_tmd.h /pair_kolmogorov_crespi_full.cpp /pair_kolmogorov_crespi_full.h /pair_kolmogorov_crespi_z.cpp @@ -1171,8 +1178,8 @@ /pair_nm_cut_coul_cut.h /pair_nm_cut_coul_long.cpp /pair_nm_cut_coul_long.h -/pait_nm_cut_split.cpp -/pait_nm_cut_split.h +/pair_nm_cut_split.cpp +/pair_nm_cut_split.h /pair_oxdna_*.cpp /pair_oxdna_*.h /pair_oxdna2_*.cpp @@ -1198,6 +1205,8 @@ /pair_rebo.h /pair_resquared.cpp /pair_resquared.h +/pair_saip_metal.cpp +/pair_saip_metal.h /pair_sdpd_taitwater_isothermal.cpp /pair_sdpd_taitwater_isothermal.h /pair_sph_heatconduction.cpp From 87613bb1069dc8472f4ad13837ac4e8b63a376ca Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Mar 2022 15:34:00 -0500 Subject: [PATCH 13/20] rename doc file to be more generic and match the examples folder --- doc/src/Commands_compute.rst | 6 +++--- doc/src/compute.rst | 6 +++--- doc/src/compute_stress_mop.rst | 4 ++-- ...pute_stress_cartesian.rst => compute_stress_profile.rst} | 0 4 files changed, 8 insertions(+), 8 deletions(-) rename doc/src/{compute_stress_cartesian.rst => compute_stress_profile.rst} (100%) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 590b3d2ea8..6ad5ed4435 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -142,11 +142,11 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`sph/t/atom ` * :doc:`spin ` * :doc:`stress/atom ` - * :doc:`stress/cartesian ` - * :doc:`stress/cylinder ` + * :doc:`stress/cartesian ` + * :doc:`stress/cylinder ` * :doc:`stress/mop ` * :doc:`stress/mop/profile ` - * :doc:`stress/spherical ` + * :doc:`stress/spherical ` * :doc:`stress/tally ` * :doc:`tdpd/cc/atom ` * :doc:`temp (k) ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6973559d16..2768543179 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -288,11 +288,11 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`sph/t/atom ` - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms * :doc:`spin ` - magnetic quantities for a system of atoms having spins * :doc:`stress/atom ` - stress tensor for each atom -* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates -* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates +* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates +* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates * :doc:`stress/mop ` - normal components of the local stress tensor using the method of planes * :doc:`stress/mop/profile ` - profile of the normal components of the local stress tensor using the method of planes -* :doc:`stress/spherical ` - stress tensor in spherical coordinates +* :doc:`stress/spherical ` - stress tensor in spherical coordinates * :doc:`stress/tally ` - stress between two groups of atoms via the tally callback mechanism * :doc:`tdpd/cc/atom ` - per-atom chemical concentration of a specified species for each tDPD particle * :doc:`temp ` - temperature of group of atoms diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 18b9e703c5..d439e18d34 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -73,7 +73,7 @@ command, since those are contributions to the global system pressure. NOTE 3: The local stress profile generated by compute *stress/mop/profile* is similar to that obtained by compute -:doc:`stress/cartesian `. +:doc:`stress/cartesian `. A key difference is that compute *stress/mop/profile* considers particles crossing a set of planes, @@ -122,7 +122,7 @@ intra-molecular interactions, and long range (kspace) interactions. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` +:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` Default """"""" diff --git a/doc/src/compute_stress_cartesian.rst b/doc/src/compute_stress_profile.rst similarity index 100% rename from doc/src/compute_stress_cartesian.rst rename to doc/src/compute_stress_profile.rst From 374917fb6f09cf7445f6b4b430bc8cc807bf0012 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 7 Mar 2022 21:42:02 -0500 Subject: [PATCH 14/20] tweak platform tests for CPU time to avoid bogus failures on windows --- unittest/utils/test_platform.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/unittest/utils/test_platform.cpp b/unittest/utils/test_platform.cpp index 37e749b9be..a809b61458 100644 --- a/unittest/utils/test_platform.cpp +++ b/unittest/utils/test_platform.cpp @@ -21,7 +21,7 @@ TEST(Platform, clock) // spend some time computing pi constexpr double known_pi = 3.141592653589793238462643; - constexpr int n = 10000000; + constexpr int n = 100000000; constexpr double h = 1.0 / (double)n; double my_pi = 0.0, x; for (int i = 0; i < n; ++i) { @@ -34,7 +34,8 @@ TEST(Platform, clock) ASSERT_NEAR(my_pi, known_pi, 1e-12); ASSERT_GT(wt_used, 1e-4); - ASSERT_GT(ct_used, 1e-4); + // windows sometimes incorrectly reports used CPU time as 0.0 + if (ct_used != 0.0) ASSERT_GT(ct_used, 1e-4); } TEST(Platform, putenv_unsetenv) From c4b23bd7e55a60665ef1baf8a6d482fdd88de581 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Mar 2022 19:21:03 -0500 Subject: [PATCH 15/20] programming style updates use newer/simpler APIs --- src/LATTE/fix_latte.cpp | 61 ++++++++++------------------------------- 1 file changed, 15 insertions(+), 46 deletions(-) diff --git a/src/LATTE/fix_latte.cpp b/src/LATTE/fix_latte.cpp index e1d0224e62..63b24b1cdd 100644 --- a/src/LATTE/fix_latte.cpp +++ b/src/LATTE/fix_latte.cpp @@ -17,20 +17,18 @@ ------------------------------------------------------------------------- */ #include "fix_latte.h" -#include -#include + #include "atom.h" #include "comm.h" -#include "update.h" -#include "neighbor.h" -#include "domain.h" -#include "force.h" -#include "neigh_request.h" -#include "neigh_list.h" -#include "modify.h" #include "compute.h" -#include "memory.h" +#include "domain.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -82,13 +80,12 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[3],"NULL") != 0) { coulomb = 1; - error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation " - "of a Coulomb potential"); + error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation of a Coulomb potential"); id_pe = utils::strdup(arg[3]); - int ipe = modify->find_compute(id_pe); - if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID"); - if (modify->compute[ipe]->peatomflag == 0) + c_pe = modify->get_compute_by_id(id_pe); + if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); + if (c_pe->peatomflag == 0) error->all(FLERR,"Fix latte compute ID does not compute pe/atom"); } @@ -105,7 +102,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) : FixLatte::~FixLatte() { - delete [] id_pe; + delete[] id_pe; memory->destroy(qpotential); memory->destroy(flatte); } @@ -134,9 +131,8 @@ void FixLatte::init() if (atom->q_flag == 0 || force->pair == nullptr || force->kspace == nullptr) error->all(FLERR,"Fix latte cannot compute Coulomb potential"); - int ipe = modify->find_compute(id_pe); - if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID"); - c_pe = modify->compute[ipe]; + c_pe = modify->get_compute_by_id(id_pe); + if (!c_pe) error->all(FLERR,"Could not find fix latte compute ID {}", id_pe); } // must be fully periodic or fully non-periodic @@ -153,33 +149,6 @@ void FixLatte::init() memory->create(qpotential,atom->nlocal,"latte:qpotential"); memory->create(flatte,atom->nlocal,3,"latte:flatte"); } - - /* - // warn if any integrate fix comes after this one - // is it actually necessary for q(n) update to come after x,v update ?? - - int after = 0; - int flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(id,modify->fix[i]->id) == 0) after = 1; - else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1; - } - if (flag && comm->me == 0) - error->warning(FLERR,"Fix latte should come after all other " - "integration fixes"); - */ - - /* - // need a full neighbor list - // could we use a half list? - // perpetual list, built whenever re-neighboring occurs - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->fix = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->full = 1; - */ } /* ---------------------------------------------------------------------- */ From 8e85481afa6139017c2da7742987b2269827a6f8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Mar 2022 14:00:34 -0500 Subject: [PATCH 16/20] use MathConst::MY_PI instead of inferring it from acos(-1). --- src/MANYBODY/pair_tersoff_table.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/MANYBODY/pair_tersoff_table.cpp b/src/MANYBODY/pair_tersoff_table.cpp index 0c11760737..12733e1c60 100644 --- a/src/MANYBODY/pair_tersoff_table.cpp +++ b/src/MANYBODY/pair_tersoff_table.cpp @@ -27,6 +27,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "math_const.h" #include "memory.h" #include "neigh_list.h" #include "neigh_request.h" @@ -37,6 +38,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; #define MAXLINE 1024 #define DELTA 4 @@ -543,7 +545,6 @@ void PairTersoffTable::allocateGrids() double deltaArgumentCutoffFunction, deltaArgumentExponential, deltaArgumentBetaZetaPower; double deltaArgumentGtetaFunction; double r, minMu, maxLambda, maxCutoff; - double const PI=acos(-1.0); deallocateGrids(); @@ -652,8 +653,8 @@ void PairTersoffTable::allocateGrids() } for (l = numGridPointsOneCutoffFunction; l < numGridPointsCutoffFunction; l++) { - cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (PI * (r - cutoffR)/(cutoffS-cutoffR)) ; - cutoffFunctionDerived[i][j][l] = -0.5 * PI * sin (PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR) ; + cutoffFunction[i][j][l] = 0.5 + 0.5 * cos (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) ; + cutoffFunctionDerived[i][j][l] = -0.5 * MY_PI * sin (MY_PI * (r - cutoffR)/(cutoffS-cutoffR)) / (cutoffS-cutoffR); r += deltaArgumentCutoffFunction; } } From 0fcf40c48e713ef50c3eb46fa207d27fa91a648b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Mar 2022 14:14:32 -0500 Subject: [PATCH 17/20] use MathConst::MY_PI more consistently and benefit from it being a constexpr --- src/DPD-REACT/pair_multi_lucy.cpp | 7 ++--- src/DPD-REACT/pair_multi_lucy_rx.cpp | 10 +++---- src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp | 35 ++++++++++++------------ src/MEAM/meam_setup_param.cpp | 1 - 4 files changed, 24 insertions(+), 29 deletions(-) diff --git a/src/DPD-REACT/pair_multi_lucy.cpp b/src/DPD-REACT/pair_multi_lucy.cpp index 0f54a57870..41c9d9fb66 100644 --- a/src/DPD-REACT/pair_multi_lucy.cpp +++ b/src/DPD-REACT/pair_multi_lucy.cpp @@ -37,6 +37,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -104,7 +105,6 @@ void PairMultiLucy::compute(int eflag, int vflag) int nlocal = atom->nlocal; int newton_pair = force->newton_pair; - double pi = MathConst::MY_PI; double A_i; double A_j; double fraction_i,fraction_j; @@ -198,7 +198,7 @@ void PairMultiLucy::compute(int eflag, int vflag) evdwl = tb->e[itable] + fraction_i*tb->de[itable]; } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy"); - evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0; + evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0; if (evflag) ev_tally(0,0,nlocal,newton_pair, evdwl,0.0,0.0,0.0,0.0,0.0); @@ -733,7 +733,6 @@ void PairMultiLucy::computeLocalDensity() numneigh = list->numneigh; firstneigh = list->firstneigh; - double pi = MathConst::MY_PI; double *rho = atom->rho; // zero out density @@ -766,7 +765,7 @@ void PairMultiLucy::computeLocalDensity() if (rsq < cutsq[itype][jtype]) { double rcut = sqrt(cutsq[itype][jtype]); - factor= (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut); + factor= (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut)*(1.0-sqrt(rsq)/rcut); rho[i] += factor; if (newton_pair || j < nlocal) { rho[j] += factor; diff --git a/src/DPD-REACT/pair_multi_lucy_rx.cpp b/src/DPD-REACT/pair_multi_lucy_rx.cpp index 64b2b4c7d5..8b348810fd 100644 --- a/src/DPD-REACT/pair_multi_lucy_rx.cpp +++ b/src/DPD-REACT/pair_multi_lucy_rx.cpp @@ -39,6 +39,7 @@ #include using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -127,7 +128,6 @@ void PairMultiLucyRX::compute(int eflag, int vflag) double *uCG = atom->uCG; double *uCGnew = atom->uCGnew; - double pi = MathConst::MY_PI; double A_i, A_j; double fraction_i,fraction_j; int jtable; @@ -276,7 +276,7 @@ void PairMultiLucyRX::compute(int eflag, int vflag) evdwl = tb->e[itable] + fraction_i*tb->de[itable]; } else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx"); - evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0; + evdwl *=(MY_PI*cutsq[itype][itype]*cutsq[itype][itype])/84.0; evdwlOld = mixWtSite1old_i*evdwl; evdwl = mixWtSite1_i*evdwl; @@ -866,15 +866,13 @@ void PairMultiLucyRX::computeLocalDensity() const int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - const double pi = MathConst::MY_PI; - const bool newton_pair = force->newton_pair; const bool one_type = (atom->ntypes == 1); // Special cut-off values for when there's only one type. const double cutsq_type11 = cutsq[1][1]; const double rcut_type11 = sqrt(cutsq_type11); - const double factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11); + const double factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11); double *rho = atom->rho; @@ -925,7 +923,7 @@ void PairMultiLucyRX::computeLocalDensity() const double rcut = sqrt(cutsq[itype][jtype]); const double tmpFactor = 1.0-sqrt(rsq)/rcut; const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; - const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; + const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; rho_i += factor; if (newton_pair || j < nlocal) rho[j] += factor; diff --git a/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp b/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp index 8c558db32c..c2942288ff 100644 --- a/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp +++ b/src/KOKKOS/pair_multi_lucy_rx_kokkos.cpp @@ -23,20 +23,24 @@ ------------------------------------------------------------------------------------------- */ #include "pair_multi_lucy_rx_kokkos.h" + +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "kokkos.h" +#include "math_const.h" +#include "math_const.h" +#include "memory_kokkos.h" +#include "neigh_list.h" +#include "neigh_request.h" + #include #include -#include "math_const.h" -#include "atom_kokkos.h" -#include "force.h" -#include "comm.h" -#include "neigh_list.h" -#include "memory_kokkos.h" -#include "error.h" -#include "atom_masks.h" -#include "neigh_request.h" -#include "kokkos.h" using namespace LAMMPS_NS; +using MathConst::MY_PI; enum{NONE,RLINEAR,RSQ}; @@ -282,7 +286,6 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXCompute::operator()(TagPairMultiLucyRXCompute()() = 3; - evdwl *=(pi*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0; + evdwl *=(MY_PI*d_cutsq(itype,itype)*d_cutsq(itype,itype))/84.0; evdwlOld = mixWtSite1old_i*evdwl; evdwl = mixWtSite1_i*evdwl; @@ -458,15 +461,13 @@ void PairMultiLucyRXKokkos::computeLocalDensity() d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; - const double pi = MathConst::MY_PI; - const bool newton_pair = force->newton_pair; const bool one_type = (atom->ntypes == 1); // Special cut-off values for when there's only one type. cutsq_type11 = cutsq[1][1]; rcut_type11 = sqrt(cutsq_type11); - factor_type11 = 84.0/(5.0*pi*rcut_type11*rcut_type11*rcut_type11); + factor_type11 = 84.0/(5.0*MY_PI*rcut_type11*rcut_type11*rcut_type11); // zero out density int m = nlocal; @@ -548,8 +549,6 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXComputeLoca const int itype = type[i]; const int jnum = d_numneigh[i]; - const double pi = MathConst::MY_PI; - for (int jj = 0; jj < jnum; jj++) { const int j = (d_neighbors(i,jj) & NEIGHMASK); const int jtype = type[j]; @@ -574,7 +573,7 @@ void PairMultiLucyRXKokkos::operator()(TagPairMultiLucyRXComputeLoca const double rcut = sqrt(d_cutsq(itype,jtype)); const double tmpFactor = 1.0-sqrt(rsq)/rcut; const double tmpFactor4 = tmpFactor*tmpFactor*tmpFactor*tmpFactor; - const double factor = (84.0/(5.0*pi*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; + const double factor = (84.0/(5.0*MY_PI*rcut*rcut*rcut))*(1.0+3.0*sqrt(rsq)/(2.0*rcut))*tmpFactor4; rho_i_contrib += factor; if ((NEIGHFLAG==HALF || NEIGHFLAG==HALFTHREAD) && (NEWTON_PAIR || j < nlocal)) a_rho[j] += factor; diff --git a/src/MEAM/meam_setup_param.cpp b/src/MEAM/meam_setup_param.cpp index b15b70e263..5612936c20 100644 --- a/src/MEAM/meam_setup_param.cpp +++ b/src/MEAM/meam_setup_param.cpp @@ -233,7 +233,6 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3 // 21 = theta // see alloyparams(void) in meam_setup_done.cpp case 21: - // const double PI = 3.141592653589793238463; meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; From 0828ae66a077a21bd623c265a14b887512d84a7e Mon Sep 17 00:00:00 2001 From: keisukefujii Date: Fri, 11 Mar 2022 05:48:34 +0900 Subject: [PATCH 18/20] Fix a typo in the doc --- doc/src/Packages_details.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index 1b4b82a5e7..396ecd6442 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -9,7 +9,7 @@ gives links to documentation, example scripts, and pictures/movies (if available) that illustrate use of the package. The majority of packages can be included in a LAMMPS build with a -single setting (``-D PGK_=on`` for CMake) or command +single setting (``-D PKG_=on`` for CMake) or command (``make yes-`` for make). See the :doc:`Build package ` page for more info. A few packages may require additional steps; this is indicated in the descriptions below. The :doc:`Build extras ` From ca76ff360fda4c1641abff5ef4924d7b8a7042fa Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Mar 2022 21:17:58 -0500 Subject: [PATCH 19/20] augment cmake library search path to include the CUDA stubs library folder this will help configuring and compiling LAMMPS with CUDA support on machines where there is no CUDA driver installed --- cmake/Modules/Packages/GPU.cmake | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/cmake/Modules/Packages/GPU.cmake b/cmake/Modules/Packages/GPU.cmake index fe15917f47..aec8887c30 100644 --- a/cmake/Modules/Packages/GPU.cmake +++ b/cmake/Modules/Packages/GPU.cmake @@ -30,7 +30,15 @@ file(GLOB GPU_LIB_SOURCES ${LAMMPS_LIB_SOURCE_DIR}/gpu/[^.]*.cpp) file(MAKE_DIRECTORY ${LAMMPS_LIB_BINARY_DIR}/gpu) if(GPU_API STREQUAL "CUDA") - find_package(CUDA REQUIRED) + find_package(CUDA QUIET) + # augment search path for CUDA toolkit libraries to include the stub versions. Needed to find libcuda.so on machines without a CUDA driver installation + if(CUDA_FOUND) + set(CMAKE_LIBRARY_PATH "${CUDA_TOOLKIT_ROOT_DIR}/lib64/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib/stubs;${CUDA_TOOLKIT_ROOT_DIR}/lib64;${CUDA_TOOLKIT_ROOT_DIR}/lib;${CMAKE_LIBRARY_PATH}") + find_package(CUDA REQUIRED) + else() + message(FATAL_ERROR "CUDA Toolkit not found") + endif() + find_program(BIN2C bin2c) if(NOT BIN2C) message(FATAL_ERROR "Could not find bin2c, use -DBIN2C=/path/to/bin2c to help cmake finding it.") From 3a1921b3ed59fa9f159e99df7a5708e0c6ed139e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 10 Mar 2022 23:01:23 -0500 Subject: [PATCH 20/20] correct SHAKE bond stats output and avoid division by zero --- src/RIGID/fix_shake.cpp | 8 ++++++-- src/RIGID/fix_shake.h | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index 735a36ce55..f012789414 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -52,7 +52,7 @@ 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_ave(nullptr), b_max(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) @@ -198,6 +198,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : 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]; @@ -291,6 +292,7 @@ FixShake::~FixShake() delete[] b_count; delete[] b_count_all; delete[] b_atom; + delete[] b_atom_all; delete[] b_ave; delete[] b_ave_all; delete[] b_max; @@ -2472,6 +2474,7 @@ void FixShake::stats() 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++) { a_count[i] = 0; @@ -2547,6 +2550,7 @@ 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); @@ -2566,7 +2570,7 @@ void FixShake::stats() 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[i]); + b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]); } for (i = 1; i < na; i++) { const auto acnt = a_count_all[i]; diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 8a076afaba..677cdfa942 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -108,7 +108,7 @@ class FixShake : public Fix { int nlist, maxlist; // size and max-size of list // stat quantities - int *b_count, *b_count_all, *b_atom; // counts for each bond type + int *b_count, *b_count_all, *b_atom, *b_atom_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