lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4048 from akohlmey/lepton_updates 

Update LEPTON styles for more flexibility
This commit is contained in:
Axel Kohlmeyer 2024-01-18 16:34:38 -05:00 committed by GitHub
parent 0ad1d29c4d 425421c1ca
commit af60285a3f
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
23 changed files with 665 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
doc/src/angle_lepton.rst
View File

@ -11,7 +11,16 @@ Syntax
.. code-block:: LAMMPS
angle_style lepton
angle_style style args
* style = *lepton*
* args = optional arguments
.. parsed-literal::
args = *auto_offset* or *no_offset*
*auto_offset* = offset the potential energy so that the value at theta0 is 0.0 (default)
*no_offset* = do not offset the potential energy
Examples
""""""""
@ -19,6 +28,7 @@ Examples
.. code-block:: LAMMPS
angle_style lepton
angle_style lepton no_offset
angle_coeff 1 120.0 "k*theta^2; k=250.0"
angle_coeff 2 90.0 "k2*theta^2 + k3*theta^3 + k4*theta^4; k2=300.0; k3=-100.0; k4=50.0"
@ -41,6 +51,13 @@ angle coefficient. For example `"200.0*theta^2"` represents a
U_{angle,i} = K (\theta_i - \theta_0)^2 = K \theta^2 \qquad \theta = \theta_i - \theta_0
.. versionchanged:: TBD
By default the potential energy U is shifted so that the value U is 0.0
for $theta = theta_0$. This is equivalent to using the optional keyword
*auto_offset*. When using the keyword *no_offset* instead, the
potential energy is not shifted.
The `Lepton library <https://simtk.org/projects/lepton>`_, that the
*lepton* angle style interfaces with, evaluates this expression string
at run time to compute the pairwise energy. It also creates an

19
doc/src/bond_lepton.rst
View File

@ -11,7 +11,16 @@ Syntax
.. code-block:: LAMMPS
bond_style lepton
bond_style style args
* style = *lepton*
* args = optional arguments
.. parsed-literal::
args = *auto_offset* or *no_offset*
*auto_offset* = offset the potential energy so that the value at r0 is 0.0 (default)
*no_offset* = do not offset the potential energy
Examples
""""""""
@ -19,6 +28,7 @@ Examples
.. code-block:: LAMMPS
bond_style lepton
bond_style lepton no_offset
bond_coeff 1 1.5 "k*r^2; k=250.0"
bond_coeff 2 1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
@ -40,6 +50,13 @@ constant *K* of 200.0 energy units:
U_{bond,i} = K (r_i - r_0)^2 = K r^2 \qquad r = r_i - r_0
.. versionchanged:: TBD
By default the potential energy U is shifted so that he value U is 0.0
for $r = r_0$. This is equivalent to using the optional keyword
*auto_offset*. When using the keyword *no_offset* instead, the
potential energy is not shifted.
The `Lepton library <https://simtk.org/projects/lepton>`_, that the
*lepton* bond style interfaces with, evaluates this expression string at
run time to compute the pairwise energy. It also creates an analytical

6
doc/src/pair_lepton.rst
View File

@ -72,7 +72,7 @@ interactions between particles which depend on the distance and have a
cutoff. The potential function must be provided as an expression string
using "r" as the distance variable. With pair style *lepton/coul* one
may additionally reference the charges of the two atoms of the pair with
"qi" and "qj", respectively. With pair style *lepton/coul* one may
"qi" and "qj", respectively. With pair style *lepton/sphere* one may
instead reference the radii of the two atoms of the pair with "radi" and
"radj", respectively; this is half of the diameter that can be set in
:doc:`data files <read_data>` or the :doc:`set command <set>`.
@ -166,8 +166,8 @@ mixing. Thus, expressions for *all* I,J pairs must be specified
explicitly.
Only pair style *lepton* supports the :doc:`pair_modify shift <pair_modify>`
option for shifting the energy of the pair interaction so that it is
0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*.
option for shifting the potential energy of the pair interaction so that
it is 0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*.
The :doc:`pair_modify table <pair_modify>` options are not relevant for
the these pair styles.

74
src/LEPTON/angle_lepton.cpp
View File

@ -44,6 +44,7 @@ AngleLepton::AngleLepton(LAMMPS *_lmp) :
{
writedata = 1;
reinitflag = 0;
auto_offset = 1;
}
/* ---------------------------------------------------------------------- */
@ -90,10 +91,21 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
{
std::vector<Lepton::CompiledExpression> angleforce;
std::vector<Lepton::CompiledExpression> anglepot;
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
has_ref.push_back(true);
try {
angleforce.back().getVariableReference("theta");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
error->all(FLERR, e.what());
}
const double *const *const x = atom->x;
@ -142,8 +154,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
const double dtheta = acos(c) - theta0[type];
const int idx = type2expression[type];
angleforce[idx].getVariableReference("theta") = dtheta;
if (has_ref[idx]) angleforce[idx].getVariableReference("theta") = dtheta;
const double a = -angleforce[idx].evaluate() * s;
const double a11 = a * c / rsq1;
const double a12 = -a / (r1 * r2);
@ -179,7 +190,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
double eangle = 0.0;
if (EFLAG) {
anglepot[idx].getVariableReference("theta") = dtheta;
try {
anglepot[idx].getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant force
}
eangle = anglepot[idx].evaluate() - offset[type];
}
if (EVFLAG)
@ -202,6 +217,24 @@ void AngleLepton::allocate()
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void AngleLepton::settings(int narg, char **arg)
{
auto_offset = 1;
if (narg > 0) {
if (strcmp(arg[0],"auto_offset") == 0) {
auto_offset = 1;
} else if (strcmp(arg[0],"no_offset") == 0) {
auto_offset = 0;
} else {
error->all(FLERR, "Unknown angle style lepton setting {}", arg[0]);
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
@ -224,9 +257,20 @@ void AngleLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto anglepot = parsed.createCompiledExpression();
auto angleforce = parsed.differentiate("theta").createCompiledExpression();
anglepot.getVariableReference("theta") = 0.0;
angleforce.getVariableReference("theta") = 0.0;
offset_one = anglepot.evaluate();
try {
anglepot.getVariableReference("theta") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'theta'", exp_one);
}
try {
angleforce.getVariableReference("theta") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'theta'",
exp_one);
}
if (auto_offset) offset_one = anglepot.evaluate();
angleforce.evaluate();
} catch (std::exception &e) {
error->all(FLERR, e.what());
@ -284,6 +328,7 @@ void AngleLepton::write_restart(FILE *fp)
fwrite(&n, sizeof(int), 1, fp);
fwrite(exp.c_str(), sizeof(char), n, fp);
}
fwrite(&auto_offset, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -323,6 +368,9 @@ void AngleLepton::read_restart(FILE *fp)
expressions.emplace_back(buf);
}
if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
delete[] buf;
}
@ -363,7 +411,11 @@ double AngleLepton::single(int type, int i1, int i2, int i3)
const auto &expr = expressions[type2expression[type]];
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
auto anglepot = parsed.createCompiledExpression();
anglepot.getVariableReference("theta") = dtheta;
try {
anglepot.getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
return anglepot.evaluate() - offset[type];
}

2
src/LEPTON/angle_lepton.h
View File

@ -29,6 +29,7 @@ class AngleLepton : public Angle {
AngleLepton(class LAMMPS *);
~AngleLepton() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
double equilibrium_angle(int) override;
void write_restart(FILE *) override;
@ -42,6 +43,7 @@ class AngleLepton : public Angle {
double *theta0;
int *type2expression;
double *offset;
int auto_offset;
virtual void allocate();

62
src/LEPTON/bond_lepton.cpp
View File

@ -37,6 +37,7 @@ BondLepton::BondLepton(LAMMPS *_lmp) :
{
writedata = 1;
reinitflag = 0;
auto_offset = 1;
}
/* ---------------------------------------------------------------------- */
@ -82,10 +83,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
{
std::vector<Lepton::CompiledExpression> bondforce;
std::vector<Lepton::CompiledExpression> bondpot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back(true);
try {
bondforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -116,7 +124,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
double fbond = 0.0;
if (r > 0.0) {
bondforce[idx].getVariableReference("r") = dr;
if (has_ref[idx]) bondforce[idx].getVariableReference("r") = dr;
fbond = -bondforce[idx].evaluate() / r;
}
@ -136,7 +144,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
double ebond = 0.0;
if (EFLAG) {
bondpot[idx].getVariableReference("r") = dr;
try {
bondpot[idx].getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
ebond = bondpot[idx].evaluate() - offset[type];
}
if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz);
@ -157,6 +169,24 @@ void BondLepton::allocate()
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void BondLepton::settings(int narg, char **arg)
{
auto_offset = 1;
if (narg > 0) {
if (strcmp(arg[0],"auto_offset") == 0) {
auto_offset = 1;
} else if (strcmp(arg[0],"no_offset") == 0) {
auto_offset = 0;
} else {
error->all(FLERR, "Unknown bond style lepton setting {}", arg[0]);
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
@ -179,9 +209,19 @@ void BondLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto bondpot = parsed.createCompiledExpression();
auto bondforce = parsed.differentiate("r").createCompiledExpression();
bondpot.getVariableReference("r") = 0.0;
bondforce.getVariableReference("r") = 0.0;
offset_one = bondpot.evaluate();
try {
bondpot.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &e) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'r'", exp_one);
}
try {
bondforce.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &e) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one);
}
if (auto_offset) offset_one = bondpot.evaluate();
bondforce.evaluate();
} catch (std::exception &e) {
error->all(FLERR, e.what());
@ -239,6 +279,7 @@ void BondLepton::write_restart(FILE *fp)
fwrite(&n, sizeof(int), 1, fp);
fwrite(exp.c_str(), sizeof(char), n, fp);
}
fwrite(&auto_offset, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -278,6 +319,9 @@ void BondLepton::read_restart(FILE *fp)
expressions.emplace_back(buf);
}
if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
delete[] buf;
}
@ -302,8 +346,12 @@ double BondLepton::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
auto bondpot = parsed.createCompiledExpression();
auto bondforce = parsed.differentiate("r").createCompiledExpression();
bondforce.getVariableReference("r") = dr;
bondpot.getVariableReference("r") = dr;
try {
bondpot.getVariableReference("r") = dr;
bondforce.getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
// force and energy

2
src/LEPTON/bond_lepton.h
View File

@ -29,6 +29,7 @@ class BondLepton : public Bond {
BondLepton(class LAMMPS *);
~BondLepton() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
double equilibrium_distance(int) override;
void write_restart(FILE *) override;
@ -42,6 +43,7 @@ class BondLepton : public Bond {
double *r0;
int *type2expression;
double *offset;
int auto_offset;
virtual void allocate();

29
src/LEPTON/dihedral_lepton.cpp
View File

@ -92,10 +92,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
{
std::vector<Lepton::CompiledExpression> dihedralforce;
std::vector<Lepton::CompiledExpression> dihedralpot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression());
has_ref.push_back(true);
try {
dihedralforce.back().getVariableReference("phi");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -278,7 +285,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
}
const int idx = type2expression[type];
dihedralforce[idx].getVariableReference("phi") = phi;
if (has_ref[idx]) dihedralforce[idx].getVariableReference("phi") = phi;
double m_du_dphi = -dihedralforce[idx].evaluate();
// ----- Step 4: Calculate the force direction in real space -----
@ -322,7 +329,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void DihedralLepton::eval()
double edihedral = 0.0;
if (EFLAG) {
dihedralpot[idx].getVariableReference("phi") = phi;
try {
dihedralpot[idx].getVariableReference("phi") = phi;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
edihedral = dihedralpot[idx].evaluate();
}
if (EVFLAG)
@ -362,8 +373,18 @@ void DihedralLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto dihedralpot = parsed.createCompiledExpression();
auto dihedralforce = parsed.differentiate("phi").createCompiledExpression();
dihedralpot.getVariableReference("phi") = 0.0;
dihedralforce.getVariableReference("phi") = 0.0;
try {
dihedralpot.getVariableReference("phi") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'phi'", exp_one);
}
try {
dihedralforce.getVariableReference("phi") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'phi'", exp_one);
}
dihedralforce.evaluate();
} catch (std::exception &e) {
error->all(FLERR, e.what());

15
src/LEPTON/fix_wall_lepton.cpp
View File

@ -13,6 +13,7 @@
#include "fix_wall_lepton.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "Lepton.h"
@ -41,8 +42,18 @@ void FixWallLepton::post_constructor()
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
auto wallpot = parsed.createCompiledExpression();
auto wallforce = parsed.differentiate("r").createCompiledExpression();
wallpot.getVariableReference("r") = 0.0;
wallforce.getVariableReference("r") = 0.0;
try {
wallpot.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Lepton potential expression {} does not depend on 'r'", exp_one);
}
try {
wallforce.getVariableReference("r") = 0.0;
} catch (Lepton::Exception &) {
if (comm->me == 0)
error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one);
}
wallpot.evaluate();
wallforce.evaluate();
} catch (std::exception &e) {

41
src/LEPTON/pair_lepton.cpp
View File

@ -27,6 +27,7 @@
#include "Lepton.h"
#include "lepton_utils.h"
#include <array>
#include <cmath>
#include <map>
@ -105,11 +106,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
pairforce.back().getVariableReference("r");
has_ref.push_back(true);
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -142,8 +149,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
double &r_for = pairforce[idx].getVariableReference("r");
r_for = r;
if (has_ref[idx]) pairforce[idx].getVariableReference("r") = r;
const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair;
@ -157,7 +163,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
double evdwl = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
evdwl *= factor_lj;
}
@ -229,8 +239,12 @@ void PairLepton::coeff(int narg, char **arg)
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp), functions);
auto pairforce = parsed.differentiate("r").createCompiledExpression();
auto pairpot = parsed.createCompiledExpression();
pairpot.getVariableReference("r") = 1.0;
pairforce.getVariableReference("r") = 1.0;
try {
pairpot.getVariableReference("r") = 1.0;
pairforce.getVariableReference("r") = 1.0;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
pairpot.evaluate();
pairforce.evaluate();
} catch (std::exception &e) {
@ -270,7 +284,11 @@ double PairLepton::init_one(int i, int j)
try {
auto expr = LeptonUtils::substitute(expressions[type2expression[i][j]], lmp);
auto pairpot = Lepton::Parser::parse(expr, functions).createCompiledExpression();
pairpot.getVariableReference("r") = cut[i][j];
try {
pairpot.getVariableReference("r") = cut[i][j];
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
offset[i][j] = pairpot.evaluate();
} catch (std::exception &) {
}
@ -429,9 +447,12 @@ double PairLepton::single(int /* i */, int /* j */, int itype, int jtype, double
auto pairforce = parsed.differentiate("r").createCompiledExpression();
const double r = sqrt(rsq);
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
try {
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
fforce = -pairforce.evaluate() / r * factor_lj;
return (pairpot.evaluate() - offset[itype][jtype]) * factor_lj;
}

50
src/LEPTON/pair_lepton_coul.cpp
View File

@ -28,6 +28,8 @@
#include "Lepton.h"
#include "lepton_utils.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS;
@ -79,25 +81,30 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_q;
std::vector<std::array<bool, 3>> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_q.emplace_back(true, true);
// check if there are references to charges
try {
pairforce.back().getVariableReference("qi");
} catch (std::exception &) {
have_q.back().first = false;
} catch (Lepton::Exception &) {
has_ref.back()[1] = false;
}
try {
pairforce.back().getVariableReference("qj");
} catch (std::exception &) {
have_q.back().second = false;
} catch (Lepton::Exception &) {
has_ref.back()[2] = false;
}
}
} catch (std::exception &e) {
@ -130,9 +137,9 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairforce[idx].getVariableReference("qj") = q2e * q[j];
if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (has_ref[idx][1]) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairforce[idx].getVariableReference("qj") = q2e * q[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_coul;
fxtmp += delx * fpair;
@ -146,9 +153,14 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
double ecoul = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairpot[idx].getVariableReference("qj") = q2e * q[j];
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("qj") = q2e * q[j];
ecoul = pairpot[idx].evaluate();
ecoul *= factor_coul;
}
@ -249,18 +261,22 @@ double PairLeptonCoul::single(int i, int j, int itype, int jtype, double rsq, do
const double r = sqrt(rsq);
const double q2e = sqrt(force->qqrd2e);
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
try {
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
try {
pairpot.getVariableReference("qi") = q2e * atom->q[i];
pairforce.getVariableReference("qi") = q2e * atom->q[i];
} catch (std::exception &) {
} catch (Lepton::Exception &) {
/* ignore */
}
try {
pairpot.getVariableReference("qj") = q2e * atom->q[j];
pairforce.getVariableReference("qj") = q2e * atom->q[j];
} catch (std::exception &) {
} catch (Lepton::Exception &) {
/* ignore */
}

55
src/LEPTON/pair_lepton_sphere.cpp
View File

@ -28,6 +28,7 @@
#include "Lepton.h"
#include "lepton_utils.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS;
@ -77,25 +78,30 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_rad;
std::vector<std::array<bool, 3>> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_rad.emplace_back(true, true);
// check if there are references to charges
// check if there are references to radii
try {
pairforce.back().getVariableReference("radi");
} catch (std::exception &) {
have_rad.back().first = false;
} catch (Lepton::Exception &) {
has_ref.back()[1] = false;
}
try {
pairforce.back().getVariableReference("radj");
} catch (std::exception &) {
have_rad.back().second = false;
} catch (Lepton::Exception &) {
has_ref.back()[2] = false;
}
}
} catch (std::exception &e) {
@ -128,9 +134,9 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairforce[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairforce[idx].getVariableReference("radj") = radius[j];
if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (has_ref[idx][1]) pairforce[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairforce[idx].getVariableReference("radj") = radius[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair;
@ -144,9 +150,14 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonSphere::eval()
double evdwl = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairpot[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairpot[idx].getVariableReference("radj") = radius[j];
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("radj") = radius[j];
evdwl = pairpot[idx].evaluate();
evdwl *= factor_lj;
}
@ -211,19 +222,23 @@ double PairLeptonSphere::single(int i, int j, int itype, int jtype, double rsq,
auto pairforce = parsed.differentiate("r").createCompiledExpression();
const double r = sqrt(rsq);
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
try {
pairpot.getVariableReference("r") = r;
pairforce.getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential or force
}
try {
pairpot.getVariableReference("radi") = atom->radius[i];
pairforce.getVariableReference("radi") = atom->radius[i];
} catch (std::exception &) {
/* ignore */
} catch (Lepton::Exception &) {
; // ignore
}
try {
pairpot.getVariableReference("radj") = atom->radius[j];
pairforce.getVariableReference("radj") = atom->radius[j];
} catch (std::exception &) {
/* ignore */
} catch (Lepton::Exception &) {
; // ignore
}
fforce = -pairforce.evaluate() / r * factor_lj;

16
src/OPENMP/angle_lepton_omp.cpp
View File

@ -91,10 +91,17 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{
std::vector<Lepton::CompiledExpression> angleforce;
std::vector<Lepton::CompiledExpression> anglepot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
has_ref.push_back(true);
try {
angleforce.back().getVariableReference("theta");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -146,8 +153,7 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
const double dtheta = acos(c) - theta0[type];
const int idx = type2expression[type];
angleforce[idx].getVariableReference("theta") = dtheta;
if (has_ref[idx]) angleforce[idx].getVariableReference("theta") = dtheta;
const double a = -angleforce[idx].evaluate() * s;
const double a11 = a * c / rsq1;
const double a12 = -a / (r1 * r2);
@ -183,7 +189,11 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double eangle = 0.0;
if (EFLAG) {
anglepot[idx].getVariableReference("theta") = dtheta;
try {
anglepot[idx].getVariableReference("theta") = dtheta;
} catch (Lepton::Exception &) {
; // ignore -> constant force
}
eangle = anglepot[idx].evaluate() - offset[type];
}
if (EVFLAG)

15
src/OPENMP/bond_lepton_omp.cpp
View File

@ -89,10 +89,17 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{
std::vector<Lepton::CompiledExpression> bondforce;
std::vector<Lepton::CompiledExpression> bondpot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back(true);
try {
bondforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -122,7 +129,7 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double fbond = 0.0;
if (r > 0.0) {
bondforce[idx].getVariableReference("r") = dr;
if (has_ref[idx]) bondforce[idx].getVariableReference("r") = dr;
fbond = -bondforce[idx].evaluate() / r;
}
@ -142,7 +149,11 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double ebond = 0.0;
if (EFLAG) {
bondpot[idx].getVariableReference("r") = dr;
try {
bondpot[idx].getVariableReference("r") = dr;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
ebond = bondpot[idx].evaluate() - offset[type];
}
if (EVFLAG)

19
src/OPENMP/dihedral_lepton_omp.cpp
View File

@ -19,9 +19,9 @@
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "math_extra.h"
#include "neighbor.h"
#include "suffix.h"
#include "math_extra.h"
#include <cmath>
@ -94,10 +94,17 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
{
std::vector<Lepton::CompiledExpression> dihedralforce;
std::vector<Lepton::CompiledExpression> dihedralpot;
std::vector<bool> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
dihedralforce.emplace_back(parsed.differentiate("phi").createCompiledExpression());
has_ref.push_back(true);
try {
dihedralforce.back().getVariableReference("phi");
} catch (Lepton::Exception &) {
has_ref.back() = false;
}
if (EFLAG) dihedralpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -106,7 +113,7 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
const double *const *const x = atom->x;
auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
const int * const * const dihedrallist = neighbor->dihedrallist;
const int *const *const dihedrallist = neighbor->dihedrallist;
const int nlocal = atom->nlocal;
// The dihedral angle "phi" is the angle between n123 and n234
@ -279,7 +286,7 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
}
const int idx = type2expression[type];
dihedralforce[idx].getVariableReference("phi") = phi;
if (has_ref[idx]) dihedralforce[idx].getVariableReference("phi") = phi;
double m_du_dphi = -dihedralforce[idx].evaluate();
// ----- Step 4: Calculate the force direction in real space -----
@ -323,7 +330,11 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
double edihedral = 0.0;
if (EFLAG) {
dihedralpot[idx].getVariableReference("phi") = phi;
try {
dihedralpot[idx].getVariableReference("phi") = phi;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
edihedral = dihedralpot[idx].evaluate();
}
if (EVFLAG)

42
src/OPENMP/pair_lepton_coul_omp.cpp
View File

@ -20,11 +20,13 @@
#include "neigh_list.h"
#include "suffix.h"
#include <cmath>
#include "Lepton.h"
#include "lepton_utils.h"
#include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -101,25 +103,30 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_q;
std::vector<std::array<bool, 3>> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_q.emplace_back(true, true);
// check if there are references to charges
try {
pairforce.back().getVariableReference("qi");
} catch (std::exception &) {
have_q.back().first = false;
} catch (Lepton::Exception &) {
has_ref.back()[1] = false;
}
try {
pairforce.back().getVariableReference("qj");
} catch (std::exception &) {
have_q.back().second = false;
} catch (Lepton::Exception &) {
has_ref.back()[2] = false;
}
}
} catch (std::exception &e) {
@ -152,9 +159,9 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairforce[idx].getVariableReference("qj") = q2e * q[j];
if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (has_ref[idx][1]) pairforce[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairforce[idx].getVariableReference("qj") = q2e * q[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_coul;
fxtmp += delx * fpair;
@ -168,9 +175,14 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
double ecoul = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
if (have_q[idx].first) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (have_q[idx].second) pairpot[idx].getVariableReference("qj") = q2e * q[j];
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("qi") = q2e * q[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("qj") = q2e * q[j];
ecoul = pairpot[idx].evaluate();
ecoul *= factor_coul;
}

20
src/OPENMP/pair_lepton_omp.cpp
View File

@ -20,11 +20,12 @@
#include "neigh_list.h"
#include "suffix.h"
#include <cmath>
#include "Lepton.h"
#include "lepton_utils.h"
#include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -96,10 +97,17 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<bool> have_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
have_ref.push_back(true);
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
have_ref.back() = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
}
} catch (std::exception &e) {
@ -132,7 +140,7 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r;
if (have_ref[idx]) pairforce[idx].getVariableReference("r") = r;
const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair;
@ -146,7 +154,11 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
double evdwl = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
evdwl *= factor_lj;
}

44
src/OPENMP/pair_lepton_sphere_omp.cpp
View File

@ -20,11 +20,13 @@
#include "neigh_list.h"
#include "suffix.h"
#include <cmath>
#include "Lepton.h"
#include "lepton_utils.h"
#include "omp_compat.h"
#include <array>
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -99,25 +101,30 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
std::vector<Lepton::CompiledExpression> pairforce;
std::vector<Lepton::CompiledExpression> pairpot;
std::vector<std::pair<bool, bool>> have_rad;
std::vector<std::array<bool, 3>> has_ref;
try {
for (const auto &expr : expressions) {
auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
has_ref.push_back({true, true, true});
try {
pairforce.back().getVariableReference("r");
} catch (Lepton::Exception &) {
has_ref.back()[0] = false;
}
if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
pairforce.back().getVariableReference("r");
have_rad.emplace_back(true, true);
// check if there are references to charges
// check if there are references to radii
try {
pairforce.back().getVariableReference("radi");
} catch (std::exception &) {
have_rad.back().first = false;
} catch (Lepton::Exception &) {
has_ref.back()[1] = false;
}
try {
pairforce.back().getVariableReference("radj");
} catch (std::exception &) {
have_rad.back().second = false;
} catch (Lepton::Exception &) {
has_ref.back()[2] = false;
}
}
} catch (std::exception &e) {
@ -150,9 +157,9 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
if (rsq < cutsq[itype][jtype]) {
const double r = sqrt(rsq);
const int idx = type2expression[itype][jtype];
pairforce[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairforce[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairforce[idx].getVariableReference("radj") = radius[j];
if (has_ref[idx][0]) pairforce[idx].getVariableReference("r") = r;
if (has_ref[idx][1]) pairforce[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairforce[idx].getVariableReference("radj") = radius[j];
const double fpair = -pairforce[idx].evaluate() / r * factor_lj;
fxtmp += delx * fpair;
@ -166,9 +173,14 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr)
double evdwl = 0.0;
if (EFLAG) {
pairpot[idx].getVariableReference("r") = r;
if (have_rad[idx].first) pairpot[idx].getVariableReference("radi") = radius[i];
if (have_rad[idx].second) pairpot[idx].getVariableReference("radj") = radius[j];
try {
pairpot[idx].getVariableReference("r") = r;
} catch (Lepton::Exception &) {
; // ignore -> constant potential
}
if (has_ref[idx][1]) pairpot[idx].getVariableReference("radi") = radius[i];
if (has_ref[idx][2]) pairpot[idx].getVariableReference("radj") = radius[j];
evdwl = pairpot[idx].evaluate();
evdwl *= factor_lj;
}

2
unittest/force-styles/tests/angle-lepton.yaml
View File

@ -9,7 +9,7 @@ prerequisites: ! |
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
angle_style: lepton
angle_style: lepton auto_offset
angle_coeff: ! |
1 110.1 "k*theta^2; k=75.0"
2 111.0 "k*theta^2; k=45.0"

88
unittest/force-styles/tests/angle-lepton_nooffset.yaml Normal file
View File

@ -0,0 +1,88 @@
---
lammps_version: 22 Dec 2022
date_generated: Fri Dec 23 15:10:29 2022
epsilon: 7.5e-13
skip_tests:
prerequisites: ! |
atom full
angle lepton
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
angle_style: lepton no_offset
angle_coeff: ! |
1 110.1 "k*theta^2; k=75.0"
2 111.0 "k*theta^2; k=45.0"
3 120.0 "k*theta^2; k=50.0"
4 108.5 "k*theta^2; k=100.0"
equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476
extract: ! |
theta0 1
natoms: 29
init_energy: 41.53081789649104
init_stress: ! |2-
8.9723357320869297e+01 -8.7188643750026529e+01 -2.5347135708427655e+00 9.2043419883119782e+01 -2.8187238090404904e+01 -1.5291148024926793e+00
init_forces: ! |2
1 4.7865489310693540e+01 7.8760925902181516e+00 -3.2694525514709866e+01
2 -1.1124882516177341e+00 -9.0075464203887403e+00 -7.2431691227364459e+00
3 -5.9057050592859328e+00 5.3263619873546261e+01 5.2353380124691469e+01
4 -1.6032230038990633e+01 -2.4560529343731403e+01 1.2891625920422307e+01
5 -4.4802331573497639e+01 -4.8300919461089379e+01 -2.3310767889219324e+01
6 4.7083124388174824e+01 -9.5212933434476312e+00 -3.2526392870546800e+01
7 -1.6208182775476303e+01 1.4458587960739102e+01 -3.5314745459502710e+00
8 -6.5664612141881040e+00 -2.5126850154274202e+01 8.2187944731423329e+01
9 -1.5504395262358301e+01 1.6121044185227817e+01 -4.2007069622477866e-01
10 9.9863759179365275e+00 4.1873540105704549e+01 -6.6085640966037403e+01
11 -2.0441876158908627e+01 -6.5186824168985984e+00 9.0023620309811072e+00
12 -1.0772126658369565e+01 -1.0807367300158219e+01 -9.6049647456797871e+00
13 2.8847886813946291e+00 7.2973241014859198e+00 -1.0414233993842981e-01
14 1.5267407478336393e+01 -9.4754911480231776e+00 -6.6307012925544200e+00
15 1.2402914209534773e+01 -6.2644630791613967e+00 1.8484576795819933e+01
16 3.8927757686508357e-01 1.0690061587911176e+01 6.1542759189377696e+00
17 1.4664194297570785e+00 -1.9971277376602425e+00 1.0776844613215999e+00
18 1.5785371874873322e-01 1.6495665212200166e+00 -6.6944747776990434e+00
19 -1.9328033033421670e+00 -2.4078805870919706e+00 2.8669575541313534e+00
20 1.7749495845934338e+00 7.5831406587195394e-01 3.8275172235676900e+00
21 3.4186149299343742e+00 4.2795410364249484e+00 -1.2789555411020650e+01
22 -6.0875600315279677e+00 -4.1504951869796605e+00 4.5212856070195766e+00
23 2.6689451015935934e+00 -1.2904584944528752e-01 8.2682698040010738e+00
24 -1.3053945393770587e+00 5.0741459325183271e+00 -3.0209518576073018e+00
25 -1.0471133765834284e+00 -3.5082261409793856e+00 5.7374874908501228e-01
26 2.3525079159604871e+00 -1.5659197915389413e+00 2.4472031085222894e+00
27 -2.8720725187343754e-01 2.3577465459557132e+00 -8.0312673032168869e-01
28 -6.2799575211500369e-01 -1.4097313073755862e+00 3.2747938980616453e-02
29 9.1520300398844123e-01 -9.4801523858012704e-01 7.7037879134107223e-01
run_energy: 41.28323739029462
run_stress: ! |2-
8.8236221596506681e+01 -8.6492260623309562e+01 -1.7439609731970940e+00 9.0601855980531312e+01 -2.8735005690484968e+01 -2.6097632235197477e+00
run_forces: ! |2
1 4.7316793853445830e+01 8.2815577813110188e+00 -3.2021703111755464e+01
2 -1.1508196824491330e+00 -9.3814982172707460e+00 -7.5761211707510139e+00
3 -5.1083163691832576e+00 5.2667553294971619e+01 5.1784852458007592e+01
4 -1.6078177452605999e+01 -2.4156048365236213e+01 1.3140924677013103e+01
5 -4.4915734474022280e+01 -4.8095168640411821e+01 -2.3331149037574161e+01
6 4.7077916942842350e+01 -9.5906213020090156e+00 -3.2570331503075487e+01
7 -1.6228599672412471e+01 1.4485102617342370e+01 -3.5441153194985300e+00
8 -6.5097893981550730e+00 -2.5117582302614530e+01 8.2131369512416001e+01
9 -1.5527440970965937e+01 1.6147270375910470e+01 -4.0812004993325646e-01
10 1.0070812216240984e+01 4.1571532807578805e+01 -6.5968810328796337e+01
11 -2.0431584971707451e+01 -6.4817395192247664e+00 8.9879981618991636e+00
12 -1.0884695976714678e+01 -1.1067390190389006e+01 -9.1551242768940568e+00
13 2.8052913970098801e+00 7.1296301666594912e+00 1.3173039168682621e-02
14 1.5254877537873529e+01 -8.9700095533297350e+00 -6.5719846903613162e+00
15 1.2392009100170984e+01 -6.0827695435257292e+00 1.7929674392339596e+01
16 4.7158712437377481e-01 1.0631038523396533e+01 6.0960085687560355e+00
17 1.4458707962589659e+00 -1.9708579331587350e+00 1.0634586790394520e+00
18 1.4201882413835909e-01 1.4265339757773337e+00 -5.7663956896747992e+00
19 -1.6609130686729365e+00 -2.0735307593211125e+00 2.4755525101127143e+00
20 1.5188942445345774e+00 6.4699678354377899e-01 3.2908431795620849e+00
21 3.2242729509516406e+00 4.0079233768386153e+00 -1.2047892238650988e+01
22 -5.7215184687399772e+00 -3.8871624402883409e+00 4.2679223469272234e+00
23 2.4972455177883366e+00 -1.2076093655027398e-01 7.7799698917237645e+00
24 -1.1661978296905471e+00 4.5271404898674854e+00 -2.6925565853370195e+00
25 -9.2712094527152167e-01 -3.1291890525017125e+00 5.1208215565053827e-01
26 2.0933187749620688e+00 -1.3979514373657731e+00 2.1804744296864813e+00
27 -2.6804542538020537e-01 2.1830651328698103e+00 -7.3931790038945400e-01
28 -5.7927072943128310e-01 -1.3052929090347909e+00 2.8365455885795865e-02
29 8.4731615481148848e-01 -8.7777222383501941e-01 7.1095244450365813e-01
...

89
unittest/force-styles/tests/bond-lepton_nooffset.yaml Normal file
View File

@ -0,0 +1,89 @@
---
lammps_version: 21 Nov 2023
date_generated: Thu Jan 18 10:15:41 2024
epsilon: 2.5e-13
skip_tests:
prerequisites: ! |
atom full
bond lepton
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
bond_style: lepton no_offset
bond_coeff: ! |
1 1.5 "k*r^2; k=250.0"
2 1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
3 1.3 "k*r^2; k=350.0"
4 1.2 "k*(r-0.2)^2; k=500.0"
5 1.0 "k*r^2; k=450.0"
equilibrium: 5 1.5 1.1 1.3 1.2 1
extract: ! |
r0 1
natoms: 29
init_energy: 38.295825321689215
init_stress: ! |-
-4.7778964706834920e+01 -9.3066674567350432e+01 3.4789470658440035e+02 -3.0023920169312170e+01 -8.0421418879842847e+01 5.8592449335969732e+01
init_forces: ! |2
1 -5.9149914305071416e+00 -3.7728809612345245e+01 -2.7769433362963369e+01
2 -9.4281609567839944e+00 -7.7586487054273015e+00 1.1096676787527940e+01
3 3.2211742366572125e+01 2.7682361264425523e+01 -7.0109911672970497e+00
4 4.9260777576375503e+00 -1.3809750102765932e+00 3.4951785613141868e+00
5 -1.2606902198593501e+00 -1.9373397933007170e+00 6.4372463095041841e+00
6 -3.8858476307965482e+01 6.8567296300319640e+01 1.9889888806614337e+02
7 7.5297927100028144e+00 -3.8622600737556944e+01 -1.9268793182212875e+02
8 1.3018665172824681e+01 -1.2902789438539877e+01 3.2406676637830003e+00
9 7.4343536239661590e-01 8.0072549738604493e-01 3.2899591078538779e+00
10 6.1558871886113291e+00 -2.2419470219698296e+00 1.0080175092279852e+01
11 -3.7020922615305768e-01 -9.1704102274126453e-01 -1.5046795827370363e+00
12 5.2437190958790678e+00 3.4225915524442998e+00 -2.5523597276998897e+00
13 -1.1277007635800260e+01 4.4610677459696646e+00 2.1195215396108269e-01
14 2.9813926585641828e+00 -6.0667387499775116e-01 7.7317115100728788e+00
15 2.5872825164662799e-01 -9.9415365173790704e+00 -3.5428115826174169e+00
16 5.2775953236493464e+01 -3.1855535724919463e+01 -1.6524229620195118e+02
17 -5.8735858023559175e+01 4.0959855098908882e+01 1.5582804819495431e+02
18 -9.0963607969319646e+00 -4.3343406270234155e+00 -1.7623055551859267e+01
19 1.2597490501067170e+01 8.0591915019111742e+00 1.5261489294231819e+01
20 -3.5011297041352050e+00 -3.7248508748877587e+00 2.3615662576274494e+00
21 -1.5332952658285048e+00 5.9630208068632040e-01 -7.4967230017303281e+00
22 4.2380253233105529e+00 1.0270453290850614e+00 6.6489894421385651e+00
23 -2.7047300574820481e+00 -1.6233474097713818e+00 8.4773355959176278e-01
24 -6.6588083188726532e+00 3.5110922792825918e+00 -6.5625174267043489e+00
25 7.9844426562464141e+00 -1.2853795683286129e+00 6.7123710742192300e+00
26 -1.3256343373737607e+00 -2.2257127109539789e+00 -1.4985364751488087e-01
27 6.6999960289138851e+00 6.3808952243186141e+00 2.0100808779497248e+00
28 -8.8466157439236681e-01 3.8018717064230995e-01 -5.9857060538593476e-01
29 -5.8153344545215182e+00 -6.7610823949609244e+00 -1.4115102725637900e+00
run_energy: 37.78424389351509
run_stress: ! |-
-4.6127506998693484e+01 -9.2129732247211749e+01 3.4548310342284810e+02 -2.9841348469661163e+01 -7.8434962689387717e+01 5.9253167412123155e+01
run_forces: ! |2
1 -5.8451208652159004e+00 -3.7483084455000643e+01 -2.7706576989352534e+01
2 -9.4646964278974774e+00 -7.8058897724822449e+00 1.1098831256058579e+01
3 3.1827086102630346e+01 2.7573911030624821e+01 -6.9576662575837211e+00
4 5.1502169659901655e+00 -1.4367546726785101e+00 3.6631301025186187e+00
5 -1.2208420775139264e+00 -1.8781699435112362e+00 6.2332639085051911e+00
6 -3.8491523409043303e+01 6.8063273218541468e+01 1.9723141045830272e+02
7 7.4838209349394775e+00 -3.8394258853636330e+01 -1.9092625515909930e+02
8 1.2676329319901857e+01 -1.2475162287097550e+01 3.3659783337736577e+00
9 6.8845241565874460e-01 7.3814593866184031e-01 3.0434095400342533e+00
10 6.2545583994797553e+00 -2.9600470917047201e+00 9.4247125735981765e+00
11 -1.9554747834212524e-01 -4.8434314068172696e-01 -7.9452259566032057e-01
12 5.2092795750960841e+00 3.1431929551776721e+00 -3.1346654851373348e+00
13 -1.1496483840617872e+01 4.5245217971580018e+00 2.1348220240918236e-01
14 3.1913399826660909e+00 -6.3760720126489068e-01 8.2740980433927742e+00
15 2.7338564489784484e-01 -9.7206665011069671e+00 -3.4841809697094543e+00
16 5.2461611410918316e+01 -3.1639255494702798e+01 -1.6483607587596811e+02
17 -5.8501866653548078e+01 4.0872194473703807e+01 1.5529162691391761e+02
18 -7.0990354207248405e+00 -2.4743922643289666e+00 -1.7824398936159682e+01
19 1.2019842510974870e+01 7.7105128268768715e+00 1.4523712108141252e+01
20 -4.9208070902500296e+00 -5.2361205625479048e+00 3.3006868280184283e+00
21 -1.8548628650934149e+00 2.7467524264262122e-01 -6.7601469408617412e+00
22 3.9136757840663186e+00 9.5561415744904055e-01 6.1181929861632272e+00
23 -2.0588129189729036e+00 -1.2302894000916618e+00 6.4195395469851357e-01
24 -5.7681973234153086e+00 2.0209144998436366e+00 -5.2864044021513967e+00
25 6.3696975292216704e+00 -1.0109756418053095e+00 5.3564043759405795e+00
26 -6.0150020580636188e-01 -1.0099388580383271e+00 -6.9999973789182365e-02
27 6.8467535469188450e+00 5.7500299184200578e+00 2.2775780974490298e+00
28 -1.3929430925479587e+00 5.9772788540443345e-01 -9.4056106886485980e-01
29 -5.4538104543708865e+00 -6.3477578038244911e+00 -1.3370170285841700e+00
...

83
unittest/force-styles/tests/mol-pair-lepton.yaml
View File

@ -1,6 +1,6 @@
---
lammps_version: 22 Dec 2022
date_generated: Thu Dec 22 09:57:30 2022
lammps_version: 21 Nov 2023
date_generated: Thu Jan 18 11:01:50 2024
epsilon: 5e-14
skip_tests: intel
prerequisites: ! |
@ -23,23 +23,24 @@ pair_coeff: ! |
2 4 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.005;sig=0.5"
2 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.00866025;sig=2.05"
3 3 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.02;sig=3.2"
3 4 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.0173205;sig=3.15"
3 4 "-eps*r;eps=0.0173205;sig=3.15"
3 5 "4.0*eps*((sig/r)^12-(sig/r)^6);eps=0.0173205;sig=3.15"
4 4 "10.0"
extract: ! ""
natoms: 29
init_vdwl: 749.2468149791969
init_vdwl: 746.1575578155301
init_coul: 0
init_stress: ! |2-
2.1793853434038242e+03 2.1988955172192768e+03 4.6653977523326257e+03 -7.5956547636050584e+02 2.4751536734032861e+01 6.6652028436400667e+02
2.1723526811665593e+03 2.1959162890293533e+03 4.6328064825512138e+03 -7.5509180369489252e+02 9.4506578600439983e+00 6.7585028859193505e+02
init_forces: ! |2
1 -2.3333390280895912e+01 2.6994567613322641e+02 3.3272827850356805e+02
1 -2.3359983837422618e+01 2.6996030011590727e+02 3.3274783233743295e+02
2 1.5828554630414899e+02 1.3025008843535872e+02 -1.8629682358935722e+02
3 -1.3528903738169066e+02 -3.8704313358319990e+02 -1.4568978437133106e+02
4 -7.8711096705893366e+00 2.1350518625373538e+00 -5.5954532185548134e+00
5 -2.5176757268228540e+00 -4.0521510681020239e+00 1.2152704057877019e+01
6 -8.3190662465252137e+02 9.6394149462625603e+02 1.1509093566509248e+03
7 5.8203388932513583e+01 -3.3608997951626793e+02 -1.7179617996573040e+03
8 1.4451392284291535e+02 -1.0927475861088995e+02 3.9990593492420442e+02
7 6.6340523101244187e+01 -3.4078810185436379e+02 -1.7003039516942540e+03
8 1.3674478037618434e+02 -1.0517874373121482e+02 3.8291074246191346e+02
9 7.9156945283097443e+01 8.5273009783986538e+01 3.5032175698445189e+02
10 5.3118875219105360e+02 -6.1040990859419412e+02 -1.8355872642619292e+02
11 -2.3530157267965532e+00 -5.9077640073819717e+00 -9.6590723955414290e+00
@ -48,8 +49,8 @@ init_forces: ! |2
14 -3.3852721292265153e+00 6.8636181241903649e-01 -8.7507190862499868e+00
15 -2.0454999188605300e-01 8.4846165523049883e+00 3.0131615419406712e+00
16 4.6326310311812108e+02 -3.3087715736498188e+02 -1.1893024561782554e+03
17 -4.5334300923766727e+02 3.1554283255882569e+02 1.2058417793481203e+03
18 -1.8862623280672661e-02 -3.3402010907951661e-02 3.1000479299095260e-02
17 -4.5371128972368928e+02 3.1609940794953951e+02 1.2052011419527653e+03
18 8.0197172683943874e-03 -2.4939258820032362e-03 -1.0571459969936936e-02
19 3.1843079640570047e-04 -2.3918627818763426e-04 1.7427252638513439e-03
20 -9.9760831209706009e-04 -1.0209184826753090e-03 3.6910972636601454e-04
21 -7.1566125273265186e+01 -8.1615678329920655e+01 2.2589561408339890e+02
@ -61,38 +62,38 @@ init_forces: ! |2
27 5.1810388677546001e+01 -2.2705458321213797e+02 9.0849111082069669e+01
28 -1.8041307121444069e+02 7.7534042932772905e+01 -1.2206956760706598e+02
29 1.2861057254925012e+02 1.4952711274394568e+02 3.1216025556267880e+01
run_vdwl: 719.4530651193046
run_vdwl: 716.5213000416621
run_coul: 0
run_stress: ! |2-
2.1330153957371017e+03 2.1547728168285516e+03 4.3976497417710125e+03 -7.3873328448298525e+02 4.1743821105370067e+01 6.2788012209191027e+02
2.1263870112744726e+03 2.1520080341389726e+03 4.3663519512361027e+03 -7.3456213833770062e+02 2.6927285459244832e+01 6.3691834104928068e+02
run_forces: ! |2
1 -2.0299419751359164e+01 2.6686193378823020e+02 3.2358785870694015e+02
2 1.5298617928491225e+02 1.2596516341409203e+02 -1.7961292655338619e+02
3 -1.3353630652439830e+02 -3.7923748696131315e+02 -1.4291839793625817e+02
4 -7.8374717836161762e+00 2.1276610789823409e+00 -5.5845014473820616e+00
5 -2.5014258630866735e+00 -4.0250131424704412e+00 1.2103512372025639e+01
6 -8.0681462887292457e+02 9.2165637136761688e+02 1.0270795806932783e+03
7 5.5780279349903523e+01 -3.1117530951561656e+02 -1.5746991292869018e+03
8 1.3452983055535049e+02 -1.0064659350255911e+02 3.8851791558207651e+02
9 7.6746213883425980e+01 8.2501469877402130e+01 3.3944351200617882e+02
10 5.2128033527695595e+02 -5.9920098848285863e+02 -1.8126029815043339e+02
11 -2.3573118090915246e+00 -5.8616944550888359e+00 -9.6049808811326205e+00
12 1.7503975847822900e+01 1.0626930310560814e+01 -8.0603160272054968e+00
13 8.0530313322973104e+00 -3.1756495170399117e+00 -1.4618315664740528e-01
14 -3.3416065168069773e+00 6.6492606336082150e-01 -8.6345131440469700e+00
15 -2.2253843262374914e-01 8.5025661635348779e+00 3.0369735873081622e+00
16 4.3476311264989465e+02 -3.1171086735551415e+02 -1.1135217194927448e+03
17 -4.2469846140777133e+02 2.9615411776780593e+02 1.1302573488400669e+03
18 -1.8849981672825908e-02 -3.3371636477421307e-02 3.0986293443778727e-02
19 3.0940277774414027e-04 -2.4634536455373044e-04 1.7433360008861016e-03
20 -9.8648131277150790e-04 -1.0112587134526946e-03 3.6932948773965417e-04
21 -7.0490745283106378e+01 -7.9749153581142139e+01 2.2171003384646431e+02
22 -1.0638717908920071e+02 -2.5949502163177968e+01 -1.6645589526812276e+02
23 1.7686797710735027e+02 1.0571018898885514e+02 -5.5243337084099387e+01
24 3.8206017656281375e+01 -2.1022820141992960e+02 1.1260711266189014e+02
25 -1.4918881473530880e+02 2.3762151395876508e+01 -1.2549188139143085e+02
26 1.1097059498808308e+02 1.8645503634228518e+02 1.2861559677865248e+01
27 5.0800844984832125e+01 -2.2296588090685469e+02 8.8607367716323253e+01
28 -1.7694190504288886e+02 7.6029945485182026e+01 -1.1950518150242071e+02
29 1.2614894925528141e+02 1.4694250820033548e+02 3.0893386672863034e+01
1 -2.0326040164905073e+01 2.6687684422507328e+02 3.2360752654223910e+02
2 1.5298608857690186e+02 1.2596506573447739e+02 -1.7961281277841888e+02
3 -1.3353631293077220e+02 -3.7923732277833739e+02 -1.4291833260989750e+02
4 -7.8374717116975035e+00 2.1276610267113969e+00 -5.5845014524498486e+00
5 -2.5014258756924157e+00 -4.0250131713717776e+00 1.2103512280982228e+01
6 -8.0714971444536457e+02 9.2203068890526424e+02 1.0274502514782534e+03
7 6.3722543724608350e+01 -3.1586173092061807e+02 -1.5580743968587681e+03
8 1.2737293861904031e+02 -9.6945064279519002e+01 3.7231518354375891e+02
9 7.6709940036396304e+01 8.2451980339096536e+01 3.3926849385746954e+02
10 5.2123408713149831e+02 -5.9914309504622599e+02 -1.8121478407355445e+02
11 -2.3573086824741427e+00 -5.8616969504300931e+00 -9.6049799947287671e+00
12 1.7504108236707797e+01 1.0626901299509713e+01 -8.0602444903747301e+00
13 8.0530313558451159e+00 -3.1756495145404533e+00 -1.4618321144421534e-01
14 -3.3416062225209915e+00 6.6492609500227240e-01 -8.6345136470911594e+00
15 -2.2253820242887132e-01 8.5025660110994483e+00 3.0369741645942137e+00
16 4.3476708820318731e+02 -3.1171425443331651e+02 -1.1135289618967258e+03
17 -4.2507048343681140e+02 2.9671384825884064e+02 1.1296230654445915e+03
18 8.0130752607770750e-03 -2.4895867517657545e-03 -1.0574351684568857e-02
19 3.0939970262803125e-04 -2.4635874092791046e-04 1.7433490521479268e-03
20 -9.8648319666298735e-04 -1.0112621691758337e-03 3.6933139856766442e-04
21 -7.0490745298133859e+01 -7.9749153568373742e+01 2.2171003384665224e+02
22 -1.0638717908973166e+02 -2.5949502162671845e+01 -1.6645589526807785e+02
23 1.7686797710711278e+02 1.0571018898899243e+02 -5.5243337084327727e+01
24 3.8206017659583978e+01 -2.1022820135505594e+02 1.1260711269986750e+02
25 -1.4918881473631544e+02 2.3762151403215309e+01 -1.2549188138812220e+02
26 1.1097059498835199e+02 1.8645503634383900e+02 1.2861559678659969e+01
27 5.0800844960383969e+01 -2.2296588092255456e+02 8.8607367714616288e+01
28 -1.7694190504410764e+02 7.6029945484553380e+01 -1.1950518150262033e+02
29 1.2614894924957088e+02 1.4694250819500266e+02 3.0893386676150566e+01
...

35
unittest/utils/test_lepton.cpp
View File

@ -542,6 +542,41 @@ TEST(Lepton, Optimize)
out.str("");
}
TEST(Lepton, Exception)
{
Lepton::CompiledExpression function, derivative;
auto parsed = Lepton::Parser::parse("x*x");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
double x = 1.5;
EXPECT_NO_THROW(function.getVariableReference("x") = x;);
EXPECT_NO_THROW(derivative.getVariableReference("x") = x;);
EXPECT_DOUBLE_EQ(function.evaluate(), 2.25);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 3.0);
parsed = Lepton::Parser::parse("x");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
x = 2.5;
EXPECT_NO_THROW(function.getVariableReference("x") = x;);
EXPECT_THROW(derivative.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_DOUBLE_EQ(function.evaluate(), 2.5);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 1.0);
parsed = Lepton::Parser::parse("1.0");
function = parsed.createCompiledExpression();
derivative = parsed.differentiate("x").createCompiledExpression();
x = 0.5;
EXPECT_THROW(function.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_THROW(derivative.getVariableReference("x") = x;, Lepton::Exception);
EXPECT_DOUBLE_EQ(function.evaluate(), 1.0);
EXPECT_DOUBLE_EQ(derivative.evaluate(), 0.0);
}
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);