diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp index 4cc05c8c90..89ce8da9a6 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/DIELECTRIC/atom_vec_dielectric.cpp @@ -38,7 +38,7 @@ static const char cite_user_dielectric_package[] = /* ---------------------------------------------------------------------- */ -AtomVecDielectric::AtomVecDielectric(LAMMPS *lmp) : AtomVec(lmp) +AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp) { if (lmp->citeme) lmp->citeme->add(cite_user_dielectric_package); diff --git a/src/DIELECTRIC/compute_efield_atom.cpp b/src/DIELECTRIC/compute_efield_atom.cpp index 38dca98cb4..24c29321ed 100644 --- a/src/DIELECTRIC/compute_efield_atom.cpp +++ b/src/DIELECTRIC/compute_efield_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/ Sandia National Laboratories @@ -39,10 +38,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -ComputeEfieldAtom::ComputeEfieldAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), efield(nullptr) +ComputeEfieldAtom::ComputeEfieldAtom(LAMMPS *_lmp, int narg, char **arg) : + Compute(_lmp, narg, arg), efield(nullptr) { - if (narg < 3) error->all(FLERR,"Illegal compute efield/atom command"); + if (narg < 3) error->all(FLERR, "Illegal compute efield/atom command"); peratom_flag = 1; size_peratom_cols = 3; @@ -58,9 +57,12 @@ ComputeEfieldAtom::ComputeEfieldAtom(LAMMPS *lmp, int narg, char **arg) : } else { int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; - else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1; - else error->all(FLERR,"Illegal compute efield/atom command"); + if (strcmp(arg[iarg], "pair") == 0) + pairflag = 1; + else if (strcmp(arg[iarg], "kspace") == 0) + kspaceflag = 1; + else + error->all(FLERR, "Illegal compute efield/atom command"); iarg++; } } @@ -81,7 +83,7 @@ ComputeEfieldAtom::~ComputeEfieldAtom() void ComputeEfieldAtom::init() { - if (!atom->q_flag) error->all(FLERR,"compute efield/atom requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "compute efield/atom requires atom attribute q"); if (!force->kspace) kspaceflag = 0; } @@ -89,28 +91,30 @@ void ComputeEfieldAtom::init() void ComputeEfieldAtom::setup() { - if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric/omp") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/msm/dielectric") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric/omp") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"coul/long/dielectric") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else if (strcmp(force->pair_style,"coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast(force->pair))->efield; - else error->all(FLERR,"Compute efield/atom not supported by pair style"); + if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "coul/long/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else + error->all(FLERR, "Compute efield/atom not supported by pair style"); if (force->kspace) { - if (strcmp(force->kspace_style,"pppm/dielectric") == 0) - efield_kspace = (dynamic_cast(force->kspace))->efield; - else if (strcmp(force->kspace_style,"msm/dielectric") == 0) - efield_kspace = (dynamic_cast(force->kspace))->efield; - else error->all(FLERR,"Compute efield/atom not supported by kspace style"); + if (strcmp(force->kspace_style, "pppm/dielectric") == 0) + efield_kspace = (dynamic_cast(force->kspace))->efield; + else if (strcmp(force->kspace_style, "msm/dielectric") == 0) + efield_kspace = (dynamic_cast(force->kspace))->efield; + else + error->all(FLERR, "Compute efield/atom not supported by kspace style"); kspaceflag = 1; } @@ -122,11 +126,11 @@ void ComputeEfieldAtom::setup() void ComputeEfieldAtom::compute_peratom() { - int i,j; + int i, j; invoked_peratom = update->ntimestep; if (update->vflag_atom != invoked_peratom) - error->all(FLERR,"Per-atom virial was not tallied on needed timestep"); + error->all(FLERR, "Per-atom virial was not tallied on needed timestep"); // grow local stress array if necessary // needs to be atom->nmax in length @@ -134,7 +138,7 @@ void ComputeEfieldAtom::compute_peratom() if (atom->nmax > nmax) { memory->destroy(efield); nmax = atom->nmax; - memory->create(efield,nmax,3,"stress/atom:efield"); + memory->create(efield, nmax, 3, "stress/atom:efield"); array_atom = efield; } @@ -144,7 +148,7 @@ void ComputeEfieldAtom::compute_peratom() // ntotal includes ghosts if either newton flag is set // KSpace includes ghosts if tip4pflag is set - double* q = atom->q; + double *q = atom->q; int nlocal = atom->nlocal; int npair = nlocal; int ntotal = nlocal; @@ -156,8 +160,7 @@ void ComputeEfieldAtom::compute_peratom() // clear local stress array for (i = 0; i < ntotal; i++) - for (j = 0; j < 3; j++) - efield[i][j] = 0.0; + for (j = 0; j < 3; j++) efield[i][j] = 0.0; // add in per-atom contributions from each force @@ -170,14 +173,12 @@ void ComputeEfieldAtom::compute_peratom() if (kspaceflag && force->kspace) { for (i = 0; i < nkspace; i++) - for (j = 0; j < 3; j++) - efield[i][j] += efield_kspace[i][j]; + for (j = 0; j < 3; j++) efield[i][j] += efield_kspace[i][j]; } // communicate ghost efield between neighbor procs - if (force->newton || (force->kspace && force->kspace->tip4pflag)) - comm->reverse_comm(this); + if (force->newton || (force->kspace && force->kspace->tip4pflag)) comm->reverse_comm(this); // zero efield of atoms not in group // only do this after comm since ghost contributions must be included @@ -192,12 +193,11 @@ void ComputeEfieldAtom::compute_peratom() } } - /* ---------------------------------------------------------------------- */ int ComputeEfieldAtom::pack_reverse_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -213,7 +213,7 @@ int ComputeEfieldAtom::pack_reverse_comm(int n, int first, double *buf) void ComputeEfieldAtom::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i, j, m; m = 0; for (i = 0; i < n; i++) { @@ -230,6 +230,6 @@ void ComputeEfieldAtom::unpack_reverse_comm(int n, int *list, double *buf) double ComputeEfieldAtom::memory_usage() { - double bytes = nmax*3 * sizeof(double); + double bytes = nmax * 3 * sizeof(double); return bytes; } diff --git a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp index 62fa7e8641..4fdb1045b6 100644 --- a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp +++ b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp @@ -60,19 +60,17 @@ using namespace LAMMPS_NS; using namespace FixConst; -using namespace MathConst; - -//#define _POLARIZE_DEBUG +using MathConst::MY_PI; /* ---------------------------------------------------------------------- */ -FixPolarizeBEMGMRES::FixPolarizeBEMGMRES(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), q_backup(nullptr), c(nullptr), g(nullptr), h(nullptr), r(nullptr), s(nullptr), v(nullptr), - y(nullptr) +FixPolarizeBEMGMRES::FixPolarizeBEMGMRES(LAMMPS *_lmp, int narg, char **arg) : + Fix(_lmp, narg, arg), q_backup(nullptr), c(nullptr), g(nullptr), h(nullptr), r(nullptr), + s(nullptr), v(nullptr), y(nullptr) { if (narg < 5) error->all(FLERR, "Illegal fix polarize/bem/gmres command"); - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Fix polarize requires atom style dielectric"); // parse required arguments @@ -236,9 +234,7 @@ void FixPolarizeBEMGMRES::init() } if (comm->me == 0) - utils::logmesg(lmp, - "GMRES solver for {} induced charges " - "using maximum {} q-vectors\n", + utils::logmesg(lmp, "GMRES solver for {} induced charges using maximum {} q-vectors\n", num_induced_charges, mr); } @@ -249,32 +245,32 @@ void FixPolarizeBEMGMRES::setup(int /*vflag*/) // check if the pair styles in use are compatible if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else error->all(FLERR, "Pair style not compatible with fix polarize/bem/gmres"); if (kspaceflag) { if (force->kspace) { if (strcmp(force->kspace_style, "pppm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else if (strcmp(force->kspace_style, "msm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else error->all(FLERR, "Kspace style not compatible with fix polarize/bem/gmres"); } else @@ -369,8 +365,7 @@ void FixPolarizeBEMGMRES::compute_induced_charges() Ey += efield_kspace[i][1]; Ez += efield_kspace[i][2]; } - double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / - epsilon[i]; + double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i]; double sigma_f = q_real[i] / area[i]; buffer[idx] = (1 - em[i]) * sigma_f - ed[i] * ndotE / (4 * MY_PI); } @@ -535,10 +530,6 @@ void FixPolarizeBEMGMRES::gmres_solve(double *x, double *r) rho = fabs(g[k]); -#ifdef _POLARIZE_DEBUG - if (comm->me == 0) - error->warning(FLERR, "itr = {}: k = {}, norm(r) = {} norm(b) = {}", itr, k, rho, normb); -#endif if (rho <= rho_tol && rho <= tol_abs) break; } @@ -568,11 +559,6 @@ void FixPolarizeBEMGMRES::gmres_solve(double *x, double *r) rho = sqrt(vec_dot(r, r, n)); -#ifdef _POLARIZE_DEBUG - if (comm->me == 0) - error->warning(FLERR, "itr = {}: norm(r) = {} norm(b) = {}", itr, rho, normb); -#endif - // Barros et al. suggested the condition: norm(r) < EPSILON norm(b) if (rho < tol_rel * normb) break; @@ -644,8 +630,7 @@ void FixPolarizeBEMGMRES::apply_operator(double *w, double *Aw, int /*n*/) Ey += efield_kspace[i][1]; Ez += efield_kspace[i][2]; } - double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / - epsilon[i]; + double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / epsilon[i]; buffer[idx] = em[i] * w[idx] + ed[i] * ndotE / (4 * MY_PI); } @@ -717,7 +702,7 @@ void FixPolarizeBEMGMRES::update_residual(double *w, double *r, int /*n*/) Ez += efield_kspace[i][2]; } double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / - epsilon[i] / (4 * MY_PI); + epsilon[i] / (4 * MY_PI); double sigma_f = q_real[i] / area[i]; buffer[idx] = (1 - em[i]) * sigma_f - em[i] * w[idx] - ed[i] * ndotE; } diff --git a/src/DIELECTRIC/fix_polarize_bem_gmres.h b/src/DIELECTRIC/fix_polarize_bem_gmres.h index 02b9a9c744..950be39868 100644 --- a/src/DIELECTRIC/fix_polarize_bem_gmres.h +++ b/src/DIELECTRIC/fix_polarize_bem_gmres.h @@ -81,7 +81,7 @@ class FixPolarizeBEMGMRES : public Fix { int randomized; // 1 if generating random induced charges, 0 otherwise double ave_charge; // average random charge int seed_charge; - double epsilon0e2q; // convert epsilon0 times efield to unit of charge per area + double epsilon0e2q; // convert epsilon0 times efield to unit of charge per area double *c, *g, *h, *r, *s, *v, *y; // vectors used by the solver double *rhs; // right-hand side vector of the equation Ax = b diff --git a/src/DIELECTRIC/fix_polarize_bem_icc.cpp b/src/DIELECTRIC/fix_polarize_bem_icc.cpp index 8cd0fb4332..02cc96e20f 100644 --- a/src/DIELECTRIC/fix_polarize_bem_icc.cpp +++ b/src/DIELECTRIC/fix_polarize_bem_icc.cpp @@ -49,17 +49,15 @@ using namespace LAMMPS_NS; using namespace FixConst; -using namespace MathConst; - -//#define _POLARIZE_DEBUG +using MathConst::MY_PI; /* ---------------------------------------------------------------------- */ -FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *_lmp, int narg, char **arg) : Fix(_lmp, narg, arg) { if (narg < 5) error->all(FLERR, "Illegal fix polarize/bem/icc command"); - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Fix polarize requires atom style dielectric"); // parse required arguments @@ -147,40 +145,38 @@ void FixPolarizeBEMICC::setup(int /*vflag*/) // check if the pair styles in use are compatible if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else error->all(FLERR, "Pair style not compatible with fix polarize/bem/icc"); // check if kspace is used for force computation if (force->kspace) { - kspaceflag = 1; if (strcmp(force->kspace_style, "pppm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else if (strcmp(force->kspace_style, "msm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else error->all(FLERR, "Kspace style not compatible with fix polarize/bem/icc"); } else { - if (kspaceflag == 1) { // users specified kspace yes error->warning(FLERR, "No Kspace style available for fix polarize/bem/icc"); kspaceflag = 0; @@ -260,7 +256,7 @@ void FixPolarizeBEMICC::compute_induced_charges() // divide (Ex,Ey,Ez) by epsilon[i] here double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / - epsilon[i] / (2 * MY_PI); + epsilon[i] / (2 * MY_PI); double q_free = q_real[i]; double q_bound = 0; q_bound = (1.0 / em[i] - 1) * q_free - (ed[i] / (2 * em[i])) * ndotE * area[i]; @@ -298,7 +294,7 @@ void FixPolarizeBEMICC::compute_induced_charges() // for direct use in force/efield compute double ndotE = epsilon0e2q * (Ex * norm[i][0] + Ey * norm[i][1] + Ez * norm[i][2]) / - (4 * MY_PI) / epsilon[i]; + (4 * MY_PI) / epsilon[i]; double q_bound = q[i] - q_free; q_bound = (1 - omega) * q_bound + omega * ((1.0 / em[i] - 1) * q_free - (ed[i] / em[i]) * ndotE * area[i]); @@ -320,18 +316,11 @@ void FixPolarizeBEMICC::compute_induced_charges() double delta = fabs(qtmp - q_bound); double r = (fabs(qtmp) > 0) ? delta / fabs(qtmp) : 0; if (tol < r) tol = r; - -#ifdef _POLARIZE_DEBUG -//printf("i = %d: q_bound = %f \n", i, q_bound); -#endif } comm->forward_comm(this); MPI_Allreduce(&tol, &rho, 1, MPI_DOUBLE, MPI_MAX, world); -#ifdef _POLARIZE_DEBUG - printf("itr = %d: rho = %f\n", itr, rho); -#endif if (itr > 0 && rho < tol_rel) break; } diff --git a/src/DIELECTRIC/fix_polarize_bem_icc.h b/src/DIELECTRIC/fix_polarize_bem_icc.h index ff59d68059..ade1438496 100644 --- a/src/DIELECTRIC/fix_polarize_bem_icc.h +++ b/src/DIELECTRIC/fix_polarize_bem_icc.h @@ -60,7 +60,7 @@ class FixPolarizeBEMICC : public Fix { int randomized; // 1 if generating random induced charges, 0 otherwise double ave_charge; // average random charge int seed_charge; - double epsilon0e2q; // convert epsilon0 times efield to unit of charge per area + double epsilon0e2q; // convert epsilon0 times efield to unit of charge per area }; } // namespace LAMMPS_NS diff --git a/src/DIELECTRIC/fix_polarize_functional.cpp b/src/DIELECTRIC/fix_polarize_functional.cpp index adf71ed54f..482bc001a9 100644 --- a/src/DIELECTRIC/fix_polarize_functional.cpp +++ b/src/DIELECTRIC/fix_polarize_functional.cpp @@ -61,18 +61,16 @@ using namespace MathSpecial; enum { REAL2SCALED = 0, SCALED2REAL = 1 }; -#define EPSILON 1e-6 - -//#define _POLARIZE_DEBUG +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -FixPolarizeFunctional::FixPolarizeFunctional(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixPolarizeFunctional::FixPolarizeFunctional(LAMMPS *_lmp, int narg, char **arg) : + Fix(_lmp, narg, arg) { if (narg < 4) error->all(FLERR, "Illegal fix polarize/functional command"); - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Fix polarize/functional requires atom style dielectric"); nevery = utils::inumeric(FLERR, arg[3], false, lmp); @@ -291,23 +289,23 @@ void FixPolarizeFunctional::setup(int /*vflag*/) // check if the pair styles in use are compatible if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; - else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; + else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0) + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/long/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0) - efield_pair = (dynamic_cast( force->pair))->efield; + efield_pair = (dynamic_cast(force->pair))->efield; else error->all(FLERR, "Pair style not compatible with fix polarize/functional"); @@ -315,9 +313,9 @@ void FixPolarizeFunctional::setup(int /*vflag*/) kspaceflag = 1; if (strcmp(force->kspace_style, "pppm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else if (strcmp(force->kspace_style, "msm/dielectric") == 0) - efield_kspace = (dynamic_cast( force->kspace))->efield; + efield_kspace = (dynamic_cast(force->kspace))->efield; else error->all(FLERR, "Kspace style not compatible with fix polarize/functional"); @@ -582,21 +580,21 @@ void FixPolarizeFunctional::set_arrays(int i) double FixPolarizeFunctional::memory_usage() { double bytes = 0; - bytes += square(num_induced_charges) * sizeof(double); // inverse_matrix - bytes += square(num_induced_charges) * sizeof(double); // Rww - bytes += square(num_induced_charges) * sizeof(double); // G1ww - bytes += square(num_induced_charges) * sizeof(double); // ndotGww - bytes += square(num_induced_charges) * sizeof(double); // G2ww - bytes += square(num_induced_charges) * sizeof(double); // G3ww - bytes += num_induced_charges * sizeof(double); // qiRqwVector - bytes += num_induced_charges * sizeof(double); // sum2G2wq - bytes += num_induced_charges * sizeof(double); // sum1G2qw - bytes += num_induced_charges * sizeof(double); // sum1G1qw_epsilon - bytes += num_induced_charges * sizeof(double); // sum2ndotGwq_epsilon - bytes += (double) num_ions * num_induced_charges * sizeof(double); // G1qw_real - bytes += nmax * sizeof(int); // induced_charge_idx - bytes += nmax * sizeof(int); // ion_idx - bytes += num_induced_charges * sizeof(double); // induced_charges + bytes += square(num_induced_charges) * sizeof(double); // inverse_matrix + bytes += square(num_induced_charges) * sizeof(double); // Rww + bytes += square(num_induced_charges) * sizeof(double); // G1ww + bytes += square(num_induced_charges) * sizeof(double); // ndotGww + bytes += square(num_induced_charges) * sizeof(double); // G2ww + bytes += square(num_induced_charges) * sizeof(double); // G3ww + bytes += num_induced_charges * sizeof(double); // qiRqwVector + bytes += num_induced_charges * sizeof(double); // sum2G2wq + bytes += num_induced_charges * sizeof(double); // sum1G2qw + bytes += num_induced_charges * sizeof(double); // sum1G1qw_epsilon + bytes += num_induced_charges * sizeof(double); // sum2ndotGwq_epsilon + bytes += (double) num_ions * num_induced_charges * sizeof(double); // G1qw_real + bytes += nmax * sizeof(int); // induced_charge_idx + bytes += nmax * sizeof(int); // ion_idx + bytes += num_induced_charges * sizeof(double); // induced_charges return bytes; } @@ -809,16 +807,6 @@ void FixPolarizeFunctional::calculate_Rww_cutoff() MPI_Allreduce(buffer1[0], Rww[0], num_induced_charges * num_induced_charges, MPI_DOUBLE, MPI_SUM, world); - -#ifdef _POLARIZE_DEBUG - if (comm->me == 0) { - FILE *fp = fopen("Rww-functional.txt", "w"); - for (int i = 0; i < num_induced_charges; i++) - fprintf(fp, "%d %g %g %g\n", i, Rww[i][i], Rww[i][num_induced_charges / 2], - Rww[num_induced_charges / 2][i]); - fclose(fp); - } -#endif } /* ---------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/msm_dielectric.cpp b/src/DIELECTRIC/msm_dielectric.cpp index d3d3da7ed7..3217c8dbad 100644 --- a/src/DIELECTRIC/msm_dielectric.cpp +++ b/src/DIELECTRIC/msm_dielectric.cpp @@ -34,7 +34,7 @@ enum{REVERSE_RHO,REVERSE_AD,REVERSE_AD_PERATOM}; enum{FORWARD_RHO,FORWARD_AD,FORWARD_AD_PERATOM}; /* ---------------------------------------------------------------------- */ -MSMDielectric::MSMDielectric(LAMMPS *lmp) : MSM(lmp) +MSMDielectric::MSMDielectric(LAMMPS *_lmp) : MSM(_lmp) { efield = nullptr; phi = nullptr; diff --git a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp index dc1b2041d7..ba8c665729 100644 --- a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -29,13 +29,13 @@ #include using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairCoulCutDielectric::PairCoulCutDielectric(LAMMPS *lmp) : PairCoulCut(lmp) +PairCoulCutDielectric::PairCoulCutDielectric(LAMMPS *_lmp) : PairCoulCut(_lmp) { efield = nullptr; nmax = 0; @@ -162,7 +162,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag) void PairCoulCutDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair coul/cut/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_coul_long_dielectric.cpp index c8e5424a92..ea705bc29e 100644 --- a/src/DIELECTRIC/pair_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_long_dielectric.cpp @@ -30,8 +30,6 @@ #include using namespace LAMMPS_NS; -using namespace MathConst; - #define EWALD_F 1.12837917 #define EWALD_P 0.3275911 #define A1 0.254829592 @@ -39,10 +37,11 @@ using namespace MathConst; #define A3 1.421413741 #define A4 -1.453152027 #define A5 1.061405429 +using MathConst::MY_PIS; /* ---------------------------------------------------------------------- */ -PairCoulLongDielectric::PairCoulLongDielectric(LAMMPS *lmp) : PairCoulLong(lmp) +PairCoulLongDielectric::PairCoulLongDielectric(LAMMPS *_lmp) : PairCoulLong(_lmp) { efield = nullptr; nmax = 0; @@ -207,7 +206,7 @@ void PairCoulLongDielectric::compute(int eflag, int vflag) void PairCoulLongDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair coul/long/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp index ee95c1c2f5..905bd89550 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -29,13 +29,13 @@ #include using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulCutDielectric::PairLJCutCoulCutDielectric(LAMMPS *lmp) : PairLJCutCoulCut(lmp) +PairLJCutCoulCutDielectric::PairLJCutCoulCutDielectric(LAMMPS *_lmp) : PairLJCutCoulCut(_lmp) { efield = nullptr; epot = nullptr; @@ -190,7 +190,7 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) void PairLJCutCoulCutDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair lj/cut/coul/cut/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp index dd80c5d6ed..814215872b 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp @@ -29,13 +29,13 @@ #include using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulDebyeDielectric::PairLJCutCoulDebyeDielectric(LAMMPS *lmp) : PairLJCutCoulDebye(lmp) +PairLJCutCoulDebyeDielectric::PairLJCutCoulDebyeDielectric(LAMMPS *_lmp) : PairLJCutCoulDebye(_lmp) { efield = nullptr; epot = nullptr; @@ -65,8 +65,8 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag) memory->destroy(efield); memory->destroy(epot); nmax = atom->nmax; - memory->create(efield,nmax,3,"pair:efield"); - memory->create(epot,nmax,"pair:epot"); + memory->create(efield, nmax, 3, "pair:efield"); + memory->create(epot, nmax, "pair:epot"); } evdwl = ecoul = 0.0; @@ -193,7 +193,7 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag) void PairLJCutCoulDebyeDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair lj/cut/coul/debye/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h index 9bd4086e2a..dbe626e669 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h +++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h @@ -40,7 +40,7 @@ class PairLJCutCoulDebyeDielectric : public PairLJCutCoulDebye { int nmax; }; -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp index 823c7de50f..a206e66208 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -40,11 +40,11 @@ using namespace MathConst; #define A4 -1.453152027 #define A5 1.061405429 -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulLongDielectric::PairLJCutCoulLongDielectric(LAMMPS *lmp) : PairLJCutCoulLong(lmp) +PairLJCutCoulLongDielectric::PairLJCutCoulLongDielectric(LAMMPS *_lmp) : PairLJCutCoulLong(_lmp) { respa_enable = 0; cut_respa = nullptr; @@ -244,7 +244,7 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag) void PairLJCutCoulLongDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair lj/cut/coul/long/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp index 4dffd81aa3..9a3c74433f 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp @@ -31,13 +31,13 @@ #include using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulMSMDielectric::PairLJCutCoulMSMDielectric(LAMMPS *lmp) : PairLJCutCoulLong(lmp) +PairLJCutCoulMSMDielectric::PairLJCutCoulMSMDielectric(LAMMPS *_lmp) : PairLJCutCoulLong(_lmp) { ewaldflag = pppmflag = 0; msmflag = 1; @@ -352,7 +352,7 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, do void PairLJCutCoulMSMDielectric::init_style() { - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair lj/cut/coul/msm/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index 0f72238e63..e88d3b99cd 100644 --- a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -44,7 +44,7 @@ using namespace MathExtra; /* ---------------------------------------------------------------------- */ -PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp) : PairLJLongCoulLong(lmp) +PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *_lmp) : PairLJLongCoulLong(_lmp) { respa_enable = 0; cut_respa = nullptr; @@ -71,7 +71,7 @@ void PairLJLongCoulLongDielectric::init_style() { PairLJLongCoulLong::init_style(); - avec = dynamic_cast( atom->style_match("dielectric")); + avec = dynamic_cast(atom->style_match("dielectric")); if (!avec) error->all(FLERR, "Pair lj/long/coul/long/dielectric requires atom style dielectric"); neighbor->add_request(this, NeighConst::REQ_FULL); diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp index 25e5936cdb..7f32a0a3f3 100644 --- a/src/DIELECTRIC/pppm_dielectric.cpp +++ b/src/DIELECTRIC/pppm_dielectric.cpp @@ -50,7 +50,7 @@ enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM}; /* ---------------------------------------------------------------------- */ -PPPMDielectric::PPPMDielectric(LAMMPS *lmp) : PPPM(lmp) +PPPMDielectric::PPPMDielectric(LAMMPS *_lmp) : PPPM(_lmp) { group_group_enable = 0; diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index 36847a87ac..e525ba7384 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -58,7 +58,7 @@ enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM, /* ---------------------------------------------------------------------- */ -PPPMDispDielectric::PPPMDispDielectric(LAMMPS *lmp) : PPPMDisp(lmp) +PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp) { dipoleflag = 0; // turned off for now, until dipole works group_group_enable = 0; diff --git a/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp index 04f2730f13..e49d033b4f 100644 --- a/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp +++ b/src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp @@ -28,14 +28,14 @@ #include "omp_compat.h" using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulCutDielectricOMP::PairLJCutCoulCutDielectricOMP(LAMMPS *lmp) : - PairLJCutCoulCutDielectric(lmp), ThrOMP(lmp, THR_PAIR) +PairLJCutCoulCutDielectricOMP::PairLJCutCoulCutDielectricOMP(LAMMPS *_lmp) : + PairLJCutCoulCutDielectric(_lmp), ThrOMP(_lmp, THR_PAIR) { } diff --git a/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp index 1b33ee8b5f..032a2b4c3c 100644 --- a/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp +++ b/src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp @@ -28,14 +28,14 @@ #include "omp_compat.h" using namespace LAMMPS_NS; -using namespace MathConst; +using MathConst::MY_PIS; -#define EPSILON 1e-6 +static constexpr double EPSILON = 1.0e-6; /* ---------------------------------------------------------------------- */ -PairLJCutCoulDebyeDielectricOMP::PairLJCutCoulDebyeDielectricOMP(LAMMPS *lmp) : - PairLJCutCoulDebyeDielectric(lmp), ThrOMP(lmp, THR_PAIR) +PairLJCutCoulDebyeDielectricOMP::PairLJCutCoulDebyeDielectricOMP(LAMMPS *_lmp) : + PairLJCutCoulDebyeDielectric(_lmp), ThrOMP(_lmp, THR_PAIR) { } diff --git a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp index 5c773891f8..4040b7d2f0 100644 --- a/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp +++ b/src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp @@ -40,8 +40,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairLJCutCoulLongDielectricOMP::PairLJCutCoulLongDielectricOMP(LAMMPS *lmp) : - PairLJCutCoulLongDielectric(lmp), ThrOMP(lmp, THR_PAIR) +PairLJCutCoulLongDielectricOMP::PairLJCutCoulLongDielectricOMP(LAMMPS *_lmp) : + PairLJCutCoulLongDielectric(_lmp), ThrOMP(_lmp, THR_PAIR) { }