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

Name change, add example

Browse Source
This commit is contained in:
Gabriel Alkuino
2024-11-12 11:40:07 -05:00
committed by GitHub
parent 92d4bec1a3
commit 2bcd6d1f73
11 changed files with 357 additions and 855 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_fix.rst
View File

@ -58,6 +58,7 @@ OPT.
* :doc:`dt/reset (k) <fix_dt_reset>`
* :doc:`edpd/source <fix_dpd_source>`
* :doc:`efield (k) <fix_efield>`
* :doc:`efield/lepton <fix_efield_lepton>`
* :doc:`efield/tip4p <fix_efield>`
* :doc:`ehex <fix_ehex>`
* :doc:`electrode/conp (i) <fix_electrode>`
@ -69,7 +70,6 @@ OPT.
* :doc:`eos/cv <fix_eos_cv>`
* :doc:`eos/table <fix_eos_table>`
* :doc:`eos/table/rx (k) <fix_eos_table_rx>`
* :doc:`epot/lepton <fix_epot_lepton>`
* :doc:`evaporate <fix_evaporate>`
* :doc:`external <fix_external>`
* :doc:`ffl <fix_ffl>`

2
doc/src/fix.rst
View File

@ -237,6 +237,7 @@ accelerated styles exist.
* :doc:`dt/reset <fix_dt_reset>` - reset the timestep based on velocity, forces
* :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
* :doc:`efield <fix_efield>` - impose electric field on system
* :doc:`efield/lepton <fix_efield_lepton>` - impose electric field on system using a Lepton expression for the potential
* :doc:`efield/tip4p <fix_efield>` - impose electric field on system with TIP4P molecules
* :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
* :doc:`electrode/conp <fix_electrode>` - impose electric potential
@ -248,7 +249,6 @@ accelerated styles exist.
* :doc:`eos/cv <fix_eos_cv>` - applies a mesoparticle equation of state to relate the particle internal energy to the particle internal temperature
* :doc:`eos/table <fix_eos_table>` - applies a tabulated mesoparticle equation of state to relate the particle internal energy to the particle internal temperature
* :doc:`eos/table/rx <fix_eos_table_rx>` - applies a tabulated mesoparticle equation of state to relate the concentration-dependent particle internal energy to the particle internal temperature
* :doc:`epot/lepton <fix_epot_lepton>` - apply electric potential on system using Lepton expression
* :doc:`evaporate <fix_evaporate>` - remove atoms from simulation periodically
* :doc:`external <fix_external>` - callback to an external driver program
* :doc:`ffl <fix_ffl>` - apply a Fast-Forward Langevin equation thermostat

8
doc/src/fix_efield.rst
View File

@ -45,8 +45,9 @@ Description
Add a force :math:`\vec{F} = q\vec{E}` to each charged atom in the group due to an
external electric field being applied to the system. If the system
contains point-dipoles, also add a torque on the dipoles due to the
external electric field.
contains point-dipoles, also add a torque :math:`\vec{T} = \vec{p} \times \vec{E}` on the dipoles due to the
external electric field. This fix does not compute the dipole force :math:`\vec{F} = (\vec{p} \cdot \nabla) \vec{E}`,
and the :doc:`fix efield/lepton <fix_efield_lepton>` command should be used instead.
.. versionadded:: 28Mar2023
@ -68,6 +69,7 @@ For point-dipoles, equal-style variables can be used, but atom-style
variables are not currently supported, since they imply a spatial
gradient in the electric field which means additional terms with
gradients of the field are required for the force and torque on dipoles.
The :doc:`fix efield/lepton <fix_efield_lepton>` command should be used instead.
Equal-style variables can specify formulas with various mathematical
functions, and include :doc:`thermo_style <thermo_style>` command
@ -229,7 +231,7 @@ Fix style *efield/tip4p* can only be used with tip4p pair styles.
Related commands
""""""""""""""""
:doc:`fix addforce <fix_addforce>`
:doc:`fix addforce <fix_addforce>`, :doc:`fix efield/lepton <fix_efield_lepton>`
Default
"""""""

31
doc/src/fix_epot_lepton.rst → doc/src/fix_efield_lepton.rst
View File

@ -1,17 +1,17 @@
.. index:: fix epot/lepton
.. index:: fix efield/lepton
fix epot/lepton command
=======================
fix efield/lepton command
=========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID epot/lepton V ...
fix ID group-ID efield/lepton V ...
* ID, group-ID are documented in the :doc:`fix <fix>` command
* style = *epot/lepton*
* style = *efield/lepton*
* V = electric potential (electric field * distance units)
* V must be a Lepton expression (see below)
* zero or more keyword/value pairs may be appended to args
@ -29,13 +29,13 @@ Examples
.. code-block:: LAMMPS
fix ex all epot/lepton "-E*x; E=1"
fix dexx all epot/lepton "-0.5*x^2" step 1
fix yukawa all epot/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1" step 1e-6
fix infp all epot/lepton "-abs(x)" step 1
fix ex all efield/lepton "-E*x; E=1"
fix dexx all efield/lepton "-0.5*x^2" step 1
fix yukawa all efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1" step 1e-6
fix infp all efield/lepton "-abs(x)" step 1
variable th equal 2*PI*ramp(0,1)
fix erot all epot/lepton "-(x*cos(v_th)+y*sin(v_th))"
fix erot all efield/lepton "-(x*cos(v_th)+y*sin(v_th))"
Description
"""""""""""
@ -47,7 +47,13 @@ and to dipoles a force :math:`\vec{F} = (\vec{p} \cdot \nabla) \vec{E}` and torq
where :math:`\vec{E} = - \nabla V`. The fix also evaluates the electrostatic energy (:math:`U_{q} = q V` and :math:`U_{p} = - \vec{p} \cdot \vec{E}`)
due to this potential when the :doc:`fix_modify energy yes <fix_modify>` command is specified (see below).
The `Lepton library <https://simtk.org/projects/lepton>`_, that the *epot/lepton* fix style interfaces with, evaluates
.. note::
This command should be used instead of :doc:`fix efield <fix_efield>` if you want to impose a non-uniform electric field on a system with dipoles
since the latter does not include the dipole force term. If you only have charges or if the electric field gradient is negligible,
:doc:`fix efield <fix_efield>` should be used since it is faster.
The `Lepton library <https://simtk.org/projects/lepton>`_, that the *efield/lepton* fix style interfaces with, evaluates
the expression string at run time to compute the energy, forces, and torques. It creates an analytical representation
of :math:`V` and :math:`\vec{E}`, while the gradient force is computed using a central difference scheme
@ -70,7 +76,6 @@ If the *region* keyword is used, the atom must also be in the specified
geometric :doc:`region <region>` in order to be affected by the potential.
The *step* keyword is required when :doc:`atom_style dipole <atom_style>` is used and the electric field is non-uniform.
If you only want to apply a uniform electric field, the :doc:`fix efield <fix_efield>` command will be faster.
----------
@ -123,7 +128,7 @@ the iteration count during the minimization.
Restrictions
""""""""""""
Fix style *epot/lepton* is part of the LEPTON package. It is only enabled if LAMMPS was built with that package.
Fix style *efield/lepton* is part of the LEPTON package. It is only enabled if LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` page for more info.

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -1020,7 +1020,6 @@ eos
epair
epcc
epl
epot
Epp
Epq
eps

49
examples/LEPTON/in.efield-lepton Normal file
View File

@ -0,0 +1,49 @@
# Point dipoles in a 3d box with an external potential (ignoring dipolar interactions)
units lj
atom_style hybrid sphere dipole
dimension 3
boundary s s s
region box block -2 2 -2 2 -2 2
create_box 1 box
create_atoms 1 random 100 12345 NULL
# need both mass settings due to hybrid atom style
mass 1 1.0
set group all mass 1.0
set group all diameter 0.1
set group all dipole/random 98934 0.01
pair_style none
velocity all create 0.0 87287 mom yes rot yes
fix 1 all nve/sphere update dipole
###############################################################################################################
## Yukawa potential
#fix 2 all efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A = 0.1; B = 5" step 1e-8
## Gradually increasing uniform field
#variable E equal ramp(0,1)
#fix 2 all efield/lepton "-v_E*(x+y+z)"
## Linear gradient field
fix 2 all efield/lepton "-1/6*x^3" step 1e-6
fix_modify 2 energy yes
###############################################################################################################
timestep 1e-3
compute erot all erotate/sphere
variable etotal equal "ke + c_erot + pe" # thermo etotal doesn't include erot
thermo_style custom step temp ke c_erot pe v_etotal
thermo 500
thermo_modify norm no
#dump 1 all custom 500 dump.dipole id x y z diameter mux muy muz fx fy fz tqx tqy tqz
run 10000

4
src/.gitignore vendored
View File

@ -95,8 +95,8 @@
/angle_lepton.h
/dihedral_lepton.cpp
/dihedral_lepton.h
/fix_epot_lepton.cpp
/fix_epot_lepton.h
/fix_efield_lepton.cpp
/fix_efield_lepton.h
/fix_wall_lepton.cpp
/fix_wall_lepton.h
/lepton_utils.cpp

495
src/LEPTON/fix_epot_lepton.cpp
View File

@ -1,495 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Gabriel Alkuino (Syracuse University) - gsalkuin@syr.edu
Modified from fix_efield
------------------------------------------------------------------------- */
#include "fix_epot_lepton.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "modify.h"
#include "region.h"
#include "respa.h"
#include "update.h"
#include <array>
#include "Lepton.h"
#include "lepton_utils.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define EPSILON 1.0e-10
/* ---------------------------------------------------------------------- */
FixEpotLepton::FixEpotLepton(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), idregion(nullptr), region(nullptr)
{
if (domain->xperiodic || domain->yperiodic || domain->zperiodic) {
error->warning(FLERR, "Fix {} uses unwrapped coordinates", style);
}
if (narg < 4) utils::missing_cmd_args(FLERR, std::string("fix ") + style, error);
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
energy_global_flag = 1;
virial_global_flag = virial_peratom_flag = 1;
respa_level_support = 1;
ilevel_respa = 0;
// optional args
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg], "region") == 0) {
if (iarg + 2 > narg)
utils::missing_cmd_args(FLERR, std::string("fix ") + style + " region", error);
region = domain->get_region_by_id(arg[iarg + 1]);
if (!region) error->all(FLERR, "Region {} for fix {} does not exist", arg[iarg + 1], style);
idregion = utils::strdup(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg], "step") == 0) {
if (iarg + 2 > narg)
utils::missing_cmd_args(FLERR, std::string("fix ") + style + "step", error);
h = utils::numeric(FLERR, arg[iarg+1], false, lmp);
iarg += 2;
} else {
error->all(FLERR, "Unknown keyword for fix {} command: {}", style, arg[iarg]);
}
}
// check validity of Lepton expression
// remove whitespace and quotes from expression string and then
// check if the expression can be parsed without error
expr = LeptonUtils::condense(arg[3]);
try {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
auto phi = parsed.createCompiledExpression();
} catch (std::exception &e) {
error->all(FLERR, e.what());
}
force_flag = 0;
fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
}
/* ---------------------------------------------------------------------- */
FixEpotLepton::~FixEpotLepton()
{
delete[] idregion;
}
/* ---------------------------------------------------------------------- */
int FixEpotLepton::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixEpotLepton::init()
{
if (!atom->q_flag && !atom->mu_flag)
error->all(FLERR, "Fix {} requires atom attribute q or mu", style);
if (atom->mu_flag && !atom->torque_flag)
error->all(FLERR, "Dipoles must be finite-sized to rotate", style);
// set index and check validity of region
if (idregion) {
region = domain->get_region_by_id(idregion);
if (!region) error->all(FLERR, "Region {} for fix {} does not exist", idregion, style);
}
if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
}
// unit conversion factors and restrictions (see issue #1377)
char *unit_style = update->unit_style;
qe2f = force->qe2f;
mue2e = qe2f;
if (strcmp(unit_style, "electron") == 0 || strcmp(unit_style, "micro") == 0 || strcmp(unit_style, "nano") == 0) {
error->all(FLERR, "Fix {} does not support {} units", style, unit_style);
}
}
/* ---------------------------------------------------------------------- */
void FixEpotLepton::setup(int vflag)
{
if (utils::strmatch(update->integrate_style, "^respa")) {
auto respa = dynamic_cast<Respa *>(update->integrate);
respa->copy_flevel_f(ilevel_respa);
post_force_respa(vflag, ilevel_respa, 0);
respa->copy_f_flevel(ilevel_respa);
} else {
post_force(vflag);
}
}
/* ---------------------------------------------------------------------- */
void FixEpotLepton::min_setup(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
Apply F = qE,
F = (mu . D) E,
T = mu x E
------------------------------------------------------------------------- */
void FixEpotLepton::post_force(int vflag)
{
double **f = atom->f;
double **x = atom->x;
int *mask = atom->mask;
imageint *image = atom->image;
int nlocal = atom->nlocal;
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)).optimize();
Lepton::CompiledExpression phi;
auto dphi_x = parsed.differentiate("x").createCompiledExpression();
auto dphi_y = parsed.differentiate("y").createCompiledExpression();
auto dphi_z = parsed.differentiate("z").createCompiledExpression();
std::vector<Lepton::CompiledExpression*> dphis = {&dphi_x, &dphi_y, &dphi_z};
// check if reference to x, y, z exist
const std::array<std::string, 3> variableNames = {"x", "y", "z"};
std::array<bool, 3> phi_has_ref = {true, true, true};
if (atom->q_flag){
phi = parsed.createCompiledExpression();
for (size_t i = 0; i < 3; i++) {
try {
phi.getVariableReference(variableNames[i]);
}
catch (Lepton::Exception &) {
phi_has_ref[i] = false;
}
}
}
std::vector<std::array<bool, 3>> dphis_has_ref;
bool e_uniform = true;
for (auto &dphi : dphis) {
dphis_has_ref.push_back({false, false, false});
for (size_t i = 0; i < 3; i++) {
try {
(*dphi).getVariableReference(variableNames[i]);
dphis_has_ref.back()[i] = true;
e_uniform = false;
}
catch (Lepton::Exception &) {
// do nothing
}
}
}
if (!e_uniform && atom->mu_flag && h < 0) {
error->all(FLERR, "Fix {} requires keyword `step' for dipoles in a non-uniform electric field", style);
}
// virial setup
v_init(vflag);
// update region if necessary
if (region) region->prematch();
// fsum[0] = "potential energy" for added force
// fsum[123] = extra force added to atoms
fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
force_flag = 0;
double ex, ey, ez;
double fx, fy, fz;
double v[6], unwrap[3];
double xf, yf, zf, xb, yb, zb;
double exf, eyf, ezf, exb, eyb, ezb;
double tx, ty, tz;
double mu_norm, h_mu;
if (atom->q_flag && atom->mu_flag) {
double *q = atom->q;
double **mu = atom->mu;
double **t = atom->torque;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
domain->unmap(x[i], image[i], unwrap);
// evaluate e-field, used by q and mu
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();
if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];
// charges
// force = q E
fx = qe2f * q[i] * ex;
fy = qe2f * q[i] * ey;
fz = qe2f * q[i] * ez;
// potential energy = q phi
fsum[0] += qe2f * q[i] * phi.evaluate();
// dipoles
mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
if (mu_norm > EPSILON) {
// torque = mu cross E
t[i][0] += mue2e * (ez * mu[i][1] - ey * mu[i][2]);
t[i][1] += mue2e * (ex * mu[i][2] - ez * mu[i][0]);
t[i][2] += mue2e * (ey * mu[i][0] - ex * mu[i][1]);
// potential energy = - mu dot E
fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);
// force = (mu dot D) E
// using central difference method
h_mu = h / mu_norm;
xf = unwrap[0] + h_mu * mu[i][0];
yf = unwrap[1] + h_mu * mu[i][1];
zf = unwrap[2] + h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
}
exf = -dphi_x.evaluate();
eyf = -dphi_y.evaluate();
ezf = -dphi_z.evaluate();
xb = unwrap[0] - h_mu * mu[i][0];
yb = unwrap[1] - h_mu * mu[i][1];
zb = unwrap[2] - h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
}
exb = -dphi_x.evaluate();
eyb = -dphi_y.evaluate();
ezb = -dphi_z.evaluate();
fx += qe2f * (exf - exb) / 2.0 / h_mu;
fy += qe2f * (eyf - eyb) / 2.0 / h_mu;
fz += qe2f * (ezf - ezb) / 2.0 / h_mu;
}
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
} else if (atom->q_flag && !atom->mu_flag) {
double *q = atom->q;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
domain->unmap(x[i], image[i], unwrap);
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();
if (phi_has_ref[0]) phi.getVariableReference("x") = unwrap[0];
if (phi_has_ref[1]) phi.getVariableReference("y") = unwrap[1];
if (phi_has_ref[2]) phi.getVariableReference("z") = unwrap[2];
// force = q E
fx = qe2f * q[i] * ex;
fy = qe2f * q[i] * ey;
fz = qe2f * q[i] * ez;
// potential energy = q phi
fsum[0] += qe2f * q[i] * phi.evaluate();
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
} else if (!atom->q_flag && atom->mu_flag) {
double **mu = atom->mu;
double **t = atom->torque;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue;
mu_norm = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
if (mu_norm > EPSILON) continue;
domain->unmap(x[i], image[i], unwrap);
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = unwrap[0];
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = unwrap[1];
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = unwrap[2];
}
ex = -dphi_x.evaluate();
ey = -dphi_y.evaluate();
ez = -dphi_z.evaluate();
// torque = mu cross E
t[i][0] += mue2e * (ez * mu[i][1] - ey * mu[i][2]);
t[i][1] += mue2e * (ex * mu[i][2] - ez * mu[i][0]);
t[i][2] += mue2e * (ey * mu[i][0] - ex * mu[i][1]);
// potential energy = - mu dot E
fsum[0] -= mue2e * (mu[i][0] * ex + mu[i][1] * ey + mu[i][2] * ez);
// force = (mu dot D) E
// using central difference method
h_mu = h / sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
xf = unwrap[0] + h_mu * mu[i][0];
yf = unwrap[1] + h_mu * mu[i][1];
zf = unwrap[2] + h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xf;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yf;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zf;
}
exf = -dphi_x.evaluate();
eyf = -dphi_y.evaluate();
ezf = -dphi_z.evaluate();
xb = unwrap[0] - h_mu * mu[i][0];
yb = unwrap[1] - h_mu * mu[i][1];
zb = unwrap[2] - h_mu * mu[i][2];
for (size_t j = 0; j < 3; j++) {
if (dphis_has_ref[j][0]) (*dphis[j]).getVariableReference("x") = xb;
if (dphis_has_ref[j][1]) (*dphis[j]).getVariableReference("y") = yb;
if (dphis_has_ref[j][2]) (*dphis[j]).getVariableReference("z") = zb;
}
exb = -dphi_x.evaluate();
eyb = -dphi_y.evaluate();
ezb = -dphi_z.evaluate();
fx = qe2f * (exf - exb) / 2.0 / h_mu;
fy = qe2f * (eyf - eyb) / 2.0 / h_mu;
fz = qe2f * (ezf - ezb) / 2.0 / h_mu;
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
fsum[1] += fx;
fsum[2] += fy;
fsum[3] += fz;
if (evflag) {
v[0] = fx * unwrap[0];
v[1] = fy * unwrap[1];
v[2] = fz * unwrap[2];
v[3] = fx * unwrap[1];
v[4] = fx * unwrap[2];
v[5] = fy * unwrap[2];
v_tally(i, v);
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixEpotLepton::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixEpotLepton::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
return energy added by fix
------------------------------------------------------------------------- */
double FixEpotLepton::compute_scalar()
{
if (force_flag == 0) {
MPI_Allreduce(fsum, fsum_all, 4, MPI_DOUBLE, MPI_SUM, world);
force_flag = 1;
}
return fsum_all[0];
}
/* ----------------------------------------------------------------------
return total extra force due to fix
------------------------------------------------------------------------- */
double FixEpotLepton::compute_vector(int n)
{
if (force_flag == 0) {
MPI_Allreduce(fsum, fsum_all, 4, MPI_DOUBLE, MPI_SUM, world);
force_flag = 1;
}
return fsum_all[n + 1];
}

60
src/LEPTON/fix_epot_lepton.h
View File

@ -1,60 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Gabriel Alkuino (Syracuse University) - gsalkuin@syr.edu
Modified from fix_efield
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(epot/lepton,FixEpotLepton);
// clang-format on
#else
#ifndef LMP_FIX_EPOT_LEPTON_H
#define LMP_FIX_EPOT_LEPTON_H
#include "fix.h"
namespace LAMMPS_NS {
class FixEpotLepton : public Fix {
public:
FixEpotLepton(class LAMMPS *, int, char **);
~FixEpotLepton() override;
int setmask() override;
void init() override;
void setup(int) override;
void min_setup(int) override;
void post_force(int) override;
void post_force_respa(int, int, int) override;
void min_post_force(int) override;
double compute_scalar() override;
double compute_vector(int) override;
protected:
char *idregion;
class Region *region;
int ilevel_respa;
std::string expr;
int force_flag;
double h = -1.0;
double qe2f, mue2e;
double fsum[4], fsum_all[4];
};
} // namespace LAMMPS_NS
#endif
#endif

7
unittest/force-styles/tests/fix-timestep-epot_lepton.yaml → unittest/force-styles/tests/fix-timestep-efield_lepton.yaml
View File

@ -1,6 +1,7 @@
---
lammps_version: 27 Jun 2024
date_generated: Tue Sep 17 08:38:51 2024
lammps_version: 29 Aug 2024
tags: generated
date_generated: Tue Nov 12 08:53:18 2024
epsilon: 2e-13
skip_tests:
prerequisites: ! |
@ -10,7 +11,7 @@ pre_commands: ! ""
post_commands: ! |
region half block 0 EDGE EDGE EDGE EDGE EDGE
fix move all nve
fix test solute epot/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1"
fix test solute efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1"
input_file: in.fourmol
natoms: 29
global_scalar: 1.6704454516566274

7
unittest/force-styles/tests/fix-timestep-epot_lepton_region.yaml → unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml
View File

@ -1,6 +1,7 @@
---
lammps_version: 27 Jun 2024
date_generated: Tue Sep 17 08:42:57 2024
lammps_version: 29 Aug 2024
tags: generated
date_generated: Tue Nov 12 08:53:32 2024
epsilon: 2e-13
skip_tests:
prerequisites: ! |
@ -10,7 +11,7 @@ pre_commands: ! ""
post_commands: ! |
region half block 0 EDGE EDGE EDGE EDGE EDGE
fix move all nve
fix test solute epot/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1" region half
fix test solute efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1" region half
input_file: in.fourmol
natoms: 29
global_scalar: 1.6224257709052041