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

programming style updates and clang-format applied

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-14 16:47:10 -04:00
parent 81e203b5fa
commit 8b31edb102
21 changed files with 157 additions and 196 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/DIELECTRIC/atom_vec_dielectric.cpp
View File

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

84
src/DIELECTRIC/compute_efield_atom.cpp
View File

@ -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<PairLJCutCoulLongDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/msm/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"coul/long/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulLongDielectric*>(force->pair))->efield;
else if (strcmp(force->pair_style,"coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulCutDielectric*>(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<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulCutDielectric *>(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<PPPMDielectric*>(force->kspace))->efield;
else if (strcmp(force->kspace_style,"msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric*>(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<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>(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;
}

59
src/DIELECTRIC/fix_polarize_bem_gmres.cpp
View File

@ -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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulCutDielectric *>(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<PPPMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<MSMDielectric *>(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;
}

2
src/DIELECTRIC/fix_polarize_bem_gmres.h
View File

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

47
src/DIELECTRIC/fix_polarize_bem_icc.cpp
View File

@ -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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulCutDielectric *>(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<PPPMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<MSMDielectric *>(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;
}

2
src/DIELECTRIC/fix_polarize_bem_icc.h
View File

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

76
src/DIELECTRIC/fix_polarize_functional.cpp
View File

@ -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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/long/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/msm/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulMSMDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/cut/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
else if (strcmp(force->pair_style,"lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairLJCutCoulCutDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "lj/cut/coul/debye/dielectric/omp") == 0)
efield_pair = (dynamic_cast<PairLJCutCoulDebyeDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/long/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulLongDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulLongDielectric *>(force->pair))->efield;
else if (strcmp(force->pair_style, "coul/cut/dielectric") == 0)
efield_pair = (dynamic_cast<PairCoulCutDielectric *>( force->pair))->efield;
efield_pair = (dynamic_cast<PairCoulCutDielectric *>(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<PPPMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<PPPMDielectric *>(force->kspace))->efield;
else if (strcmp(force->kspace_style, "msm/dielectric") == 0)
efield_kspace = (dynamic_cast<MSMDielectric *>( force->kspace))->efield;
efield_kspace = (dynamic_cast<MSMDielectric *>(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
}
/* ---------------------------------------------------------------------- */

2
src/DIELECTRIC/msm_dielectric.cpp
View File

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

8
src/DIELECTRIC/pair_coul_cut_dielectric.cpp
View File

@ -29,13 +29,13 @@
#include <cmath>
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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(atom->style_match("dielectric"));
if (!avec) error->all(FLERR, "Pair coul/cut/dielectric requires atom style dielectric");
neighbor->add_request(this, NeighConst::REQ_FULL);

7
src/DIELECTRIC/pair_coul_long_dielectric.cpp
View File

@ -30,8 +30,6 @@
#include <cmath>
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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(atom->style_match("dielectric"));
if (!avec) error->all(FLERR, "Pair coul/long/dielectric requires atom style dielectric");
neighbor->add_request(this, NeighConst::REQ_FULL);

8
src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
View File

@ -29,13 +29,13 @@
#include <cmath>
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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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);

12
src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
View File

@ -29,13 +29,13 @@
#include <cmath>
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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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);

2
src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h
View File

@ -40,7 +40,7 @@ class PairLJCutCoulDebyeDielectric : public PairLJCutCoulDebye {
int nmax;
};
}
} // namespace LAMMPS_NS
#endif
#endif

6
src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
View File

@ -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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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);

8
src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
View File

@ -31,13 +31,13 @@
#include <cstring>
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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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);

4
src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
View File

@ -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<AtomVecDielectric *>( atom->style_match("dielectric"));
avec = dynamic_cast<AtomVecDielectric *>(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);

2
src/DIELECTRIC/pppm_dielectric.cpp
View File

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

2
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

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

8
src/OPENMP/pair_lj_cut_coul_cut_dielectric_omp.cpp
View File

@ -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)
{
}

8
src/OPENMP/pair_lj_cut_coul_debye_dielectric_omp.cpp
View File

@ -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)
{
}

4
src/OPENMP/pair_lj_cut_coul_long_dielectric_omp.cpp
View File

@ -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)
{
}