From 5304c43fef26d1f113eb37a4971518551ed2e14e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 18 Jan 2024 04:32:20 -0500 Subject: [PATCH 1/5] add test for exceptions and evaluation of compiled expressions --- unittest/utils/test_lepton.cpp | 35 ++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/unittest/utils/test_lepton.cpp b/unittest/utils/test_lepton.cpp index a9fa6e3543..55d3bf8351 100644 --- a/unittest/utils/test_lepton.cpp +++ b/unittest/utils/test_lepton.cpp @@ -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); From 57db9be64fc578a97e15e84d48861d62f350d6d5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 18 Jan 2024 09:46:23 -0500 Subject: [PATCH 2/5] add support for Lepton expressions with constant force or constant potential --- src/LEPTON/angle_lepton.cpp | 49 +++++++++++++++++++----- src/LEPTON/bond_lepton.cpp | 37 +++++++++++++++--- src/LEPTON/dihedral_lepton.cpp | 29 ++++++++++++-- src/LEPTON/fix_wall_lepton.cpp | 15 +++++++- src/LEPTON/pair_lepton.cpp | 41 +++++++++++++++----- src/LEPTON/pair_lepton_coul.cpp | 50 +++++++++++++++--------- src/LEPTON/pair_lepton_sphere.cpp | 55 +++++++++++++++++---------- src/OPENMP/angle_lepton_omp.cpp | 16 ++++++-- src/OPENMP/bond_lepton_omp.cpp | 15 +++++++- src/OPENMP/dihedral_lepton_omp.cpp | 19 +++++++-- src/OPENMP/pair_lepton_coul_omp.cpp | 42 ++++++++++++-------- src/OPENMP/pair_lepton_omp.cpp | 20 ++++++++-- src/OPENMP/pair_lepton_sphere_omp.cpp | 44 +++++++++++++-------- 13 files changed, 319 insertions(+), 113 deletions(-) diff --git a/src/LEPTON/angle_lepton.cpp b/src/LEPTON/angle_lepton.cpp index 59310f5637..70116287d0 100644 --- a/src/LEPTON/angle_lepton.cpp +++ b/src/LEPTON/angle_lepton.cpp @@ -90,10 +90,21 @@ template void AngleLepton::eval() { std::vector angleforce; std::vector 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 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 +153,7 @@ template 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 +189,11 @@ template 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) @@ -224,8 +238,19 @@ 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; + 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); + } offset_one = anglepot.evaluate(); angleforce.evaluate(); } catch (std::exception &e) { @@ -363,7 +388,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]; } diff --git a/src/LEPTON/bond_lepton.cpp b/src/LEPTON/bond_lepton.cpp index 773607782d..119981fab6 100644 --- a/src/LEPTON/bond_lepton.cpp +++ b/src/LEPTON/bond_lepton.cpp @@ -82,10 +82,17 @@ template void BondLepton::eval() { std::vector bondforce; std::vector bondpot; + std::vector 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 +123,7 @@ template 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 +143,11 @@ template 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); @@ -179,8 +190,18 @@ 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; + 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); + } offset_one = bondpot.evaluate(); bondforce.evaluate(); } catch (std::exception &e) { @@ -302,8 +323,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 diff --git a/src/LEPTON/dihedral_lepton.cpp b/src/LEPTON/dihedral_lepton.cpp index 6470e43033..069ff13d74 100644 --- a/src/LEPTON/dihedral_lepton.cpp +++ b/src/LEPTON/dihedral_lepton.cpp @@ -92,10 +92,17 @@ template void DihedralLepton::eval() { std::vector dihedralforce; std::vector dihedralpot; + std::vector 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 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 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()); diff --git a/src/LEPTON/fix_wall_lepton.cpp b/src/LEPTON/fix_wall_lepton.cpp index a81d3c4edb..7530188c00 100644 --- a/src/LEPTON/fix_wall_lepton.cpp +++ b/src/LEPTON/fix_wall_lepton.cpp @@ -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) { diff --git a/src/LEPTON/pair_lepton.cpp b/src/LEPTON/pair_lepton.cpp index a8af0ce576..adc07cbfa8 100644 --- a/src/LEPTON/pair_lepton.cpp +++ b/src/LEPTON/pair_lepton.cpp @@ -27,6 +27,7 @@ #include "Lepton.h" #include "lepton_utils.h" +#include #include #include @@ -105,11 +106,17 @@ template void PairLepton::eval() std::vector pairforce; std::vector pairpot; + std::vector 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 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 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; } diff --git a/src/LEPTON/pair_lepton_coul.cpp b/src/LEPTON/pair_lepton_coul.cpp index 841565e874..f7d2042874 100644 --- a/src/LEPTON/pair_lepton_coul.cpp +++ b/src/LEPTON/pair_lepton_coul.cpp @@ -28,6 +28,8 @@ #include "Lepton.h" #include "lepton_utils.h" + +#include #include using namespace LAMMPS_NS; @@ -79,25 +81,30 @@ template void PairLeptonCoul::eval() std::vector pairforce; std::vector pairpot; - std::vector> have_q; + std::vector> 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 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 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 */ } diff --git a/src/LEPTON/pair_lepton_sphere.cpp b/src/LEPTON/pair_lepton_sphere.cpp index 29514aed38..72d0e85d0b 100644 --- a/src/LEPTON/pair_lepton_sphere.cpp +++ b/src/LEPTON/pair_lepton_sphere.cpp @@ -28,6 +28,7 @@ #include "Lepton.h" #include "lepton_utils.h" +#include #include using namespace LAMMPS_NS; @@ -77,25 +78,30 @@ template void PairLeptonSphere::eval() std::vector pairforce; std::vector pairpot; - std::vector> have_rad; + std::vector> 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 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 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; diff --git a/src/OPENMP/angle_lepton_omp.cpp b/src/OPENMP/angle_lepton_omp.cpp index 7e86a9e9bb..f57cf916a2 100644 --- a/src/OPENMP/angle_lepton_omp.cpp +++ b/src/OPENMP/angle_lepton_omp.cpp @@ -91,10 +91,17 @@ void AngleLeptonOMP::eval(int nfrom, int nto, ThrData *const thr) { std::vector angleforce; std::vector anglepot; + std::vector 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) diff --git a/src/OPENMP/bond_lepton_omp.cpp b/src/OPENMP/bond_lepton_omp.cpp index 0029062366..d9982b08f8 100644 --- a/src/OPENMP/bond_lepton_omp.cpp +++ b/src/OPENMP/bond_lepton_omp.cpp @@ -89,10 +89,17 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr) { std::vector bondforce; std::vector bondpot; + std::vector 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) diff --git a/src/OPENMP/dihedral_lepton_omp.cpp b/src/OPENMP/dihedral_lepton_omp.cpp index 13a1328058..37748ce9d5 100644 --- a/src/OPENMP/dihedral_lepton_omp.cpp +++ b/src/OPENMP/dihedral_lepton_omp.cpp @@ -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 @@ -94,10 +94,17 @@ void DihedralLeptonOMP::eval(int nfrom, int nto, ThrData *const thr) { std::vector dihedralforce; std::vector dihedralpot; + std::vector 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) diff --git a/src/OPENMP/pair_lepton_coul_omp.cpp b/src/OPENMP/pair_lepton_coul_omp.cpp index bc34bc00af..532c16d797 100644 --- a/src/OPENMP/pair_lepton_coul_omp.cpp +++ b/src/OPENMP/pair_lepton_coul_omp.cpp @@ -20,11 +20,13 @@ #include "neigh_list.h" #include "suffix.h" -#include - #include "Lepton.h" #include "lepton_utils.h" #include "omp_compat.h" + +#include +#include + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -101,25 +103,30 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr) std::vector pairforce; std::vector pairpot; - std::vector> have_q; + std::vector> 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; } diff --git a/src/OPENMP/pair_lepton_omp.cpp b/src/OPENMP/pair_lepton_omp.cpp index b57b0fe11e..58692e52d6 100644 --- a/src/OPENMP/pair_lepton_omp.cpp +++ b/src/OPENMP/pair_lepton_omp.cpp @@ -20,11 +20,12 @@ #include "neigh_list.h" #include "suffix.h" -#include - #include "Lepton.h" #include "lepton_utils.h" #include "omp_compat.h" +#include +#include + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -96,10 +97,17 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr) std::vector pairforce; std::vector pairpot; + std::vector 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; } diff --git a/src/OPENMP/pair_lepton_sphere_omp.cpp b/src/OPENMP/pair_lepton_sphere_omp.cpp index 6d3a4827b3..79afe27717 100644 --- a/src/OPENMP/pair_lepton_sphere_omp.cpp +++ b/src/OPENMP/pair_lepton_sphere_omp.cpp @@ -20,11 +20,13 @@ #include "neigh_list.h" #include "suffix.h" -#include - #include "Lepton.h" #include "lepton_utils.h" #include "omp_compat.h" + +#include +#include + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -99,25 +101,30 @@ void PairLeptonSphereOMP::eval(int iifrom, int iito, ThrData *const thr) std::vector pairforce; std::vector pairpot; - std::vector> have_rad; + std::vector> 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; } From 73194764e9b25a20607fadaa5e8b5a22f26cbe28 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 18 Jan 2024 11:05:52 -0500 Subject: [PATCH 3/5] add optional keywords "auto_offset" and "no_offset" to bond/angle style lepton --- doc/src/angle_lepton.rst | 19 +++- doc/src/bond_lepton.rst | 19 +++- src/LEPTON/angle_lepton.cpp | 25 +++++- src/LEPTON/angle_lepton.h | 2 + src/LEPTON/bond_lepton.cpp | 25 +++++- src/LEPTON/bond_lepton.h | 2 + unittest/force-styles/tests/angle-lepton.yaml | 2 +- .../tests/angle-lepton_nooffset.yaml | 88 ++++++++++++++++++ .../tests/bond-lepton_nooffset.yaml | 89 +++++++++++++++++++ 9 files changed, 266 insertions(+), 5 deletions(-) create mode 100644 unittest/force-styles/tests/angle-lepton_nooffset.yaml create mode 100644 unittest/force-styles/tests/bond-lepton_nooffset.yaml diff --git a/doc/src/angle_lepton.rst b/doc/src/angle_lepton.rst index 20fa5b1fee..c60975d740 100644 --- a/doc/src/angle_lepton.rst +++ b/doc/src/angle_lepton.rst @@ -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 he 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 `_, that the *lepton* angle style interfaces with, evaluates this expression string at run time to compute the pairwise energy. It also creates an diff --git a/doc/src/bond_lepton.rst b/doc/src/bond_lepton.rst index adfd30627d..9429535af8 100644 --- a/doc/src/bond_lepton.rst +++ b/doc/src/bond_lepton.rst @@ -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 `_, that the *lepton* bond style interfaces with, evaluates this expression string at run time to compute the pairwise energy. It also creates an analytical diff --git a/src/LEPTON/angle_lepton.cpp b/src/LEPTON/angle_lepton.cpp index 70116287d0..9fe565f8ee 100644 --- a/src/LEPTON/angle_lepton.cpp +++ b/src/LEPTON/angle_lepton.cpp @@ -44,6 +44,7 @@ AngleLepton::AngleLepton(LAMMPS *_lmp) : { writedata = 1; reinitflag = 0; + auto_offset = 1; } /* ---------------------------------------------------------------------- */ @@ -216,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 ------------------------------------------------------------------------- */ @@ -251,7 +270,7 @@ void AngleLepton::coeff(int narg, char **arg) error->warning(FLERR, "Force from Lepton expression {} does not depend on 'theta'", exp_one); } - offset_one = anglepot.evaluate(); + if (auto_offset) offset_one = anglepot.evaluate(); angleforce.evaluate(); } catch (std::exception &e) { error->all(FLERR, e.what()); @@ -309,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); } /* ---------------------------------------------------------------------- @@ -348,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; } diff --git a/src/LEPTON/angle_lepton.h b/src/LEPTON/angle_lepton.h index 67d2718fb6..4f0e5729ed 100644 --- a/src/LEPTON/angle_lepton.h +++ b/src/LEPTON/angle_lepton.h @@ -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(); diff --git a/src/LEPTON/bond_lepton.cpp b/src/LEPTON/bond_lepton.cpp index 119981fab6..8679d0ed62 100644 --- a/src/LEPTON/bond_lepton.cpp +++ b/src/LEPTON/bond_lepton.cpp @@ -37,6 +37,7 @@ BondLepton::BondLepton(LAMMPS *_lmp) : { writedata = 1; reinitflag = 0; + auto_offset = 1; } /* ---------------------------------------------------------------------- */ @@ -168,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 ------------------------------------------------------------------------- */ @@ -202,7 +221,7 @@ void BondLepton::coeff(int narg, char **arg) if (comm->me == 0) error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one); } - offset_one = bondpot.evaluate(); + if (auto_offset) offset_one = bondpot.evaluate(); bondforce.evaluate(); } catch (std::exception &e) { error->all(FLERR, e.what()); @@ -260,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); } /* ---------------------------------------------------------------------- @@ -299,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; } diff --git a/src/LEPTON/bond_lepton.h b/src/LEPTON/bond_lepton.h index 9e693298a7..e59648a3f0 100644 --- a/src/LEPTON/bond_lepton.h +++ b/src/LEPTON/bond_lepton.h @@ -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(); diff --git a/unittest/force-styles/tests/angle-lepton.yaml b/unittest/force-styles/tests/angle-lepton.yaml index ea108cfdb1..b4d6c0516f 100644 --- a/unittest/force-styles/tests/angle-lepton.yaml +++ b/unittest/force-styles/tests/angle-lepton.yaml @@ -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" diff --git a/unittest/force-styles/tests/angle-lepton_nooffset.yaml b/unittest/force-styles/tests/angle-lepton_nooffset.yaml new file mode 100644 index 0000000000..711f0cbdd5 --- /dev/null +++ b/unittest/force-styles/tests/angle-lepton_nooffset.yaml @@ -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 +... diff --git a/unittest/force-styles/tests/bond-lepton_nooffset.yaml b/unittest/force-styles/tests/bond-lepton_nooffset.yaml new file mode 100644 index 0000000000..b39288640e --- /dev/null +++ b/unittest/force-styles/tests/bond-lepton_nooffset.yaml @@ -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 +... From 90a79d9a4b819da7210b4df65f943d91273ab827 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 18 Jan 2024 11:06:52 -0500 Subject: [PATCH 4/5] change unit test to include expressions with constant force or potential --- .../force-styles/tests/mol-pair-lepton.yaml | 83 ++++++++++--------- 1 file changed, 42 insertions(+), 41 deletions(-) diff --git a/unittest/force-styles/tests/mol-pair-lepton.yaml b/unittest/force-styles/tests/mol-pair-lepton.yaml index 03117a9aa5..33576e81c2 100644 --- a/unittest/force-styles/tests/mol-pair-lepton.yaml +++ b/unittest/force-styles/tests/mol-pair-lepton.yaml @@ -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 ... From 425421c1ca52129c197dd487c6220dd3efcc952d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 18 Jan 2024 14:32:40 -0500 Subject: [PATCH 5/5] fix typo, cut-n-paste error, and add clarification --- doc/src/angle_lepton.rst | 2 +- doc/src/pair_lepton.rst | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/angle_lepton.rst b/doc/src/angle_lepton.rst index c60975d740..22873f5765 100644 --- a/doc/src/angle_lepton.rst +++ b/doc/src/angle_lepton.rst @@ -53,7 +53,7 @@ angle coefficient. For example `"200.0*theta^2"` represents a .. versionchanged:: TBD -By default the potential energy U is shifted so that he value U is 0.0 +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. diff --git a/doc/src/pair_lepton.rst b/doc/src/pair_lepton.rst index 21e619a3d9..5b5dc698e7 100644 --- a/doc/src/pair_lepton.rst +++ b/doc/src/pair_lepton.rst @@ -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 ` or the :doc:`set command `. @@ -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 ` -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 ` options are not relevant for the these pair styles.