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

Merge pull request #3793 from srtee/reaxff-varstyle-efield

Browse Source 
Enable `fix qeq/reaxff` with variable `fix efield`
This commit is contained in:
Axel Kohlmeyer
2023-06-01 10:01:56 -04:00
committed by GitHub
parent 5f8b0f56a5 90000ca166
commit bda933e0db
6 changed files with 138 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
doc/src/fix_efield.rst
View File

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

9
doc/src/fix_qeq_reaxff.rst
View File

@ -128,9 +128,12 @@ 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
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
""""""""""""""""

1
src/EXTRA-FIX/fix_efield_tip4p.cpp
View File

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

40
src/REAXFF/fix_qeq_reaxff.cpp
View File

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

107
src/fix_efield.cpp
View File

@ -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<Respa *>(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);
}

6
src/fix_efield.h
View File

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