diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index 5dc739d5e8..b25ce479cd 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -405,6 +405,9 @@ void FixQEqReaxFF::init() if (strcmp(update->unit_style,"real") != 0) error->all(FLERR,"Must use unit_style real with fix {} and external fields", style); + if (efield->varflag == FixEfield::ATOM && efield->pstyle != FixEfield::ATOM) + error->all(FLERR,"Atom-style external electric field requires atom-style " + "potential variable when used with fix {}", style); if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || ((fabs(efield->ey) > SMALL) && domain->yperiodic) || ((fabs(efield->ez) > SMALL) && domain->zperiodic)) @@ -1125,19 +1128,11 @@ void FixQEqReaxFF::get_chi_field() chi_field[i] = factor*(ex*unwrap[0] + ey*unwrap[1] + ez*unwrap[2]); } } - } else { + } else { // must use atom-style potential from FixEfield for (int i = 0; i < nlocal; i++) { if (mask[i] & efgroupbit) { if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - domain->unmap(x[i],image[i],unwrap); - double edisp = 0; // accumulate E dot displacement - edisp += unwrap[0]*( - (efield->xstyle == FixEfield::ATOM) ? qe2f*efield->efield[i][0] : ex); - edisp += unwrap[1]*( - (efield->ystyle == FixEfield::ATOM) ? qe2f*efield->efield[i][1] : ey); - edisp += unwrap[2]*( - (efield->zstyle == FixEfield::ATOM) ? qe2f*efield->efield[i][2] : ez); - chi_field[i] = factor*edisp; + chi_field[i] = -efield->efield[i][3]; } } } diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index f1adc722ae..f335d8a765 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -41,7 +41,7 @@ using namespace FixConst; FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), estr(nullptr), - idregion(nullptr), region(nullptr), efield(nullptr) + pstr(nullptr), idregion(nullptr), region(nullptr), efield(nullptr) { if (narg < 6) utils::missing_cmd_args(FLERR, std::string("fix ") + style, error); @@ -100,6 +100,14 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR, "Unsupported argument for fix {} energy command: {}", style, arg[iarg]); iarg += 2; + } else if (strcmp(arg[iarg], "potential") == 0) { + if (iarg + 2 > narg) + utils::missing_cmd_args(FLERR, std::string("fix ") + style + "potential", error); + if (utils::strmatch(arg[iarg + 1], "^v_")) { + pstr = utils::strdup(arg[iarg + 1] + 2); + } else + error->all(FLERR, "Unsupported argument for fix {} energy command: {}", style, arg[iarg]); + iarg += 2; } else { error->all(FLERR, "Unknown keyword for fix {} command: {}", style, arg[iarg]); } @@ -122,6 +130,7 @@ FixEfield::~FixEfield() delete[] ystr; delete[] zstr; delete[] estr; + delete[] pstr; delete[] idregion; memory->destroy(efield); } @@ -157,43 +166,55 @@ void FixEfield::init() if (xstr) { xvar = input->variable->find(xstr); - if (xvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", xstr, style); + if (xvar < 0) error->all(FLERR, "Variable {} for x-field in fix {} does not exist", xstr, style); if (input->variable->equalstyle(xvar)) xstyle = EQUAL; else if (input->variable->atomstyle(xvar)) xstyle = ATOM; else - error->all(FLERR, "Variable {} for fix {} is invalid style", xstr, style); + error->all(FLERR, "Variable {} for x-field in fix {} is invalid style", xstr, style); } if (ystr) { yvar = input->variable->find(ystr); - if (yvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", ystr, style); + if (yvar < 0) error->all(FLERR, "Variable {} for y-field in fix {} does not exist", ystr, style); if (input->variable->equalstyle(yvar)) ystyle = EQUAL; else if (input->variable->atomstyle(yvar)) ystyle = ATOM; else - error->all(FLERR, "Variable {} for fix {} is invalid style", ystr, style); + error->all(FLERR, "Variable {} for y-field in fix {} is invalid style", ystr, style); } if (zstr) { zvar = input->variable->find(zstr); - if (zvar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", zstr, style); + if (zvar < 0) error->all(FLERR, "Variable {} for z-field in fix {} does not exist", zstr, style); if (input->variable->equalstyle(zvar)) zstyle = EQUAL; else if (input->variable->atomstyle(zvar)) zstyle = ATOM; else - error->all(FLERR, "Variable {} for fix {} is invalid style", zstr, style); + error->all(FLERR, "Variable {} for z-field in fix {} is invalid style", zstr, style); } if (estr) { evar = input->variable->find(estr); - if (evar < 0) error->all(FLERR, "Variable {} for fix {} does not exist", estr, style); + if (evar < 0) error->all(FLERR, "Variable {} for energy in fix {} does not exist", estr, style); if (input->variable->atomstyle(evar)) estyle = ATOM; else - error->all(FLERR, "Variable {} for fix {} is invalid style", estr, style); + error->all(FLERR, "Variable {} for energy in fix {} must be atom-style", estr, style); } else estyle = NONE; + if (pstr) { + pvar = input->variable->find(pstr); + if (pvar < 0) error->all(FLERR, "Variable {} for potential in fix {} does not exist", pstr, style); + if (input->variable->atomstyle(pvar)) + pstyle = ATOM; + else + error->all(FLERR, "Variable {} for potential in fix {} must be atom-style", pstr, style); + if (estyle != NONE) + error->warning(FLERR, "fix {} will ignore variable {} for energy " + "because atom-style potential has been specified", estr, style); + } else + pstyle = NONE; // set index and check validity of region @@ -376,7 +397,8 @@ void FixEfield::post_force(int vflag) } f[i][2] += fz; fsum[3] += fz; - if (estyle == ATOM) fsum[0] += efield[i][3]; + if (pstyle == ATOM) fsum[0] += qe2f * q[i] * efield[i][3]; + else if (estyle == ATOM) fsum[0] += efield[i][3]; } } @@ -476,7 +498,8 @@ void FixEfield::update_efield_variables() } else if (zstyle == ATOM) { input->variable->compute_atom(zvar, igroup, &efield[0][2], 4, 0); } - if (estyle == ATOM) input->variable->compute_atom(evar, igroup, &efield[0][3], 4, 0); + if (pstyle == ATOM) input->variable->compute_atom(pvar, igroup, &efield[0][3], 4, 0); + else if (estyle == ATOM) input->variable->compute_atom(evar, igroup, &efield[0][3], 4, 0); modify->addstep_compute(update->ntimestep + 1); } diff --git a/src/fix_efield.h b/src/fix_efield.h index 2bab26e140..72fd204898 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -46,10 +46,11 @@ class FixEfield : public Fix { protected: double ex, ey, ez; int varflag; - char *xstr, *ystr, *zstr, *estr; + char *xstr, *ystr, *zstr, *estr, *pstr; char *idregion; class Region *region; - int xvar, yvar, zvar, evar, xstyle, ystyle, zstyle, estyle; + int xvar, yvar, zvar, xstyle, ystyle, zstyle; + int evar, pvar, estyle, pstyle; int ilevel_respa; double qe2f; int qflag, muflag;