diff --git a/src/LEPTON/bond_lepton.cpp b/src/LEPTON/bond_lepton.cpp new file mode 100644 index 0000000000..ad2ac0af12 --- /dev/null +++ b/src/LEPTON/bond_lepton.cpp @@ -0,0 +1,322 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "bond_lepton.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neighbor.h" + +#include +#include +#include + +#include "LMP_Lepton.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondLepton::BondLepton(LAMMPS *_lmp) : Bond(_lmp), r0(nullptr), type2expression(nullptr) +{ + writedata = 1; + reinitflag = 0; +} + +/* ---------------------------------------------------------------------- */ + +BondLepton::~BondLepton() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(r0); + memory->destroy(type2expression); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondLepton::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + ev_init(eflag, vflag); + if (evflag) { + if (eflag) { + if (force->newton_bond) + eval<1, 1, 1>(); + else + eval<1, 1, 0>(); + } else { + if (force->newton_bond) + eval<1, 0, 1>(); + else + eval<1, 0, 0>(); + } + } else { + if (force->newton_bond) + eval<0, 0, 1>(); + else + eval<0, 0, 0>(); + } +} + +/* ---------------------------------------------------------------------- */ +template void BondLepton::eval() +{ + std::vector bondforce; + std::vector bondpot; + for (const auto &expr : expressions) { + auto parsed = LMP_Lepton::Parser::parse(expr); + bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); + if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression()); + } + + const double *const *const x = atom->x; + double *const *const f = atom->f; + const int *const *const bondlist = neighbor->bondlist; + const int nbondlist = neighbor->nbondlist; + const int nlocal = atom->nlocal; + + for (int n = 0; n < nbondlist; n++) { + const int i1 = bondlist[n][0]; + const int i2 = bondlist[n][1]; + const int type = bondlist[n][2]; + + const double delx = x[i1][0] - x[i2][0]; + const double dely = x[i1][1] - x[i2][1]; + const double delz = x[i1][2] - x[i2][2]; + + const double rsq = delx * delx + dely * dely + delz * delz; + const double r = sqrt(rsq); + const double dr = r - r0[type]; + const int idx = type2expression[type]; + + // force and energy + + double fbond = 0.0; + if (r > 0.0) { + double &r_for = bondforce[idx].getVariableReference("r"); + r_for = dr; + fbond = -bondforce[idx].evaluate() / r; + } + + // apply force to each of 2 atoms + + if (NEWTON_BOND || i1 < nlocal) { + f[i1][0] += delx * fbond; + f[i1][1] += dely * fbond; + f[i1][2] += delz * fbond; + } + + if (NEWTON_BOND || i2 < nlocal) { + f[i2][0] -= delx * fbond; + f[i2][1] -= dely * fbond; + f[i2][2] -= delz * fbond; + } + + double ebond = 0.0; + if (EFLAG) { + double &r_pot = bondpot[idx].getVariableReference("r"); + r_pot = dr; + ebond = bondpot[idx].evaluate(); + } + if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondLepton::allocate() +{ + allocated = 1; + const int np1 = atom->nbondtypes + 1; + + memory->create(r0, np1, "bond:r0"); + memory->create(type2expression, np1, "bond:type2expression"); + memory->create(setflag, np1, "bond:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void BondLepton::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR, "Incorrect number of args for bond coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); + + double r0_one = utils::numeric(FLERR, arg[1], false, lmp); + + // remove whitespace and quotes from expression string and then + // check if the expression can be parsed and evaluated without error + std::string exp_one; + for (const auto &c : std::string(arg[2])) + if (!isspace(c) && (c != '"') && (c != '\'')) exp_one.push_back(c); + + try { + auto parsed = LMP_Lepton::Parser::parse(exp_one); + auto bondpot = parsed.createCompiledExpression(); + auto bondforce = parsed.differentiate("r").createCompiledExpression(); + double &r_pot = bondpot.getVariableReference("r"); + double &r_for = bondforce.getVariableReference("r"); + r_for = r_pot = 1.0; + bondpot.evaluate(); + bondforce.evaluate(); + } catch (std::exception &e) { + error->all(FLERR, e.what()); + } + + std::size_t idx = 0; + for (const auto &exp : expressions) { + if (exp == exp_one) break; + ++idx; + } + + // not found, add to list + if ((expressions.size() == 0) || (idx == expressions.size())) expressions.push_back(exp_one); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + r0[i] = r0_one; + type2expression[i] = idx; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + return an equilbrium bond length +------------------------------------------------------------------------- */ + +double BondLepton::equilibrium_distance(int i) +{ + return r0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void BondLepton::write_restart(FILE *fp) +{ + fwrite(&r0[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&type2expression[1], sizeof(int), atom->nbondtypes, fp); + + int num = expressions.size(); + int maxlen = 0; + for (const auto &exp : expressions) maxlen = MAX(maxlen, (int) exp.size()); + ++maxlen; + + fwrite(&num, sizeof(int), 1, fp); + fwrite(&maxlen, sizeof(int), 1, fp); + for (const auto &exp : expressions) { + int n = exp.size() + 1; + fwrite(&n, sizeof(int), 1, fp); + fwrite(exp.c_str(), sizeof(char), n, fp); + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void BondLepton::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &r0[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &type2expression[1], sizeof(int), atom->nbondtypes, fp, nullptr, error); + } + MPI_Bcast(&r0[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&type2expression[1], atom->nbondtypes, MPI_INT, 0, world); + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; + + int num, maxlen, len; + if (comm->me == 0) { + utils::sfread(FLERR, &num, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &maxlen, sizeof(int), 1, fp, nullptr, error); + } + MPI_Bcast(&num, 1, MPI_INT, 0, world); + MPI_Bcast(&maxlen, 1, MPI_INT, 0, world); + + char *buf = new char[maxlen]; + + for (int i = 0; i < num; ++i) { + if (comm->me == 0) { + utils::sfread(FLERR, &len, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, buf, sizeof(char), len, fp, nullptr, error); + } + MPI_Bcast(buf, maxlen, MPI_CHAR, 0, world); + expressions.push_back(buf); + } + + delete[] buf; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondLepton::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) + fprintf(fp, "%d %g %s\n", i, r0[i], expressions[type2expression[i]].c_str()); +} + +/* ---------------------------------------------------------------------- */ + +double BondLepton::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce) +{ + const double r = sqrt(rsq); + const double dr = r - r0[type]; + + auto parsed = LMP_Lepton::Parser::parse(expressions[type2expression[type]]); + auto bondpot = parsed.createCompiledExpression(); + auto bondforce = parsed.differentiate("r").createCompiledExpression(); + double &r_for = bondforce.getVariableReference("r"); + double &r_pot = bondpot.getVariableReference("r"); + r_for = r_pot = dr; + + // force and energy + + fforce = 0.0; + double ebond = 0.0; + if (r > 0.0) { + fforce = -bondforce.evaluate() / r; + ebond = bondpot.evaluate(); + } + return ebond; +} + +/* ---------------------------------------------------------------------- */ + +void *BondLepton::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "r0") == 0) return (void *) r0; + return nullptr; +} diff --git a/src/LEPTON/bond_lepton.h b/src/LEPTON/bond_lepton.h new file mode 100644 index 0000000000..5c430b7e63 --- /dev/null +++ b/src/LEPTON/bond_lepton.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef BOND_CLASS +// clang-format off +BondStyle(lepton,BondLepton); +// clang-format on +#else + +#ifndef LMP_BOND_LEPTON_H +#define LMP_BOND_LEPTON_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondLepton : public Bond { + public: + BondLepton(class LAMMPS *); + ~BondLepton() override; + void compute(int, int) override; + void coeff(int, char **) override; + double equilibrium_distance(int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + double single(int, double, int, int, double &) override; + void *extract(const char *, int &) override; + + protected: + std::vector expressions; + double *r0; + int *type2expression; + + virtual void allocate(); + + private: + template void eval(); +}; +} // namespace LAMMPS_NS +#endif +#endif diff --git a/src/OPENMP/bond_lepton_omp.cpp b/src/OPENMP/bond_lepton_omp.cpp new file mode 100644 index 0000000000..3171aaa51c --- /dev/null +++ b/src/OPENMP/bond_lepton_omp.cpp @@ -0,0 +1,148 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "bond_lepton_omp.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "omp_compat.h" + +#include + +#include "LMP_Lepton.h" +#include "suffix.h" +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondLeptonOMP::BondLeptonOMP(class LAMMPS *lmp) : BondLepton(lmp), ThrOMP(lmp, THR_BOND) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void BondLeptonOMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = neighbor->nbondlist; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (inum > 0) { + if (evflag) { + if (eflag) { + if (force->newton_bond) + eval<1, 1, 1>(ifrom, ito, thr); + else + eval<1, 1, 0>(ifrom, ito, thr); + } else { + if (force->newton_bond) + eval<1, 0, 1>(ifrom, ito, thr); + else + eval<1, 0, 0>(ifrom, ito, thr); + } + } else { + if (force->newton_bond) + eval<0, 0, 1>(ifrom, ito, thr); + else + eval<0, 0, 0>(ifrom, ito, thr); + } + } + thr->timer(Timer::BOND); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +/* ---------------------------------------------------------------------- */ + +template +void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr) +{ + std::vector bondforce; + std::vector bondpot; + for (const auto &expr : expressions) { + auto parsed = LMP_Lepton::Parser::parse(expr); + bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); + if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression()); + } + + const auto *_noalias const x = (dbl3_t *) atom->x[0]; + auto *_noalias const f = (dbl3_t *) thr->get_f()[0]; + const int3_t *_noalias const bondlist = (int3_t *) neighbor->bondlist[0]; + const int nlocal = atom->nlocal; + + for (int n = nfrom; n < nto; n++) { + const int i1 = bondlist[n].a; + const int i2 = bondlist[n].b; + const int type = bondlist[n].t; + + const double delx = x[i1].x - x[i2].x; + const double dely = x[i1].y - x[i2].y; + const double delz = x[i1].z - x[i2].z; + + const double rsq = delx * delx + dely * dely + delz * delz; + const double r = sqrt(rsq); + const double dr = r - r0[type]; + const int idx = type2expression[type]; + + // force and energy + + double fbond = 0.0; + if (r > 0.0) { + double &r_for = bondforce[idx].getVariableReference("r"); + r_for = dr; + fbond = -bondforce[idx].evaluate() / r; + } + + // apply force to each of 2 atoms + + if (NEWTON_BOND || i1 < nlocal) { + f[i1].x += delx * fbond; + f[i1].y += dely * fbond; + f[i1].z += delz * fbond; + } + + if (NEWTON_BOND || i2 < nlocal) { + f[i2].x -= delx * fbond; + f[i2].y -= dely * fbond; + f[i2].z -= delz * fbond; + } + + double ebond = 0.0; + if (EFLAG) { + double &r_pot = bondpot[idx].getVariableReference("r"); + r_pot = dr; + ebond = bondpot[idx].evaluate(); + } + if (EVFLAG) + ev_tally_thr(this, i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz, thr); + } +} diff --git a/src/OPENMP/bond_lepton_omp.h b/src/OPENMP/bond_lepton_omp.h new file mode 100644 index 0000000000..bdcc36434e --- /dev/null +++ b/src/OPENMP/bond_lepton_omp.h @@ -0,0 +1,44 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef BOND_CLASS +// clang-format off +BondStyle(lepton/omp,BondLeptonOMP); +// clang-format on +#else + +#ifndef LMP_BOND_LEPTON_OMP_H +#define LMP_BOND_LEPTON_OMP_H + +#include "bond_lepton.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class BondLeptonOMP : public BondLepton, public ThrOMP { + + public: + BondLeptonOMP(class LAMMPS *lmp); + void compute(int, int) override; + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; +} // namespace LAMMPS_NS +#endif +#endif diff --git a/unittest/force-styles/tests/bond-lepton.yaml b/unittest/force-styles/tests/bond-lepton.yaml new file mode 100644 index 0000000000..a220dd8b0c --- /dev/null +++ b/unittest/force-styles/tests/bond-lepton.yaml @@ -0,0 +1,90 @@ +--- +lammps_version: 3 Nov 2022 +tags: generated +date_generated: Thu Dec 22 02:03:59 2022 +epsilon: 2.5e-13 +skip_tests: +prerequisites: ! | + atom full + bond lepton +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +bond_style: lepton +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 +...