add documentation and checks

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 (overrides *energy*)
 
 Examples
 """"""""
@@ -112,7 +114,8 @@ one or more variables, and if you are performing dynamics via the
 :doc:`run <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 <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,23 @@ 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.
 
+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 times 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<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 is ignored (the simulation
+will proceed but with a warning issued). As with *energy*, the *potential*
+keyword is not allowed if the added field is a constant vector.
+
 ----------
 
 Restart, fix_modify, output, run start/stop, minimize info

doc/src/fix_qeq_reaxff.rst

@@ -128,9 +128,13 @@ periodic cell dimensions less than 10 Angstroms.
 
 This fix may be used in combination with :doc:`fix efield <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
+but there may be only one fix efield instance used and the electric field
-constant electric field, and the electric field vector may only have
+vector may only have components in non-periodic directions. Equal-style
-components in non-periodic directions.
+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 (new) *potential* keyword for 
+`fix efield`.
 
 Related commands
 """"""""""""""""

src/REAXFF/fix_qeq_reaxff.cpp

@@ -408,6 +408,11 @@ void FixQEqReaxFF::init()
     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))

src/fix_efield.cpp

@@ -238,8 +238,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)
+  if (varflag == CONSTANT && pstyle != NONE)
-    error->all(FLERR, "Must use variable energy with fix {}", style);
+    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<Respa *>(update->integrate))->nlevels - 1;