diff --git a/doc/src/fix_efield.rst b/doc/src/fix_efield.rst index e72510e4da..4cf92cf946 100644 --- a/doc/src/fix_efield.rst +++ b/doc/src/fix_efield.rst @@ -19,7 +19,7 @@ Syntax * ex,ey,ez = E-field component values (electric field units) * any of ex,ey,ez can be a variable (see below) * zero or more keyword/value pairs may be appended to args -* keyword = *region* or *energy* +* keyword = *region* or *energy* or *potential* .. parsed-literal:: @@ -27,6 +27,8 @@ Syntax region-ID = ID of region atoms must be in to have added force *energy* value = v_name v_name = variable with name that calculates the potential energy of each atom in the added E-field + *potential* value = v_name + v_name = variable with name that calculates the electric potential of each atom in the added E-field Examples """""""" @@ -112,7 +114,8 @@ one or more variables, and if you are performing dynamics via the :doc:`run ` command. If the keyword is not used, LAMMPS will set the energy to 0.0, which is typically fine for dynamics. -The *energy* keyword is required if the added force is defined with +The *energy* keyword (or *potential* keyword, described below) +is required if the added force is defined with one or more variables, and you are performing energy minimization via the "minimize" command for charged particles. It is not required for point-dipoles, but a warning is issued since the minimizer in LAMMPS @@ -122,7 +125,7 @@ minimize the orientation of dipoles in an applied electric field. The *energy* keyword specifies the name of an atom-style :doc:`variable ` which is used to compute the energy of each atom as function of its position. Like variables used for *ex*, -*ey*, *ez*, the energy variable is specified as v_name, where name +*ey*, *ez*, the energy variable is specified as "v_name", where "name" is the variable name. Note that when the *energy* keyword is used during an energy @@ -133,6 +136,27 @@ due to the electric field were a spring-like F = kx, then the energy formula should be E = -0.5kx\^2. If you don't do this correctly, the minimization will not converge properly. +.. versionadded:: TBD + +The *potential* keyword can be used as an alternative to the *energy* keyword +to specify the name of an atom-style variable, which is used to compute the +added electric potential to each atom as a function of its position. The +variable should have units of electric field multiplied by distance (that is, +in `units real`, the potential should be in volts). As with the *energy* +keyword, the variable name is specified as "v_name". The energy added by this +fix is then calculated as the electric potential multiplied by charge. + +The *potential* keyword is mainly intended for correct charge +equilibration in simulations with :doc:`fix qeq/reaxff`, +since with variable charges the electric potential can be known +beforehand but the energy cannot. A small additional benefit is that +the *energy* keyword requires an additional conversion to energy units +which the *potential* keyword avoids. Thus, when the *potential* +keyword is specified, the *energy* keyword must not be used. As with +*energy*, the *potential* keyword is not allowed if the added field is a +constant vector. The *potential* keyword is not supported by *fix +efield/tip4p*. + ---------- Restart, fix_modify, output, run start/stop, minimize info diff --git a/doc/src/fix_qeq_reaxff.rst b/doc/src/fix_qeq_reaxff.rst index db9015f187..00cfb6d3ce 100644 --- a/doc/src/fix_qeq_reaxff.rst +++ b/doc/src/fix_qeq_reaxff.rst @@ -128,9 +128,12 @@ periodic cell dimensions less than 10 Angstroms. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, -but there may be only one fix efield instance used, it may only use a -constant electric field, and the electric field vector may only have -components in non-periodic directions. +but there may be only one fix efield instance used and the electric field +vector may only have components in non-periodic directions. Equal-style +variables can be used for electric field vector components without any further +settings. Atom-style variables can be used for spatially-varying electric field +vector components, but the resulting electric potential must be specified +as an atom-style variable using the *potential* keyword for `fix efield`. Related commands """""""""""""""" diff --git a/src/EXTRA-FIX/fix_efield_tip4p.cpp b/src/EXTRA-FIX/fix_efield_tip4p.cpp index 3d6a2f66b7..47b1d9e27a 100644 --- a/src/EXTRA-FIX/fix_efield_tip4p.cpp +++ b/src/EXTRA-FIX/fix_efield_tip4p.cpp @@ -47,6 +47,7 @@ void FixEfieldTIP4P::init() if (atom->tag_enable == 0) error->all(FLERR, "Fix efield/tip4p requires atom IDs"); if (!atom->q_flag) error->all(FLERR, "Fix efield/tip4p requires atom attribute q"); if (!force->pair) error->all(FLERR, "A TIP4P pair style must be defined fix efield/tip4p"); + if (pstr) error->all(FLERR, "Fix efield/tip4p does not support the potential keyword"); int itmp; double *p_qdist = (double *) force->pair->extract("qdist", itmp); diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index 704b43e642..554b911151 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -404,9 +404,15 @@ void FixQEqReaxFF::init() efield->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::CONSTANT) - error->all(FLERR,"Cannot (yet) use fix {} with variable efield", 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 (((efield->xstyle != FixEfield::CONSTANT) && domain->xperiodic) || + ((efield->ystyle != FixEfield::CONSTANT) && domain->yperiodic) || + ((efield->zstyle != FixEfield::CONSTANT) && domain->zperiodic)) + error->all(FLERR,"Must not have electric field component in direction of periodic " + "boundary when using charge equilibration with ReaxFF."); if (((fabs(efield->ex) > SMALL) && domain->xperiodic) || ((fabs(efield->ey) > SMALL) && domain->yperiodic) || ((fabs(efield->ez) > SMALL) && domain->zperiodic)) @@ -1101,26 +1107,36 @@ void FixQEqReaxFF::get_chi_field() // efield energy is in real units of kcal/mol/angstrom, need to convert to eV - const double factor = -1.0/force->qe2f; + const double qe2f = force->qe2f; + const double factor = -1.0/qe2f; + + + if (efield->varflag != FixEfield::CONSTANT) + efield->update_efield_variables(); - // currently we only support constant efield // atom selection is for the group of fix efield - if (efield->varflag == FixEfield::CONSTANT) { - double unwrap[3]; - const double fx = efield->ex; - const double fy = efield->ey; - const double fz = efield->ez; - const int efgroupbit = efield->groupbit; + double unwrap[3]; + const double ex = efield->ex; + const double ey = efield->ey; + const double ez = efield->ez; + const int efgroupbit = efield->groupbit; // charge interactions // force = qE, potential energy = F dot x in unwrapped coords - + if (efield->varflag != FixEfield::ATOM) { 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); - chi_field[i] = factor*(fx*unwrap[0] + fy*unwrap[1] + fz*unwrap[2]); + chi_field[i] = factor*(ex*unwrap[0] + ey*unwrap[1] + ez*unwrap[2]); + } + } + } 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; + chi_field[i] = -efield->efield[i][3]; } } } diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 7880655973..d01a498d39 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); @@ -58,7 +58,7 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : virial_global_flag = virial_peratom_flag = 1; qe2f = force->qe2f; - xstr = ystr = zstr = nullptr; + xstyle = ystyle = zstyle = estyle = pstyle = NONE; if (utils::strmatch(arg[3], "^v_")) { xstr = utils::strdup(arg[3] + 2); @@ -100,11 +100,22 @@ 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]); } } + if (estr && pstr) + error->all(FLERR, "Must not use energy and potential keywords at the same time with fix efield"); + force_flag = 0; fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; @@ -122,6 +133,7 @@ FixEfield::~FixEfield() delete[] ystr; delete[] zstr; delete[] estr; + delete[] pstr; delete[] idregion; memory->destroy(efield); } @@ -157,43 +169,54 @@ 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); - } else - estyle = NONE; + error->all(FLERR, "Variable {} for energy in fix {} must be atom-style", estr, style); + } + + 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); + } // set index and check validity of region @@ -217,8 +240,10 @@ void FixEfield::init() if (varflag == CONSTANT && estyle != NONE) error->all(FLERR, "Cannot use variable energy with constant efield in fix {}", style); - if ((varflag == EQUAL || varflag == ATOM) && update->whichflag == 2 && estyle == NONE) - error->all(FLERR, "Must use variable energy with fix {}", style); + if (varflag == CONSTANT && pstyle != NONE) + error->all(FLERR, "Cannot use variable potential with constant efield in fix {}", style); + if ((varflag == EQUAL || varflag == ATOM) && update->whichflag == 2 && estyle == NONE && pstyle == NONE) + error->all(FLERR, "Must use variable energy or potential with fix {} during minimization", style); if (utils::strmatch(update->integrate_style, "^respa")) { ilevel_respa = (dynamic_cast(update->integrate))->nlevels - 1; @@ -346,26 +371,7 @@ void FixEfield::post_force(int vflag) } else { - modify->clearstep_compute(); - - if (xstyle == EQUAL) { - ex = qe2f * input->variable->compute_equal(xvar); - } else if (xstyle == ATOM) { - input->variable->compute_atom(xvar, igroup, &efield[0][0], 4, 0); - } - if (ystyle == EQUAL) { - ey = qe2f * input->variable->compute_equal(yvar); - } else if (ystyle == ATOM) { - input->variable->compute_atom(yvar, igroup, &efield[0][1], 4, 0); - } - if (zstyle == EQUAL) { - ez = qe2f * input->variable->compute_equal(zvar); - } 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); - - modify->addstep_compute(update->ntimestep + 1); + update_efield_variables(); // charge interactions // force = qE @@ -395,7 +401,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]; } } @@ -470,3 +477,33 @@ double FixEfield::compute_vector(int n) } return fsum_all[n + 1]; } + +/* ---------------------------------------------------------------------- + update efield variables without doing anything else + called by fix_qeq_reaxff +------------------------------------------------------------------------- */ + +void FixEfield::update_efield_variables() +{ + modify->clearstep_compute(); + + if (xstyle == EQUAL) { + ex = qe2f * input->variable->compute_equal(xvar); + } else if (xstyle == ATOM) { + input->variable->compute_atom(xvar, igroup, &efield[0][0], 4, 0); + } + if (ystyle == EQUAL) { + ey = qe2f * input->variable->compute_equal(yvar); + } else if (ystyle == ATOM) { + input->variable->compute_atom(yvar, igroup, &efield[0][1], 4, 0); + } + if (zstyle == EQUAL) { + ez = qe2f * input->variable->compute_equal(zvar); + } else if (zstyle == ATOM) { + input->variable->compute_atom(zvar, igroup, &efield[0][2], 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 52c827bb50..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; @@ -59,6 +60,7 @@ class FixEfield : public Fix { int force_flag; double fsum[4], fsum_all[4]; + void update_efield_variables(); }; } // namespace LAMMPS_NS #endif