tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4433 from farrelljd-iop/angle-mwlc

Browse Source 
Angle style MWLC (meltable wormlike chain)
This commit is contained in:
Axel Kohlmeyer
2025-01-11 18:07:36 -05:00
committed by GitHub
parent 38c3477d2f e3f930b33f
commit a443fe4142
9 changed files with 574 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_bond.rst
View File

@ -90,6 +90,7 @@ OPT.
* :doc:`lepton (o) <angle_lepton>`
* :doc:`mesocnt <angle_mesocnt>`
* :doc:`mm3 <angle_mm3>`
* :doc:`mwlc <angle_mwlc>`
* :doc:`quartic (o) <angle_quartic>`
* :doc:`spica (ko) <angle_spica>`
* :doc:`table (o) <angle_table>`

94
doc/src/angle_mwlc.rst Normal file
View File

@ -0,0 +1,94 @@
.. index:: angle_style mwlc
angle_style mwlc command
==========================
Syntax
""""""
.. code-block:: LAMMPS
angle_style mwlc
Examples
""""""""
.. code-block:: LAMMPS
angle_style mwlc
angle_coeff * 25 1 10 1
Description
"""""""""""
.. versionadded:: TBD
The *mwlc* angle style models a meltable wormlike chain and can be used
to model non-linear bending elasticity of polymers, e.g. DNA. *mwlc*
uses a potential that is a canonical-ensemble superposition of a
non-melted and a melted state :ref:`(Farrell) <Farrell>`. The potential
is
.. math::
E = -k_{B}T\,\log [q + q^{m}] + E_{0},
where the non-melted and melted partition functions are
.. math::
q = \exp [-k_{1}(1+\cos{\theta})/k_{B}T]; \\
q^{m} = \exp [-(\mu+k_{2}(1+\cos{\theta}))/k_{B}T].
:math:`k_1` is the bending elastic constant of the non-melted state,
:math:`k_2` is the bending elastic constant of the melted state,
:math:`\mu` is the melting energy, and
:math:`T` is the reference temperature.
The reference energy,
.. math::
E_{0} = -k_{B}T\,\log [1 + \exp[-\mu/k_{B}T]],
ensures that E is zero for a fully extended chain.
This potential is a continuous version of the two-state potential
introduced by :ref:`(Yan) <Yan>`.
The following coefficients must be defined for each angle type via the
:doc:`angle_coeff <angle_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
or :doc:`read_restart <read_restart>` commands:
* :math:`k_1` (energy)
* :math:`k_2` (energy)
* :math:`\mu` (energy)
* :math:`T` (temperature)
----------
Restrictions
""""""""""""
This angle style can only be used if LAMMPS was built with the
EXTRA-MOLECULE package. See the :doc:`Build package <Build_package>`
doc page for more info.
Related commands
""""""""""""""""
:doc:`angle_coeff <angle_coeff>`
Default
"""""""
none
----------
.. _Farrell:
**(Farrell)** `Farrell, Dobnikar, Podgornik, Curk, Phys Rev Lett, 133, 148101 (2024). <https://doi.org/10.1103/PhysRevLett.133.148101>`_
.. _Yan:
**(Yan)** `Yan, Marko, Phys Rev Lett, 93, 108108 (2004). <https://doi.org/10.1103/PhysRevLett.93.108108>`_

1
doc/src/angle_style.rst
View File

@ -94,6 +94,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist.
* :doc:`lepton <angle_lepton>` - angle potential from evaluating a string
* :doc:`mesocnt <angle_mesocnt>` - piecewise harmonic and linear angle for bending-buckling of nanotubes
* :doc:`mm3 <angle_mm3>` - anharmonic angle
* :doc:`mwlc <angle_mwlc>` - meltable wormlike chain
* :doc:`quartic <angle_quartic>` - angle with cubic and quartic terms
* :doc:`spica <angle_spica>` - harmonic angle with repulsive SPICA pair style between 1-3 atoms
* :doc:`table <angle_table>` - tabulated by angle

2
doc/src/fix_adapt.rst
View File

@ -406,6 +406,8 @@ sub-style name. The angle styles that currently work with fix adapt are:
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`mm3 <angle_mm3>` | k,theta0 | type angles |
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`mwlc <angle_mwlc>` | k1,k2,mu,T | type angles |
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`quartic <angle_quartic>` | k2,k3,k4,theta0 | type angles |
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`spica <angle_spica>` | k,theta0 | type angles |

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

@ -2116,6 +2116,7 @@ Marchi
Mariella
Marinica
Markland
Marko
Marrink
Marroquin
Marsaglia
@ -2198,6 +2199,7 @@ Meissner
Melchor
Meloni
Melrose
meltable
mem
Mem
memalign
@ -2408,6 +2410,7 @@ mV
Mvapich
mvh
mvv
mwlc
MxN
myCompute
myIndex
@ -2929,6 +2932,7 @@ Pmoltrans
pN
png
podd
Podgornik
Podhorszki
Poiseuille
poisson
@ -4113,6 +4117,7 @@ workflow
workflows
Workum
Worley
wormlike
Wriggers
writedata
Wuppertal
@ -4181,6 +4186,7 @@ yaff
YAFF
Yamada
yaml
Yan
Yanxon
Yaser
Yazdani

2
src/.gitignore vendored
View File

@ -507,6 +507,8 @@
/angle_harmonic.h
/angle_mm3.cpp
/angle_mm3.h
/angle_mwlc.cpp
/angle_mwlc.h
/angle_quartic.cpp
/angle_quartic.h
/angle_spica.cpp

326
src/EXTRA-MOLECULE/angle_mwlc.cpp Normal file
View File

@ -0,0 +1,326 @@
/* ----------------------------------------------------------------------
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: James D. Farrell (IoP CAS) j.d.farrell at gmail.com
[ based on angle_cosine.cpp ]
------------------------------------------------------------------------- */
#include "angle_mwlc.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
/* ---------------------------------------------------------------------- */
AngleMWLC::AngleMWLC(LAMMPS *_lmp) : Angle(_lmp)
{
born_matrix_enable = 1;
}
/* ---------------------------------------------------------------------- */
AngleMWLC::~AngleMWLC()
{
if (allocated && !copymode) {
memory->destroy(setflag);
memory->destroy(k1);
memory->destroy(k2);
memory->destroy(mu);
memory->destroy(temp);
}
}
/* ---------------------------------------------------------------------- */
void AngleMWLC::compute(int eflag, int vflag)
{
int i1, i2, i3, n, type;
double delx1, dely1, delz1, delx2, dely2, delz2;
double eangle, f1[3], f3[3];
double rsq1, rsq2, r1, r2, c, a, a11, a12, a22;
double q, qm, Q;
eangle = 0.0;
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
int **anglelist = neighbor->anglelist;
int nanglelist = neighbor->nanglelist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
double kbt, v_min;
for (n = 0; n < nanglelist; n++) {
i1 = anglelist[n][0];
i2 = anglelist[n][1];
i3 = anglelist[n][2];
type = anglelist[n][3];
kbt = temp[type] * force->boltz;
v_min = -kbt * log(1 + exp(-mu[type] / kbt));
// 1st bond
delx1 = x[i1][0] - x[i2][0];
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
r2 = sqrt(rsq2);
// c = cosine of angle
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
// force & energy
q = exp(-k1[type] * (1.0 + c) / kbt);
qm = exp((-k2[type] * (1.0 + c) - mu[type]) / kbt);
Q = q + qm;
if (eflag) eangle = -kbt * log(Q) - v_min;
a = (k1[type] * q + k2[type] * qm) / Q;
a11 = a * c / rsq1;
a12 = -a / (r1 * r2);
a22 = a * c / rsq2;
f1[0] = a11 * delx1 + a12 * delx2;
f1[1] = a11 * dely1 + a12 * dely2;
f1[2] = a11 * delz1 + a12 * delz2;
f3[0] = a22 * delx2 + a12 * delx1;
f3[1] = a22 * dely2 + a12 * dely1;
f3[2] = a22 * delz2 + a12 * delz1;
// apply force to each of 3 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= f1[0] + f3[0];
f[i2][1] -= f1[1] + f3[1];
f[i2][2] -= f1[2] + f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
}
if (evflag)
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
delz2);
}
}
/* ---------------------------------------------------------------------- */
void AngleMWLC::allocate()
{
allocated = 1;
const int np1 = atom->nangletypes + 1;
memory->create(k1, np1, "angle:k1");
memory->create(k2, np1, "angle:k2");
memory->create(mu, np1, "angle:mu");
memory->create(temp, np1, "angle:temp");
memory->create(setflag, np1, "angle:setflag");
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one type
------------------------------------------------------------------------- */
void AngleMWLC::coeff(int narg, char **arg)
{
if (narg != 5) error->all(FLERR, "Incorrect args for angle coefficients");
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
double k1_one = utils::numeric(FLERR, arg[1], false, lmp);
double k2_one = utils::numeric(FLERR, arg[2], false, lmp);
double mu_one = utils::numeric(FLERR, arg[3], false, lmp);
double temp_one = utils::numeric(FLERR, arg[4], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k1[i] = k1_one;
k2[i] = k2_one;
mu[i] = mu_one;
temp[i] = temp_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
}
/* ---------------------------------------------------------------------- */
double AngleMWLC::equilibrium_angle(int /*i*/)
{
return MY_PI;
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void AngleMWLC::write_restart(FILE *fp)
{
fwrite(&k1[1], sizeof(double), atom->nangletypes, fp);
fwrite(&k2[1], sizeof(double), atom->nangletypes, fp);
fwrite(&mu[1], sizeof(double), atom->nangletypes, fp);
fwrite(&temp[1], sizeof(double), atom->nangletypes, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void AngleMWLC::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
utils::sfread(FLERR, &k1[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &k2[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &mu[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &temp[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
}
MPI_Bcast(&k1[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&k2[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&mu[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&temp[1], atom->nangletypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void AngleMWLC::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp, "%d %g %g %g %g\n", i, k1[i], k2[i], mu[i], temp[i]);
}
/* ---------------------------------------------------------------------- */
double AngleMWLC::single(int type, int i1, int i2, int i3)
{
double **x = atom->x;
double kbt = temp[type] * force->boltz;
double v_min = -kbt * log(1 + exp(-mu[type] / kbt));
double delx1 = x[i1][0] - x[i2][0];
double dely1 = x[i1][1] - x[i2][1];
double delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1, dely1, delz1);
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
double delx2 = x[i3][0] - x[i2][0];
double dely2 = x[i3][1] - x[i2][1];
double delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2, dely2, delz2);
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
double q = exp(-k1[type] * (1.0 + c) / kbt);
double qm = exp((-k2[type] * (1.0 + c) - mu[type]) / kbt);
return -kbt * log(q + qm) - v_min;
}
/* ---------------------------------------------------------------------- */
void AngleMWLC::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2)
{
double **x = atom->x;
double kbt = temp[type] * force->boltz;
double delx1 = x[i1][0] - x[i2][0];
double dely1 = x[i1][1] - x[i2][1];
double delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1, dely1, delz1);
double delx2 = x[i3][0] - x[i2][0];
double dely2 = x[i3][1] - x[i2][1];
double delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2, dely2, delz2);
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= sqrt((delx1 * delx1 + dely1 * dely1 + delz1 * delz1) *
(delx2 * delx2 + dely2 * dely2 + delz2 * delz2));
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
const double q = exp(-k1[type] * (1.0 + c) / kbt);
const double qm = exp((-k2[type] * (1.0 + c) - mu[type]) / kbt);
const double Q = q + qm;
du = (k1[type] * q + k2[type] * qm) / Q;
du2 = (k1[type] - k2[type]) / Q;
du2 *= -du2 * q * qm / kbt;
}
/* ----------------------------------------------------------------------
return ptr to internal members upon request
------------------------------------------------------------------------ */
void *AngleMWLC::extract(const char *str, int &dim)
{
dim = 1;
if (strcmp(str, "k1") == 0) return (void *) k1;
if (strcmp(str, "k2") == 0) return (void *) k2;
if (strcmp(str, "mu") == 0) return (void *) mu;
if (strcmp(str, "temp") == 0) return (void *) temp;
return nullptr;
}

53
src/EXTRA-MOLECULE/angle_mwlc.h Normal file
View File

@ -0,0 +1,53 @@
/* -*- 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 ANGLE_CLASS
// clang-format off
AngleStyle(mwlc,AngleMWLC);
// clang-format on
#else
#ifndef LMP_ANGLE_MWLC_H
#define LMP_ANGLE_MWLC_H
#include "angle.h"
namespace LAMMPS_NS {
class AngleMWLC : public Angle {
public:
AngleMWLC(class LAMMPS *);
~AngleMWLC() override;
void compute(int, int) override;
void coeff(int, char **) override;
double equilibrium_angle(int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
void *extract(const char *, int &) override;
protected:
double *k1;
double *k2;
double *mu;
double *temp;
virtual void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif

89
unittest/force-styles/tests/angle-mwlc.yaml Normal file
View File

@ -0,0 +1,89 @@
---
lammps_version: 19 Nov 2024
tags: generated
date_generated: Fri Jan 10 13:54:41 2025
epsilon: 2.5e-13
skip_tests:
prerequisites: ! |
atom full
angle mwlc
pre_commands: ! ""
post_commands: ! ""
input_file: in.fourmol
angle_style: mwlc
angle_coeff: ! |
* 6.0 0.6 15.0 300.0
equilibrium: 4 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589793
extract: ! |
k1 1
k2 1
mu 1
temp 1
natoms: 29
init_energy: 115.97311255470154
init_stress: ! |2-
1.0785144954758131e+01 -5.1561329771224154e+00 -5.6290119776357184e+00 1.7509753155291005e+01 8.4620050528315804e+00 5.4460179571966592e+00
init_forces: ! |2
1 2.0049292277552344e+00 -3.4967838538707023e+00 -7.7679149207041780e+00
2 5.2379550897725657e-01 4.2410446627574743e+00 3.4103186723411198e+00
3 -2.6364713466403744e+00 1.7646356480410195e+00 6.7422761915790446e+00
4 5.9780828166032607e-01 -1.0977711474862293e+00 -1.2762851805838284e+00
5 2.7519993963547762e-02 -1.3655215729181727e+00 -4.0557483263144212e-01
6 -1.4407532935169898e-01 -1.5931923885520087e+00 -1.8663167177332292e+00
7 -4.1794796128320089e-01 4.1379756743889873e-01 -9.9274508603132139e-02
8 8.6885282425298183e-01 2.9273693794680500e-01 3.9875384936490477e+00
9 -7.4989162989673552e-01 1.9120936074864758e-02 1.6480002237326197e-01
10 4.5851731938743301e+00 5.1383619725146579e-01 -7.8986229128522627e+00
11 -1.3361238263371962e+00 -1.6154590769980759e+00 1.3132946273952864e+00
12 -9.8353444089007258e-01 5.7660344986117273e-01 -1.6760444984003993e+00
13 4.8213326840430115e-02 7.5476102608997231e-02 9.7662630040726750e-01
14 1.7020627463903515e-01 -4.9418713351985666e-02 -6.9510195676192543e-02
15 5.7466130382222369e-01 -3.0791486600816609e-01 9.0601148952989252e-01
16 -5.6604535775268596e+00 5.0708111379184277e+00 1.7013150451815964e+00
17 2.5273381761407760e+00 -3.4420010207137852e+00 1.8573629247281453e+00
18 2.0828420935010161e-01 2.1765636018344789e+00 -8.8331994782253780e+00
19 -2.5502877667816710e+00 -3.1771408888391202e+00 3.7828819753880119e+00
20 2.3420035574315694e+00 1.0005772870046410e+00 5.0503175028373661e+00
21 2.2520112617177057e+00 2.8191460010961373e+00 -8.4251146760580085e+00
22 -4.0101778142198494e+00 -2.7341370977288624e+00 2.9783951433938527e+00
23 1.7581665525021437e+00 -8.5008903367275068e-02 5.4467195326641562e+00
24 -1.9792115156439898e+00 7.6933124498060685e+00 -4.5803031378051688e+00
25 -1.5876110942732853e+00 -5.3190980720841390e+00 8.6990568523236744e-01
26 3.5668226099172751e+00 -2.3742143777219300e+00 3.7103974525728014e+00
27 -1.0535801895693373e+00 8.6490680045193269e+00 -2.9461596365031912e+00
28 -2.3037157984215133e+00 -5.1714048596555164e+00 1.2013129729191907e-01
29 3.3572959879908506e+00 -3.4776631448638109e+00 2.8260283392112724e+00
run_energy: 115.88929660500146
run_stress: ! |2-
1.0741962021125900e+01 -5.1494210571452070e+00 -5.5925409639806976e+00 1.7499135911989786e+01 8.4401379697723158e+00 5.4365488644723161e+00
run_forces: ! |2
1 2.0024084967323734e+00 -3.4852712888973301e+00 -7.7681384109278948e+00
2 5.1802520314744682e-01 4.2290273746409701e+00 3.4156643023763298e+00
3 -2.6176139409651311e+00 1.7481727755506558e+00 6.7276571726932186e+00
4 5.9646423832275675e-01 -1.0896302596160039e+00 -1.2658946496090273e+00
5 2.0626395368878558e-02 -1.3553597325783064e+00 -4.0428709993242351e-01
6 -1.5043849290340794e-01 -1.6020934387436645e+00 -1.8709909377854954e+00
7 -4.1644858135120844e-01 4.1712949378864117e-01 -1.0014611176127630e-01
8 8.7761514222742898e-01 2.9932822072179732e-01 3.9914298835888857e+00
9 -7.5185528312840066e-01 1.9244616790871305e-02 1.6542886887136499e-01
10 4.5736254271229742e+00 4.9828467262363496e-01 -7.8986022120685373e+00
11 -1.3323136450232624e+00 -1.6103729207690813e+00 1.3097303435864336e+00
12 -1.0040757700004965e+00 5.7836875530999343e-01 -1.6448472694236678e+00
13 4.8392108795294986e-02 7.7447222505411639e-02 9.6494421976551426e-01
14 1.8086042252760026e-01 -5.1422366669672037e-02 -7.3717779732260791e-02
15 5.8324025586810402e-01 -3.0468600966273196e-01 8.9577956208287191e-01
16 -5.6548207755179911e+00 5.0752603376691656e+00 1.6979261789431477e+00
17 2.5263087987770385e+00 -3.4434274526643489e+00 1.8580639393328193e+00
18 2.1733330742608858e-01 2.1830182266117371e+00 -8.8242777257965024e+00
19 -2.5415995943828809e+00 -3.1730376716173803e+00 3.7883421709467422e+00
20 2.3242662869567923e+00 9.9001944500564321e-01 5.0359355548497611e+00
21 2.2580025853773256e+00 2.8068061574356138e+00 -8.4373139893370634e+00
22 -4.0070465500208554e+00 -2.7223210888728504e+00 2.9888102357206301e+00
23 1.7490439646435296e+00 -8.4485068562763388e-02 5.4485037536164329e+00
24 -1.9860712475490458e+00 7.7098710540424751e+00 -4.5855114458364374e+00
25 -1.5787959452376497e+00 -5.3290411597315357e+00 8.7215941857315382e-01
26 3.5648671927866955e+00 -2.3808298943109389e+00 3.7133520272632836e+00
27 -1.0631738916816893e+00 8.6588340987013801e+00 -2.9324082303496679e+00
28 -2.2973201808293515e+00 -5.1771902142248960e+00 1.1264225202087563e-01
29 3.3604940725110408e+00 -3.4816438844764841e+00 2.8197659783287921e+00
...