diff --git a/doc/src/pair_lepton.rst b/doc/src/pair_lepton.rst index 303bc13bb9..34e7b20107 100644 --- a/doc/src/pair_lepton.rst +++ b/doc/src/pair_lepton.rst @@ -37,6 +37,8 @@ Examples pair_coeff * * "k*((r-r0)^2*step(r0-r)); k=200; r0=1.5" 2.0 pair_coeff 1 2 "4.0*eps*((sig/r)^12 - (sig/r)^6);eps=1.0;sig=1.0" 1.12246204830937 pair_coeff 2 2 "eps*(2.0*(sig/r)^9 - 3.0*(sig/r)^6);eps=1.0;sig=1.0" + pair_coeff 1 3 "zbl(13,6,r)" + pair_coeff 3 3 "(1.0-switch)*zbl(6,6,r)-switch*4.0*eps*((sig/r)^6);switch=0.5*(tanh(10.0*(r-sig))+1.0);eps=0.05;sig=3.20723" pair_style lepton/coul 2.5 pair_coeff 1 1 "qi*qj/r" 4.0 @@ -103,6 +105,15 @@ For pair style *lepton/coul* only the "coul" value of the :doc:`special_bonds `. The arguments of the function are the atomic numbers of +atom i (zi), atom j (zj) and the distance r. Please see the examples +above. + ---------- .. include:: lepton_expression.rst diff --git a/src/LEPTON/lepton_utils.cpp b/src/LEPTON/lepton_utils.cpp index d0cfa576cb..89e69beddd 100644 --- a/src/LEPTON/lepton_utils.cpp +++ b/src/LEPTON/lepton_utils.cpp @@ -20,14 +20,73 @@ #include "error.h" #include "input.h" #include "lammps.h" +#include "pair_zbl_const.h" #include "variable.h" #include "fmt/args.h" #include +#include #include #include +using namespace LAMMPS_NS; +using namespace PairZBLConstants; + +namespace Lepton { +class DerivativeException : public std::exception { + std::string message; + + public: + // remove unused default constructor + DerivativeException() = delete; + + explicit DerivativeException(int deg, const std::string &fn, const std::string &vn) + { + message = fmt::format("Order {} derivative of function {} in {} is not supported", deg, fn, vn); + } + const char *what() const noexcept override { return message.c_str(); } +}; +} // namespace Lepton + +double Lepton::ZBLFunction::evaluate(const double *args) const +{ + const double zi = args[0]; + const double zj = args[1]; + const double r = args[2]; + + const double rbya = r * (pow(zi, pzbl) + pow(zj, pzbl)) / (a0 * angstrom); + return zi * zj * qqr2e * qelectron * qelectron / r * + (c4 * exp(-d4 * rbya) + c3 * exp(-d3 * rbya) + c2 * exp(-d2 * rbya) + c1 * exp(-d1 * rbya)); +} + +double Lepton::ZBLFunction::evaluateDerivative(const double *args, const int *order) const +{ + if (order[0] > 0) + throw DerivativeException(order[0], "zbl()", "'zi'"); + if (order[1] > 0) + throw DerivativeException(order[0], "zbl()", "'zj'"); + if (order[2] > 1) + throw DerivativeException(order[0], "zbl()", "'r'"); + + if (order[2] == 1) { + const double zi = args[0]; + const double zj = args[1]; + const double r = args[2]; + + const double ainv = (pow(zi, pzbl) + pow(zj, pzbl)) / (a0 * angstrom); + const double e1 = exp(-d1 * ainv * r); + const double e2 = exp(-d2 * ainv * r); + const double e3 = exp(-d3 * ainv * r); + const double e4 = exp(-d4 * ainv * r); + + const double sum1 = c1 * e1 + c2 * e2 + c3 * e3 + c4 * e4; + const double sum2 = ainv * (-c1 * d1 * e1 - c2 * d2 * e2 - c3 * d3 * e3 - c4 * d4 * e4); + return (zi * zj * qqr2e * qelectron * qelectron) * (sum2 - sum1 / r) / r; + } + return 0.0; +} + namespace LeptonUtils { class VariableException : public std::exception { std::string message; diff --git a/src/LEPTON/lepton_utils.h b/src/LEPTON/lepton_utils.h index 89512e7552..caa1ce56ab 100644 --- a/src/LEPTON/lepton_utils.h +++ b/src/LEPTON/lepton_utils.h @@ -11,6 +11,8 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "Lepton.h" + #include // forward declarations @@ -19,6 +21,25 @@ namespace LAMMPS_NS { class LAMMPS; } +// custom function zbl(zi,zj,r) + +namespace Lepton { +class ZBLFunction : public CustomFunction { + public: + ZBLFunction(double _qqr2e, double _angstrom, double _qelectron) : + qqr2e(_qqr2e), angstrom(_angstrom), qelectron(_qelectron){}; + ZBLFunction() = delete; + + int getNumArguments() const override { return 3; } + CustomFunction *clone() const override { return new ZBLFunction(qqr2e, angstrom, qelectron); } + double evaluate(const double *) const override; + double evaluateDerivative(const double *, const int *) const override; + + private: + double qqr2e, angstrom, qelectron; +}; +} // namespace Lepton + // utility functions and classes namespace LeptonUtils { diff --git a/src/LEPTON/pair_lepton.cpp b/src/LEPTON/pair_lepton.cpp index f0ccebbb44..84049ea291 100644 --- a/src/LEPTON/pair_lepton.cpp +++ b/src/LEPTON/pair_lepton.cpp @@ -28,6 +28,7 @@ #include "Lepton.h" #include "lepton_utils.h" #include +#include using namespace LAMMPS_NS; @@ -43,12 +44,15 @@ PairLepton::PairLepton(LAMMPS *lmp) : reinitflag = 0; cut_global = 0.0; centroidstressflag = CENTROID_SAME; + + functions["zbl"] = new Lepton::ZBLFunction(force->qqr2e, force->angstrom, force->qelectron); } /* ---------------------------------------------------------------------- */ PairLepton::~PairLepton() { + for (auto &f : functions) delete f.second; if (allocated) { memory->destroy(cut); memory->destroy(cutsq); @@ -103,7 +107,7 @@ template void PairLepton::eval() std::vector pairpot; try { for (const auto &expr : expressions) { - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); pairforce.back().getVariableReference("r"); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); @@ -222,7 +226,7 @@ void PairLepton::coeff(int narg, char **arg) // check if the expression can be parsed and evaluated without error auto exp_one = LeptonUtils::condense(arg[2]); try { - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp)); + 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; @@ -265,7 +269,7 @@ double PairLepton::init_one(int i, int j) if (offset_flag) { try { auto expr = LeptonUtils::substitute(expressions[type2expression[i][j]], lmp); - auto pairpot = Lepton::Parser::parse(expr).createCompiledExpression(); + auto pairpot = Lepton::Parser::parse(expr, functions).createCompiledExpression(); pairpot.getVariableReference("r") = cut[i][j]; offset[i][j] = pairpot.evaluate(); } catch (std::exception &) { @@ -420,7 +424,7 @@ double PairLepton::single(int /* i */, int /* j */, int itype, int jtype, double double /* factor_coul */, double factor_lj, double &fforce) { auto expr = expressions[type2expression[itype][jtype]]; - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); auto pairpot = parsed.createCompiledExpression(); auto pairforce = parsed.differentiate("r").createCompiledExpression(); diff --git a/src/LEPTON/pair_lepton.h b/src/LEPTON/pair_lepton.h index e8454ce80e..8ec03bf28b 100644 --- a/src/LEPTON/pair_lepton.h +++ b/src/LEPTON/pair_lepton.h @@ -22,6 +22,12 @@ PairStyle(lepton,PairLepton); #include "pair.h" +#include + +namespace Lepton { +class CustomFunction; +} + namespace LAMMPS_NS { class PairLepton : public Pair { @@ -42,6 +48,8 @@ class PairLepton : public Pair { protected: std::vector expressions; + std::map functions; + double **cut; int **type2expression; double **offset; diff --git a/src/LEPTON/pair_lepton_coul.cpp b/src/LEPTON/pair_lepton_coul.cpp index 81d2ff7a81..863daf57c6 100644 --- a/src/LEPTON/pair_lepton_coul.cpp +++ b/src/LEPTON/pair_lepton_coul.cpp @@ -82,7 +82,7 @@ template void PairLeptonCoul::eval() std::vector> have_q; try { for (const auto &expr : expressions) { - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); pairforce.back().getVariableReference("r"); @@ -243,7 +243,7 @@ double PairLeptonCoul::single(int i, int j, int itype, int jtype, double rsq, do double /* factor_lj */, double &fforce) { auto expr = expressions[type2expression[itype][jtype]]; - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); auto pairpot = parsed.createCompiledExpression(); auto pairforce = parsed.differentiate("r").createCompiledExpression(); diff --git a/src/OPENMP/pair_lepton_coul_omp.cpp b/src/OPENMP/pair_lepton_coul_omp.cpp index a67779aeb8..14d49c308f 100644 --- a/src/OPENMP/pair_lepton_coul_omp.cpp +++ b/src/OPENMP/pair_lepton_coul_omp.cpp @@ -104,7 +104,7 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr) std::vector> have_q; try { for (const auto &expr : expressions) { - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); pairforce.back().getVariableReference("r"); diff --git a/src/OPENMP/pair_lepton_omp.cpp b/src/OPENMP/pair_lepton_omp.cpp index 2c96b63a7f..b57b0fe11e 100644 --- a/src/OPENMP/pair_lepton_omp.cpp +++ b/src/OPENMP/pair_lepton_omp.cpp @@ -98,7 +98,7 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr) std::vector pairpot; try { for (const auto &expr : expressions) { - auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp)); + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions); pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); } diff --git a/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml b/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml new file mode 100644 index 0000000000..883021c26b --- /dev/null +++ b/unittest/force-styles/tests/mol-pair-lepton_zbl.yaml @@ -0,0 +1,102 @@ +--- +lammps_version: 22 Dec 2022 +date_generated: Sun Jan 8 01:13:49 2023 +epsilon: 5e-10 +skip_tests: +prerequisites: ! | + atom full + pair lepton +pre_commands: ! | + variable bond_factor index 1.0 + variable angle_factor index 1.0 + variable dihedral_factor index 1.0 +post_commands: ! "" +input_file: in.fourmol +pair_style: lepton 8.0 +pair_coeff: ! | + 1 1 zbl(6,6,r) + 1 2 zbl(6,7,r) + 1 3 zbl(6,2,r) + 1 4 zbl(6,8,r) + 1 5 zbl(6,8,r) + 2 2 zbl(7,7,r) + 2 3 zbl(7,2,r) + 2 4 zbl(7,8,r) + 2 5 zbl(7,8,r) + 3 3 zbl(2,2,r) + 3 4 zbl(2,8,r) + 3 5 zbl(2,8,r) + 4 4 zbl(8,8,r) + 4 5 zbl(8,8,r) + 5 5 zbl(8,8,r) +extract: ! "" +natoms: 29 +init_vdwl: 7661.0923034839425 +init_coul: 0 +init_stress: ! |2- + 1.0415723081438835e+04 7.6654640475740289e+03 1.0439568803573100e+04 2.5596340687548091e+03 3.0124073227376857e+03 3.1819159178174647e+02 +init_forces: ! |2 + 1 -2.2542878982886191e+02 -5.1056064396745825e+01 4.0347728160298351e+02 + 2 3.4577199265194707e+02 3.3685870288327618e+02 -3.5822722442776552e+02 + 3 6.7803530899610132e+02 3.5041867763592933e+01 -3.5973378440582275e+02 + 4 -8.1435522635968925e+02 2.2040390876529949e+02 -5.9600610062358885e+02 + 5 -2.1098162371264527e+02 -3.3211428613105602e+02 1.0436615327229824e+03 + 6 -8.2893212979873752e+01 1.2663096062116030e+02 4.4419304268422735e+02 + 7 -3.4934572486178489e+01 -1.2862622365409325e+02 -7.5679463440025643e+02 + 8 -1.4121290372363833e+02 -1.1033884290619251e+02 -3.6118050221036594e+02 + 9 9.6547187735733928e+01 2.0268265922256009e+02 6.9609964660581829e+02 + 10 8.2081753362417444e+01 2.1490541800870224e+02 6.7445631928792000e+02 + 11 -2.7501706654968336e+02 -6.4133882479619604e+02 -8.2554240753015847e+02 + 12 -2.3615126531081731e+02 -3.7888094263940093e+02 4.2996552265607966e+02 + 13 9.8935659849196156e+02 -3.9189784587066578e+02 1.6979303593289488e+00 + 14 -3.2831677005156325e+02 9.7498291303146885e+01 -1.0618689865846111e+03 + 15 2.3059516739411034e+01 1.0120552702611462e+03 3.8891358018891088e+02 + 16 3.2227110665984651e+02 -3.2696296221307608e+02 -3.9821264616502327e+02 + 17 -1.9854928611614176e+02 1.0231604600726870e+02 6.4707382129491782e+02 + 18 -1.1566833968062178e+02 -5.0500912216600113e+02 1.7898146720451105e+03 + 19 -1.0355744892944338e+03 -6.8474702051102486e+02 -1.1421521343728634e+03 + 20 1.1567664608091418e+03 1.2055463642684399e+03 -6.6488505245382839e+02 + 21 -5.1058201416833657e+02 -6.2261916138266497e+02 1.8034425299414422e+03 + 22 -1.0303537305440429e+03 -2.8540172556947783e+02 -1.4390471053430481e+03 + 23 1.5447065970289084e+03 9.0276100147850616e+02 -3.6857879996771908e+02 + 24 3.7449763655418099e+02 -1.6372799599084626e+03 9.3993003918438728e+02 + 25 -1.3650956291199398e+03 1.2266826749083317e+02 -1.1043020934470092e+03 + 26 9.9428586164565161e+02 1.5187270928735547e+03 1.7246869495941340e+02 + 27 2.8070997771034445e+02 -1.7944696249727249e+03 6.4639058889384091e+02 + 28 -1.4822994833123780e+03 5.2787642883033504e+02 -9.7293790834683261e+02 + 29 1.1993244048531997e+03 1.2647703273399616e+03 3.2788417785152956e+02 +run_vdwl: 7355.7708040367215 +run_coul: 0 +run_stress: ! |2- + 1.0014026118319458e+04 7.4146110545302845e+03 1.0129433554690771e+04 2.4014185997980703e+03 2.8618718144358500e+03 3.1547123090533302e+02 +run_forces: ! |2 + 1 -2.1963909767996165e+02 -4.7230048342286665e+01 3.9447474330057327e+02 + 2 3.3823602012753361e+02 3.3060841467248804e+02 -3.4909191972075922e+02 + 3 6.4685420887736268e+02 3.4155483339278817e+01 -3.4174127433776528e+02 + 4 -7.8636648191524102e+02 2.1116614600799022e+02 -5.7402164047826784e+02 + 5 -2.0592655691757219e+02 -3.1948425705952934e+02 1.0031851850554129e+03 + 6 -8.1885702477873082e+01 1.2476420007196047e+02 4.3622385152526408e+02 + 7 -3.5214786990923372e+01 -1.2723182110324109e+02 -7.4725777446879317e+02 + 8 -1.3797341261617603e+02 -1.0667961036533134e+02 -3.4719859501013519e+02 + 9 9.3157498402761519e+01 1.9853026083255594e+02 6.7970067058747793e+02 + 10 7.4694873724181420e+01 1.9549844828167281e+02 6.3962599563799461e+02 + 11 -2.6505567174539726e+02 -6.1810691652881007e+02 -7.8981506764583025e+02 + 12 -2.2346480012776817e+02 -3.6198854997039064e+02 4.0604938064590891e+02 + 13 9.5587946001324235e+02 -3.7492976806829023e+02 4.0412877209292042e-01 + 14 -3.1069307401641237e+02 9.5118496292191480e+01 -1.0209459166022623e+03 + 15 2.3074811498468790e+01 9.7760720729835271e+02 3.7462072751151828e+02 + 16 3.2004057723494378e+02 -3.2580441129415715e+02 -3.9361639626975176e+02 + 17 -1.9645839108898824e+02 1.0114353674483996e+02 6.4131023697143405e+02 + 18 -9.0905016804868524e+01 -4.5649801918985156e+02 1.6740121681099276e+03 + 19 -9.7189809015843355e+02 -6.4245306822849329e+02 -1.0761208024504872e+03 + 20 1.0683309313831251e+03 1.1147511786417067e+03 -6.1516522889263013e+02 + 21 -4.7823440562179218e+02 -5.8004068421284194e+02 1.6903616731086115e+03 + 22 -9.6232775232513484e+02 -2.6660807378911699e+02 -1.3469354884213956e+03 + 23 1.4443362128008671e+03 8.4138206482572116e+02 -3.4761279951083094e+02 + 24 3.5961673762879963e+02 -1.5399088881652940e+03 8.8844237731370015e+02 + 25 -1.2856259100984250e+03 1.1769150558030140e+02 -1.0419813334649532e+03 + 26 9.2973443424366496e+02 1.4263851950698765e+03 1.6174119101406066e+02 + 27 2.6191415819111842e+02 -1.6858783749946833e+03 6.0207717487186846e+02 + 28 -1.3837052677973970e+03 4.9654194958937097e+02 -9.0783879641816793e+02 + 29 1.1195044942562959e+03 1.1874984040640106e+03 3.0711352926618406e+02 +... diff --git a/unittest/utils/test_lepton.cpp b/unittest/utils/test_lepton.cpp index 91532b385a..415752c70d 100644 --- a/unittest/utils/test_lepton.cpp +++ b/unittest/utils/test_lepton.cpp @@ -85,6 +85,47 @@ TEST_F(LeptonUtilsTest, substitute) } } +// zbl() custom function + +TEST(LeptonCustomFunction, zbl) +{ + Lepton::ZBLFunction zbl(1.0, 1.0, 1.0); + std::map functions = {std::make_pair("zbl", &zbl)}; + std::map variables = {std::make_pair("zi", 6), std::make_pair("zj", 6), + std::make_pair("r", 2.0)}; + + auto parsed = Lepton::Parser::parse("zbl(zi, zj, r)", functions); + auto zbldzi = parsed.differentiate("zi"); + auto zbldzj = parsed.differentiate("zj"); + auto zbldr = parsed.differentiate("r"); + auto zbld2r = zbldr.differentiate("r"); + + double value = parsed.evaluate(variables); + ASSERT_DOUBLE_EQ(value, 0.065721538245489763); + value = zbldr.evaluate(variables); + ASSERT_DOUBLE_EQ(value, -0.15481915325334394); + variables["r"] = 1.0; + value = parsed.evaluate(variables); + ASSERT_DOUBLE_EQ(value, 1.0701488641432269); + value = zbldr.evaluate(variables); + ASSERT_DOUBLE_EQ(value, -3.6376386525054412); + variables["zi"] = 13.0; + value = parsed.evaluate(variables); + ASSERT_DOUBLE_EQ(value, 1.8430432789454971); + value = zbldr.evaluate(variables); + ASSERT_DOUBLE_EQ(value, -6.5373118484557642); + variables["zj"] = 13.0; + value = parsed.evaluate(variables); + ASSERT_DOUBLE_EQ(value, 3.1965196467438446); + value = zbldr.evaluate(variables); + ASSERT_DOUBLE_EQ(value, -11.804490148948526); + + // check for unsupported derivatives + ASSERT_ANY_THROW(value = zbldzi.evaluate(variables)); + ASSERT_ANY_THROW(value = zbldzj.evaluate(variables)); + ASSERT_ANY_THROW(value = zbld2r.evaluate(variables)); +} + /** * This is a custom function equal to f(x,y) = 2*x*y. */