tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

use utils::numeric() instead of atof and improve error messages in QEQ package

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-12-07 13:34:09 -05:00
parent 92e6c6ea9d
commit 739537930d
6 changed files with 134 additions and 143 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/QEQ/fix_qeq.cpp
View File

@ -60,7 +60,7 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
b_t(nullptr), p(nullptr), q(nullptr), r(nullptr), d(nullptr),
qf(nullptr), q1(nullptr), q2(nullptr), qv(nullptr)
{
if (narg < 8) error->all(FLERR,"Illegal fix qeq command");
if (narg < 8) utils::missing_cmd_args(FLERR, "fix " + std::string(style), error);
scalar_flag = 1;
extscalar = 0;

82
src/QEQ/fix_qeq_dynamic.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
@ -34,8 +33,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg)
FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg)
{
qdamp = 0.10;
qstep = 0.02;
@ -43,19 +41,20 @@ FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) :
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"qdamp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/dynamic command");
qdamp = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "qdamp") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/dynamic qdamp", error);
qdamp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"qstep") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/dynamic command");
qstep = atof(arg[iarg+1]);
} else if (strcmp(arg[iarg], "qstep") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/dynamic qstep", error);
qstep = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"warn") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/dynamic command");
maxwarn = utils::logical(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "warn") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/dynamic warn", error);
maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/dynamic command");
} else
error->all(FLERR, "Unknown fix qeq/dynamic keyword: {}", arg[iarg]);
}
}
@ -69,16 +68,17 @@ void FixQEqDynamic::init()
if (tolerance < 1e-4)
if (comm->me == 0)
error->warning(FLERR,"Fix qeq/dynamic tolerance may be too small for damped dynamics");
error->warning(FLERR, "Fix qeq/dynamic tolerance {} may be too small for damped dynamics",
tolerance);
}
/* ---------------------------------------------------------------------- */
void FixQEqDynamic::pre_force(int /*vflag*/)
{
int i,ii,iloop,inum,*ilist;
double qmass,dtq2;
double enegchkall,enegmaxall;
int i, ii, iloop, inum, *ilist;
double qmass, dtq2;
double enegchkall, enegmaxall;
double *q = atom->q;
int *mask = atom->mask;
@ -94,20 +94,20 @@ void FixQEqDynamic::pre_force(int /*vflag*/)
inum = list->inum;
ilist = list->ilist;
qmass = 0.016;
dtq2 = 0.5*qstep*qstep/qmass;
qmass = 0.016;
dtq2 = 0.5 * qstep * qstep / qmass;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
q1[i] = q2[i] = qf[i] = 0.0;
}
for (iloop = 0; iloop < maxiter; iloop ++) {
for (iloop = 0; iloop < maxiter; iloop++) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
q1[i] += qf[i]*dtq2 - qdamp*q1[i];
q[i] += q1[i];
q1[i] += qf[i] * dtq2 - qdamp * q1[i];
q[i] += q1[i];
}
}
@ -119,34 +119,32 @@ void FixQEqDynamic::pre_force(int /*vflag*/)
enegchk = enegmax = 0.0;
for (ii = 0; ii < inum ; ii++) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
q2[i] = enegtot-qf[i];
enegmax = MAX(enegmax,fabs(q2[i]));
q2[i] = enegtot - qf[i];
enegmax = MAX(enegmax, fabs(q2[i]));
enegchk += fabs(q2[i]);
qf[i] = q2[i];
}
}
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
enegchk = enegchkall/ngroup;
MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(&enegchk, &enegchkall, 1, MPI_DOUBLE, MPI_SUM, world);
enegchk = enegchkall / ngroup;
MPI_Allreduce(&enegmax, &enegmaxall, 1, MPI_DOUBLE, MPI_MAX, world);
enegmax = enegmaxall;
if ((enegchk <= tolerance) && (enegmax <= 100.0*tolerance)) break;
if ((enegchk <= tolerance) && (enegmax <= 100.0 * tolerance)) break;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
q1[i] += qf[i]*dtq2 - qdamp*q1[i];
if (mask[i] & groupbit) q1[i] += qf[i] * dtq2 - qdamp * q1[i];
}
}
matvecs = iloop;
if ((comm->me == 0) && maxwarn && (iloop >= maxiter))
error->warning(FLERR,"Charges did not converge at step {}: {}",
update->ntimestep,enegchk);
error->warning(FLERR, "Charges did not converge at step {}: {}", update->ntimestep, enegchk);
if (force->kspace) force->kspace->qsum_qsq();
}
@ -172,8 +170,7 @@ double FixQEqDynamic::compute_eneg()
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
qf[i] = 0.0;
if (mask[i] & groupbit) qf[i] = 0.0;
}
// communicating charge force to all nodes, first forward then reverse
@ -198,12 +195,12 @@ double FixQEqDynamic::compute_eneg()
delr[0] = x[i][0] - x[j][0];
delr[1] = x[i][1] - x[j][1];
delr[2] = x[i][2] - x[j][2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
rsq = delr[0] * delr[0] + delr[1] * delr[1] + delr[2] * delr[2];
if (rsq > cutoff_sq) continue;
r = sqrt(rsq);
rinv = 1.0/r;
rinv = 1.0 / r;
qf[i] += q[j] * rinv;
qf[j] += q[i] * rinv;
}
@ -218,20 +215,17 @@ double FixQEqDynamic::compute_eneg()
eneg = enegtot = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
eneg += qf[i];
if (mask[i] & groupbit) eneg += qf[i];
}
MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&eneg, &enegtot, 1, MPI_DOUBLE, MPI_SUM, world);
return enegtot;
}
/* ---------------------------------------------------------------------- */
int FixQEqDynamic::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int FixQEqDynamic::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int m=0;
int m = 0;
if (pack_flag == 1)
for (m = 0; m < n; m++) buf[m] = atom->q[list[m]];

98
src/QEQ/fix_qeq_fire.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
@ -45,7 +44,7 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixQEqFire::FixQEqFire(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg), comb(nullptr), comb3(nullptr)
FixQEq(lmp, narg, arg), comb(nullptr), comb3(nullptr)
{
qdamp = 0.20;
qstep = 0.20;
@ -53,19 +52,20 @@ FixQEqFire::FixQEqFire(LAMMPS *lmp, int narg, char **arg) :
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"qdamp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/fire command");
qdamp = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "qdamp") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/fire qdamp", error);
qdamp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"qstep") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/fire command");
qstep = atof(arg[iarg+1]);
} else if (strcmp(arg[iarg], "qstep") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/fire qstep", error);
qstep = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"warn") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/fire command");
maxwarn = utils::logical(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "warn") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/fire warn", error);
maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/fire command");
} else
error->all(FLERR, "Unknown fix qeq/fire keyword: {}", arg[iarg]);
}
}
@ -79,10 +79,11 @@ void FixQEqFire::init()
if (tolerance < 1e-4)
if (comm->me == 0)
error->warning(FLERR,"Fix qeq/fire tolerance may be too small for damped fires");
error->warning(FLERR, "Fix qeq/fire tolerance {} may be too small for damped fires",
tolerance);
comb3 = dynamic_cast<PairComb3 *>(force->pair_match("^comb3",0));
if (!comb3) comb = dynamic_cast<PairComb *>(force->pair_match("^comb",0));
comb3 = dynamic_cast<PairComb3 *>(force->pair_match("^comb3", 0));
if (!comb3) comb = dynamic_cast<PairComb *>(force->pair_match("^comb", 0));
}
/* ---------------------------------------------------------------------- */
@ -90,13 +91,13 @@ void FixQEqFire::init()
void FixQEqFire::pre_force(int /*vflag*/)
{
int inum, *ilist;
int i,ii,iloop;
int i, ii, iloop;
double *q = atom->q;
double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
double scale1,scale2;
double dtvone,dtv;
double enegtot,enegchk;
double vmax, vdotf, vdotfall, vdotv, vdotvall, fdotf, fdotfall;
double scale1, scale2;
double dtvone, dtv;
double enegtot, enegchk;
double alpha = qdamp;
double dt, dtmax;
double enegchkall;
@ -118,15 +119,15 @@ void FixQEqFire::pre_force(int /*vflag*/)
dt = qstep;
dtmax = TMAX * dt;
for (iloop = 0; iloop < maxiter; iloop ++) {
for (iloop = 0; iloop < maxiter; iloop++) {
pack_flag = 1;
comm->forward_comm(this);
if (comb) {
comb->yasu_char(qf,igroup);
comb->yasu_char(qf, igroup);
enegtot = comb->enegtot / ngroup;
} else if (comb3) {
comb3->combqeq(qf,igroup);
comb3->combqeq(qf, igroup);
enegtot = comb3->enegtot / ngroup;
} else {
enegtot = compute_eneg();
@ -135,7 +136,7 @@ void FixQEqFire::pre_force(int /*vflag*/)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qf[i] -= enegtot; // Enforce adiabatic
qf[i] -= enegtot; // Enforce adiabatic
}
// FIRE minimization algorithm
@ -143,30 +144,32 @@ void FixQEqFire::pre_force(int /*vflag*/)
vdotf = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
vdotf += (qv[i]*qf[i]);
vdotf += (qv[i] * qf[i]);
}
MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&vdotf, &vdotfall, 1, MPI_DOUBLE, MPI_SUM, world);
if (vdotfall > 0.0) {
vdotv = fdotf = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
vdotv += qv[i]*qv[i];
fdotf += qf[i]*qf[i];
vdotv += qv[i] * qv[i];
fdotf += qf[i] * qf[i];
}
MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&vdotv, &vdotvall, 1, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(&fdotf, &fdotfall, 1, MPI_DOUBLE, MPI_SUM, world);
scale1 = 1.0 - alpha;
if (fdotfall == 0.0) scale2 = 0.0;
else scale2 = alpha * sqrt(vdotvall/fdotfall);
if (fdotfall == 0.0)
scale2 = 0.0;
else
scale2 = alpha * sqrt(vdotvall / fdotfall);
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qv[i] = scale1*qv[i] + scale2*qf[i];
qv[i] = scale1 * qv[i] + scale2 * qf[i];
}
if (ntimestep - last_negative > DELAYSTEP) {
dt = MIN(dt*DT_GROW,dtmax);
dt = MIN(dt * DT_GROW, dtmax);
alpha *= ALPHA_SHRINK;
}
} else {
@ -184,10 +187,10 @@ void FixQEqFire::pre_force(int /*vflag*/)
double dmax = 0.1;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
vmax = MAX(fabs(qv[i]),0);
if (dtvone*vmax > dmax) dtvone = dmax/vmax;
vmax = MAX(fabs(qv[i]), 0);
if (dtvone * vmax > dmax) dtvone = dmax / vmax;
}
MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
MPI_Allreduce(&dtvone, &dtv, 1, MPI_DOUBLE, MPI_MIN, world);
//dtv = dt;
// Euler integration step
@ -198,7 +201,7 @@ void FixQEqFire::pre_force(int /*vflag*/)
qv[i] += dtv * qf[i];
enegchk += fabs(qf[i]);
}
MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&enegchk, &enegchkall, 1, MPI_DOUBLE, MPI_SUM, world);
enegchk = enegchkall / ngroup;
if (enegchk < tolerance) break;
@ -206,8 +209,7 @@ void FixQEqFire::pre_force(int /*vflag*/)
matvecs = iloop;
if ((comm->me == 0) && maxwarn && (iloop >= maxiter))
error->warning(FLERR,"Charges did not converge at step {}: {}",
update->ntimestep,enegchk);
error->warning(FLERR, "Charges did not converge at step {}: {}", update->ntimestep, enegchk);
if (force->kspace) force->kspace->qsum_qsq();
}
@ -233,8 +235,7 @@ double FixQEqFire::compute_eneg()
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
qf[i] = 0.0;
if (mask[i] & groupbit) qf[i] = 0.0;
}
// communicating charge force to all nodes, first forward then reverse
@ -259,12 +260,12 @@ double FixQEqFire::compute_eneg()
delr[0] = x[i][0] - x[j][0];
delr[1] = x[i][1] - x[j][1];
delr[2] = x[i][2] - x[j][2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
rsq = delr[0] * delr[0] + delr[1] * delr[1] + delr[2] * delr[2];
if (rsq > cutoff_sq) continue;
r = sqrt(rsq);
rinv = 1.0/r;
rinv = 1.0 / r;
qf[i] += q[j] * rinv;
qf[j] += q[i] * rinv;
}
@ -279,18 +280,15 @@ double FixQEqFire::compute_eneg()
eneg = enegtot = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
eneg += qf[i];
if (mask[i] & groupbit) eneg += qf[i];
}
MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&eneg, &enegtot, 1, MPI_DOUBLE, MPI_SUM, world);
return enegtot;
}
/* ---------------------------------------------------------------------- */
int FixQEqFire::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int FixQEqFire::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int m = 0;

19
src/QEQ/fix_qeq_point.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
@ -35,13 +34,15 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqPoint::FixQEqPoint(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg) {
FixQEqPoint::FixQEqPoint(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg)
{
if (narg == 10) {
if (strcmp(arg[8],"warn") == 0) {
maxwarn = utils::logical(FLERR,arg[9],false,lmp);
} else error->all(FLERR,"Illegal fix qeq/point command");
} else if (narg > 8) error->all(FLERR,"Illegal fix qeq/point command");
if (strcmp(arg[8], "warn") == 0) {
maxwarn = utils::logical(FLERR, arg[9], false, lmp);
} else
error->all(FLERR, "Illegal fix qeq/point command");
} else if (narg > 8)
error->all(FLERR, "Illegal fix qeq/point command");
}
/* ---------------------------------------------------------------------- */
@ -53,9 +54,11 @@ void FixQEqPoint::init()
neighbor->add_request(this, NeighConst::REQ_FULL);
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding");
memory->create(shld, ntypes + 1, ntypes + 1, "qeq:shielding");
}
// clang-format off
/* ---------------------------------------------------------------------- */
void FixQEqPoint::pre_force(int /*vflag*/)

32
src/QEQ/fix_qeq_shielded.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
@ -36,13 +35,15 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg) {
FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg)
{
if (narg == 10) {
if (strcmp(arg[8],"warn") == 0) {
maxwarn = utils::logical(FLERR,arg[9],false,lmp);
} else error->all(FLERR,"Illegal fix qeq/shielded command");
} else if (narg > 8) error->all(FLERR,"Illegal fix qeq/shielded command");
if (strcmp(arg[8], "warn") == 0) {
maxwarn = utils::logical(FLERR, arg[9], false, lmp);
} else
error->all(FLERR, "Illegal fix qeq/shielded command");
} else if (narg > 8)
error->all(FLERR, "Illegal fix qeq/shielded command");
if (reax_flag) extract_reax();
}
@ -55,14 +56,13 @@ void FixQEqShielded::init()
neighbor->add_request(this, NeighConst::REQ_FULL);
int ntypes = atom->ntypes;
memory->create(shld,ntypes+1,ntypes+1,"qeq:shielding");
memory->create(shld, ntypes + 1, ntypes + 1, "qeq:shielding");
init_shielding();
int i;
for (i = 1; i <= ntypes; i++) {
if (gamma[i] == 0.0)
error->all(FLERR,"Invalid param file for fix qeq/shielded");
if (gamma[i] == 0.0) error->all(FLERR, "Invalid param file for fix qeq/shielded");
}
}
@ -70,17 +70,17 @@ void FixQEqShielded::init()
void FixQEqShielded::extract_reax()
{
Pair *pair = force->pair_match("^reax..",0);
if (pair == nullptr) error->all(FLERR,"No pair reaxff for fix qeq/shielded");
Pair *pair = force->pair_match("^reax..", 0);
if (pair == nullptr) error->all(FLERR, "No pair reaxff for fix qeq/shielded");
int tmp;
chi = (double *) pair->extract("chi",tmp);
eta = (double *) pair->extract("eta",tmp);
gamma = (double *) pair->extract("gamma",tmp);
chi = (double *) pair->extract("chi", tmp);
eta = (double *) pair->extract("eta", tmp);
gamma = (double *) pair->extract("gamma", tmp);
if (chi == nullptr || eta == nullptr || gamma == nullptr)
error->all(FLERR, "Fix qeq/shielded could not extract params from pair reaxff");
}
// clang-format off
/* ---------------------------------------------------------------------- */
void FixQEqShielded::init_shielding()

44
src/QEQ/fix_qeq_slater.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
@ -38,23 +37,23 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) :
FixQEq(lmp, narg, arg)
FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg)
{
alpha = 0.20;
// optional arg
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"alpha") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/slater command");
alpha = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "alpha") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/slater alpha", error);
alpha = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"warn") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/slater command");
maxwarn = utils::logical(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "warn") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/slater warn", error);
maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/slater command");
} else
error->all(FLERR, "Unknown fix qeq/slater keyword: {}", arg[iarg]);
}
if (streitz_flag) extract_streitz();
@ -70,8 +69,7 @@ void FixQEqSlater::init()
int ntypes = atom->ntypes;
for (int i = 1; i <= ntypes; i++) {
if (zeta[i] == 0.0)
error->all(FLERR,"Invalid param file for fix qeq/slater");
if (zeta[i] == 0.0) error->all(FLERR, "Invalid parameter file values for fix qeq/slater");
}
}
@ -79,21 +77,19 @@ void FixQEqSlater::init()
void FixQEqSlater::extract_streitz()
{
Pair *pair = force->pair_match("coul/streitz",1);
if (pair == nullptr) error->all(FLERR,"No pair coul/streitz for fix qeq/slater");
Pair *pair = force->pair_match("coul/streitz", 1);
if (pair == nullptr) error->all(FLERR, "No pair style coul/streitz for fix qeq/slater");
int tmp;
chi = (double *) pair->extract("chi",tmp);
eta = (double *) pair->extract("eta",tmp);
gamma = (double *) pair->extract("gamma",tmp);
zeta = (double *) pair->extract("zeta",tmp);
zcore = (double *) pair->extract("zcore",tmp);
if (chi == nullptr || eta == nullptr || gamma == nullptr
|| zeta == nullptr || zcore == nullptr)
error->all(FLERR,
"Fix qeq/slater could not extract params from pair coul/streitz");
chi = (double *) pair->extract("chi", tmp);
eta = (double *) pair->extract("eta", tmp);
gamma = (double *) pair->extract("gamma", tmp);
zeta = (double *) pair->extract("zeta", tmp);
zcore = (double *) pair->extract("zcore", tmp);
if (chi == nullptr || eta == nullptr || gamma == nullptr || zeta == nullptr || zcore == nullptr)
error->all(FLERR, "Fix qeq/slater could not extract parameters from pair coul/streitz");
}
// clang-format off
/* ---------------------------------------------------------------------- */
void FixQEqSlater::pre_force(int /*vflag*/)