ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4191 from akohlmey/pair-hybrid-molecule

Browse Source 
Add new pair style hybrid/molecular
This commit is contained in:
Axel Kohlmeyer
2024-06-19 00:49:37 -04:00
committed by GitHub
parent b9f8a5b811 71b38521c5
commit 37572225f4
22 changed files with 577 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
doc/src/Commands_pair.rst
View File

@ -25,16 +25,16 @@ OPT.
* :doc:`none <pair_none>`
* :doc:`zero <pair_zero>`
* :doc:`hybrid (k) <pair_hybrid>`
* :doc:`hybrid/overlay (k) <pair_hybrid>`
* :doc:`hybrid/scaled <pair_hybrid>`
* :doc:`hybrid (ko) <pair_hybrid>`
* :doc:`hybrid/molecular (o) <pair_hybrid>`
* :doc:`hybrid/overlay (ko) <pair_hybrid>`
* :doc:`hybrid/scaled (o) <pair_hybrid>`
* :doc:`kim <pair_kim>`
* :doc:`list <pair_list>`
* :doc:`tracker <pair_tracker>`
*
*
*
*
* :doc:`adp (ko) <pair_adp>`
* :doc:`agni (o) <pair_agni>`
* :doc:`aip/water/2dm (t) <pair_aip_water_2dm>`

54
doc/src/pair_hybrid.rst
View File

@ -1,28 +1,41 @@
.. index:: pair_style hybrid
.. index:: pair_style hybrid/kk
.. index:: pair_style hybrid/omp
.. index:: pair_style hybrid/molecular
.. index:: pair_style hybrid/molecular/omp
.. index:: pair_style hybrid/overlay
.. index:: pair_style hybrid/overlay/omp
.. index:: pair_style hybrid/overlay/kk
.. index:: pair_style hybrid/scaled
.. index:: pair_style hybrid/scaled/omp
pair_style hybrid command
=========================
Accelerator Variants: *hybrid/kk*
Accelerator Variants: *hybrid/kk*, *hybrid/omp*
pair_style hybrid/molecular command
===================================
Accelerator Variant: *hybrid/molecular/omp*
pair_style hybrid/overlay command
=================================
Accelerator Variants: *hybrid/overlay/kk*
Accelerator Variants: *hybrid/overlay/kk*, *hybrid/overlay/omp*
pair_style hybrid/scaled command
==================================
Accelerator Variant: *hybrid/scaled/omp*
Syntax
""""""
.. code-block:: LAMMPS
pair_style hybrid style1 args style2 args ...
pair_style hybrid/molecular factor1 style1 args factor2 style 2 args
pair_style hybrid/overlay style1 args style2 args ...
pair_style hybrid/scaled factor1 style1 args factor2 style 2 args ...
@ -47,6 +60,10 @@ Examples
pair_coeff * * tersoff Si.tersoff Si
pair_coeff * * sw Si.sw Si
pair_style hybrid/molecular lj/cut 2.5 lj/cut 2.5
pair_coeff * * lj/cut 1 1.0 1.0
pair_coeff * * lj/cut 2 1.5 1.0
variable one equal ramp(1.0,0.0)
variable two equal 1.0-v_one
pair_style hybrid/scaled v_one lj/cut 2.5 v_two morse 2.5
@ -56,17 +73,26 @@ Examples
Description
"""""""""""
The *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles enable the
use of multiple pair styles in one simulation. With the *hybrid* style,
exactly one pair style is assigned to each pair of atom types. With the
*hybrid/overlay* and *hybrid/scaled* styles, one or more pair styles can
be assigned to each pair of atom types. The assignment of pair styles
to type pairs is made via the :doc:`pair_coeff <pair_coeff>` command.
The major difference between the *hybrid/overlay* and *hybrid/scaled*
styles is that the *hybrid/scaled* adds a scale factor for each
sub-style contribution to forces, energies and stresses. Because of the
added complexity, the *hybrid/scaled* style has more overhead and thus
may be slower than *hybrid/overlay*.
The *hybrid*, *hybrid/overlay*, *hybrid/molecular*, and *hybrid/scaled*
styles enable the use of multiple pair styles in one simulation. With
the *hybrid* style, exactly one pair style is assigned to each pair of
atom types. With the *hybrid/overlay* and *hybrid/scaled* styles, one
or more pair styles can be assigned to each pair of atom types. With
the hybrid/molecular style, pair styles are assigned to either intra-
or inter-molecular interactions.
The assignment of pair styles to type pairs is made via the
:doc:`pair_coeff <pair_coeff>` command. The major difference between
the *hybrid/overlay* and *hybrid/scaled* styles is that the
*hybrid/scaled* adds a scale factor for each sub-style contribution to
forces, energies and stresses. Because of the added complexity, the
*hybrid/scaled* style has more overhead and thus may be slower than
*hybrid/overlay*.
The *hybrid/molecular* pair style accepts *only* two sub-styles: the
first is assigned to intra-molecular interactions (i.e. both atoms
have the same molecule ID), the second to inter-molecular interactions
(i.e. interacting atoms have different molecule IDs).
Here are two examples of hybrid simulations. The *hybrid* style could
be used for a simulation of a metal droplet on a LJ surface. The metal
@ -476,6 +502,8 @@ the same or else LAMMPS will generate an error.
Pair style *hybrid/scaled* currently only works for non-accelerated
pair styles and pair styles from the OPT package.
Pair style *hybrid/molecular* is not compatible with manybody potentials.
When using pair styles from the GPU package they must not be listed
multiple times. LAMMPS will detect this and abort.

1
doc/src/pair_style.rst
View File

@ -108,6 +108,7 @@ accelerated styles exist.
* :doc:`none <pair_none>` - turn off pairwise interactions
* :doc:`hybrid <pair_hybrid>` - multiple styles of pairwise interactions
* :doc:`hybrid/molecular <pair_hybrid>` - different pair styles for intra- and inter-molecular interactions
* :doc:`hybrid/overlay <pair_hybrid>` - multiple styles of superposed pairwise interactions
* :doc:`hybrid/scaled <pair_hybrid>` - multiple styles of scaled superposed pairwise interactions
* :doc:`zero <pair_zero>` - neighbor list but no interactions

22
src/INTEL/npair_skip_intel.cpp
View File

@ -74,9 +74,13 @@ void NPairSkipIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
const int nlocal = atom->nlocal;
const int e_nall = nlocal + atom->nghost;
const int * _noalias const type = atom->type;
const tagint * _noalias const molecule = atom->molecule;
int * _noalias const ilist = list->ilist;
int * _noalias const numneigh = list->numneigh;
int ** _noalias const firstneigh = (int ** const)list->firstneigh; // NOLINT
const int molskip = list->molskip;
const int * _noalias const ilist_skip = list->listskip->ilist;
const int * _noalias const numneigh_skip = list->listskip->numneigh;
const int ** _noalias const firstneigh_skip = (const int ** const)list->listskip->firstneigh; // NOLINT
@ -110,7 +114,7 @@ void NPairSkipIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
for (int ii = ifrom; ii < ito; ii++) {
const int i = ilist_skip[ii];
const int itype = type[i];
if (iskip[itype]) continue;
if (!molskip && iskip[itype]) continue;
int n = 0;
int *neighptr = ipage.vget();
@ -142,7 +146,11 @@ void NPairSkipIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
for (int jj = 0; jj < jnum; jj++) {
const int joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
if (!ijskip[itype][type[j]]) neighptr[n++] = joriginal;
if (!molskip && ijskip[itype][type[j]]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
neighptr[n++] = joriginal;
}
}
@ -269,9 +277,13 @@ void NPairSkipTrimIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
const int e_nall = nlocal + atom->nghost;
const ATOM_T * _noalias const x = buffers->get_x();
const int * _noalias const type = atom->type;
const tagint * _noalias const molecule = atom->molecule;
int * _noalias const ilist = list->ilist;
int * _noalias const numneigh = list->numneigh;
int ** _noalias const firstneigh = (int ** const)list->firstneigh; // NOLINT
const int molskip = list->molskip;
const int * _noalias const ilist_skip = list->listskip->ilist;
const int * _noalias const numneigh_skip = list->listskip->numneigh;
const int ** _noalias const firstneigh_skip = (const int ** const)list->listskip->firstneigh; // NOLINT
@ -306,7 +318,7 @@ void NPairSkipTrimIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
for (int ii = ifrom; ii < ito; ii++) {
const int i = ilist_skip[ii];
const int itype = type[i];
if (iskip[itype]) continue;
if (!molskip && iskip[itype]) continue;
const flt_t xtmp = x[i].x;
const flt_t ytmp = x[i].y;
@ -370,7 +382,9 @@ void NPairSkipTrimIntel::build_t(NeighList *list, int *numhalf, int *cnumneigh,
const int j = joriginal & NEIGHMASK;
int addme = 1;
if (ijskip[itype][type[j]]) addme = 0;
if (!molskip && ijskip[itype][type[j]]) addme = 0;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) addme = 0;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) addme = 0;
// trim to shorter cutoff

8
src/OPENMP/fix_omp.cpp
View File

@ -229,7 +229,13 @@ void FixOMP::init()
check_hybrid = 0; \
if (force->name) { \
if ( (strcmp(force->name ## _style,"hybrid") == 0) || \
(strcmp(force->name ## _style,"hybrid/overlay") == 0) ) \
(strcmp(force->name ## _style,"hybrid/overlay") == 0) || \
(strcmp(force->name ## _style,"hybrid/scaled") == 0) || \
(strcmp(force->name ## _style,"hybrid/molecular") == 0) || \
(strcmp(force->name ## _style,"hybrid/omp") == 0) || \
(strcmp(force->name ## _style,"hybrid/overlay/omp") == 0) || \
(strcmp(force->name ## _style,"hybrid/scaled/omp") == 0) || \
(strcmp(force->name ## _style,"hybrid/molecular/omp") == 0) ) \
check_hybrid=1; \
if (force->name->suffix_flag & Suffix::OMP) { \
last_force_name = (const char *) #name; \

34
src/OPENMP/npair_skip_omp.h
View File

@ -19,14 +19,12 @@
NPairStyle(skip/omp,
NPairSkip,
NP_SKIP | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_OMP);
NPairStyle(skip/half/respa/omp,
NPairSkipRespa,
NP_SKIP | NP_RESPA | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_OMP);
NPairStyle(skip/half/size/omp,
@ -36,32 +34,27 @@ NPairStyle(skip/half/size/omp,
NPairStyle(skip/size/off2on/omp,
NPairSkipSizeOff2on,
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_MULTI_OLD |
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_OMP);
NPairStyle(skip/size/off2on/oneside/omp,
NPairSkipSizeOff2onOneside,
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF |
NP_ORTHO | NP_TRI | NP_OMP);
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_OMP);
NPairStyle(skip/ghost/omp,
NPairSkip,
NP_SKIP | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_OMP | NP_GHOST);
NPairStyle(skip/trim/omp,
NPairSkipTrim,
NP_SKIP | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP);
NPairStyle(skip/trim/half/respa/omp,
NPairSkipTrimRespa,
NP_SKIP | NP_RESPA | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP);
NPairStyle(skip/trim/half/size/omp,
@ -71,20 +64,17 @@ NPairStyle(skip/trim/half/size/omp,
NPairStyle(skip/trim/size/off2on/omp,
NPairSkipTrimSizeOff2on,
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_MULTI_OLD |
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP);
NPairStyle(skip/trim/size/off2on/oneside/omp,
NPairSkipTrimSizeOff2onOneside,
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF |
NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP);
NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI |
NP_MULTI_OLD | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP);
NPairStyle(skip/trim/ghost/omp,
NPairSkipTrim,
NP_SKIP | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM | NP_OMP | NP_GHOST);
// clang-format off
#endif

3
src/neigh_list.cpp
View File

@ -148,6 +148,7 @@ void NeighList::post_constructor(NeighRequest *nq)
copy = nq->copy;
trim = nq->trim;
id = nq->id;
molskip = nq->molskip;
if (nq->copy) {
listcopy = neighbor->lists[nq->copylist];
@ -157,6 +158,7 @@ void NeighList::post_constructor(NeighRequest *nq)
if (nq->skip) {
listskip = neighbor->lists[nq->skiplist];
if (!molskip) {
int ntypes = atom->ntypes;
iskip = new int[ntypes+1];
memory->create(ijskip,ntypes+1,ntypes+1,"neigh_list:ijskip");
@ -166,6 +168,7 @@ void NeighList::post_constructor(NeighRequest *nq)
for (j = 1; j <= ntypes; j++)
ijskip[i][j] = nq->ijskip[i][j];
}
}
if (nq->halffull)
listfull = neighbor->lists[nq->halffulllist];

1
src/neigh_list.h
View File

@ -46,6 +46,7 @@ class NeighList : protected Pointers {
int kk2cpu; // 1 if this list is copied from Kokkos to CPU
int copymode; // 1 if this is a Kokkos on-device copy
int id; // copied from neighbor list request
int molskip; // 1/2 if this is an intra-/inter-molecular skip list
// data structs to store neighbor pairs I,J and associated values

9
src/neigh_request.cpp
View File

@ -76,6 +76,7 @@ NeighRequest::NeighRequest(LAMMPS *_lmp) : Pointers(_lmp)
skip = 0;
iskip = nullptr;
ijskip = nullptr;
molskip = REGULAR;
// only set when command = 1;
@ -183,6 +184,8 @@ int NeighRequest::identical(NeighRequest *other)
int NeighRequest::same_skip(NeighRequest *other)
{
if (molskip != other->molskip) return 0;
const int ntypes = atom->ntypes;
int same = 1;
@ -307,6 +310,12 @@ void NeighRequest::set_skip(int *_iskip, int **_ijskip)
ijskip = _ijskip;
}
void NeighRequest::set_molskip(int _molskip)
{
skip = 1;
molskip = _molskip;
}
void NeighRequest::enable_full()
{
half = 0;

5
src/neigh_request.h
View File

@ -29,6 +29,9 @@ class NeighRequest : protected Pointers {
friend class NPairSkipTrimIntel;
friend class FixIntel;
public:
enum { REGULAR, INTRA, INTER };
protected:
void *requestor; // class that made request
int requestor_instance; // instance of that class (only Fix, Compute, Pair)
@ -88,6 +91,7 @@ class NeighRequest : protected Pointers {
int skip; // 1 if this list skips atom types from another list
int *iskip; // iskip[i] if atoms of type I are not in list
int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list
int molskip; // 0 reqular list, 1 keep only intra-molecular entries, 2 keep inter-molecular
// command_style only set if command = 1
// allows print_pair_info() to access command name
@ -137,6 +141,7 @@ class NeighRequest : protected Pointers {
void set_kokkos_device(int);
void set_kokkos_host(int);
void set_skip(int *, int **);
void set_molskip(int);
void enable_full();
void enable_ghost();
void enable_intel();

7
src/neighbor.cpp
View File

@ -1800,11 +1800,18 @@ void Neighbor::print_pairwise_info()
else
out += fmt::format(", half/full from ({})",rq->halffulllist+1);
} else if (rq->skip) {
if (rq->molskip) {
if (rq->trim)
out += fmt::format(", molskip trim from ({})",rq->skiplist+1);
else
out += fmt::format(", molskip from ({})",rq->skiplist+1);
} else {
if (rq->trim)
out += fmt::format(", skip trim from ({})",rq->skiplist+1);
else
out += fmt::format(", skip from ({})",rq->skiplist+1);
}
}
out += "\n";
// list of neigh list attributes

10
src/npair_skip.cpp
View File

@ -17,6 +17,7 @@
#include "error.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
using namespace LAMMPS_NS;
@ -41,11 +42,13 @@ void NPairSkipTemp<TRIM>::build(NeighList *list)
int *type = atom->type;
int nlocal = atom->nlocal;
tagint *molecule = atom->molecule;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
int molskip = list->molskip;
int *ilist_skip = list->listskip->ilist;
int *numneigh_skip = list->listskip->numneigh;
@ -71,7 +74,8 @@ void NPairSkipTemp<TRIM>::build(NeighList *list)
for (ii = 0; ii < num_skip; ii++) {
i = ilist_skip[ii];
itype = type[i];
if (iskip[itype]) continue;
if (!molskip && iskip[itype]) continue;
if (TRIM) {
xtmp = x[i][0];
@ -90,7 +94,9 @@ void NPairSkipTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (ijskip[itype][type[j]]) continue;
if (!molskip && ijskip[itype][type[j]]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
if (TRIM) {
delx = xtmp - x[j][0];

6
src/npair_skip.h
View File

@ -30,7 +30,8 @@ NPairStyle(skip/ghost,
typedef NPairSkipTemp<0> NPairSkipSize;
NPairStyle(skip/half/size,
NPairSkipSize,
NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_SIZE | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI);
typedef NPairSkipTemp<1> NPairSkipTrim;
@ -50,7 +51,8 @@ NPairStyle(skip/ghost/trim,
typedef NPairSkipTemp<1> NPairSkipTrimSize;
NPairStyle(skip/trim/half/size,
NPairSkipTrimSize,
NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_SKIP | NP_SIZE | NP_HALF | NP_FULL |
NP_NSQ | NP_BIN | NP_MULTI | NP_MULTI_OLD |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_TRIM);
// clang-format on

17
src/npair_skip_respa.cpp
View File

@ -18,6 +18,7 @@
#include "error.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
using namespace LAMMPS_NS;
@ -39,11 +40,13 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
int *neighptr, *jlist, *neighptr_inner, *neighptr_middle;
int *type = atom->type;
tagint *molecule = atom->molecule;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
int molskip = list->molskip;
int *ilist_skip = list->listskip->ilist;
int *numneigh_skip = list->listskip->numneigh;
@ -90,7 +93,7 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (ii = 0; ii < inum_skip; ii++) {
i = ilist_skip[ii];
itype = type[i];
if (iskip[itype]) continue;
if (!molskip && iskip[itype]) continue;
if (TRIM) {
xtmp = x[i][0];
@ -114,7 +117,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (ijskip[itype][type[j]]) continue;
if (!molskip && ijskip[itype][type[j]]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
if (TRIM) {
delx = xtmp - x[j][0];
@ -135,7 +140,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (ijskip[itype][type[j]]) continue;
if (!molskip && ijskip[itype][type[j]]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
if (TRIM) {
delx = xtmp - x[j][0];
@ -157,7 +164,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (ijskip[itype][type[j]]) continue;
if (!molskip && ijskip[itype][type[j]]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
if (TRIM) {
delx = xtmp - x[j][0];

11
src/pair_hybrid.cpp
View File

@ -321,11 +321,12 @@ void PairHybrid::settings(int narg, char **arg)
iarg = 0;
nstyles = 0;
const std::string mystyle = force->pair_style;
while (iarg < narg) {
if (utils::strmatch(arg[iarg],"^hybrid"))
error->all(FLERR,"Pair style hybrid cannot have hybrid as a sub-style");
error->all(FLERR,"Pair style {} cannot have hybrid as a sub-style", mystyle);
if (strcmp(arg[iarg],"none") == 0)
error->all(FLERR,"Pair style hybrid cannot have none as a sub-style");
error->all(FLERR,"Pair style {} cannot have none as a sub-style", mystyle);
styles[nstyles] = force->new_pair(arg[iarg],1,dummy);
keywords[nstyles] = force->store_style(arg[iarg],0);
@ -345,6 +346,9 @@ void PairHybrid::settings(int narg, char **arg)
nstyles++;
}
if (utils::strmatch(mystyle,"^hybrid/molecular") && (nstyles != 2))
error->all(FLERR, "Pair style {} must have exactly two sub-styles", mystyle);
delete[] cutmax_style;
cutmax_style = new double[nstyles];
memset(cutmax_style, 0, nstyles*sizeof(double));
@ -394,8 +398,7 @@ void PairHybrid::flags()
for (m = 0; m < nstyles; m++) {
if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
if (styles[m]) comm_reverse_off = MAX(comm_reverse_off,
styles[m]->comm_reverse_off);
if (styles[m]) comm_reverse_off = MAX(comm_reverse_off,styles[m]->comm_reverse_off);
}
// single_enable = 1 if all sub-styles are set

1
src/pair_hybrid.h
View File

@ -14,6 +14,7 @@
#ifdef PAIR_CLASS
// clang-format off
PairStyle(hybrid,PairHybrid);
PairStyle(hybrid/omp,PairHybrid);
// clang-format on
#else

162
src/pair_hybrid_molecular.cpp Normal file
View File

@ -0,0 +1,162 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "pair_hybrid_molecular.h"
#include "atom.h"
#include "error.h"
#include "neigh_request.h"
#include "neighbor.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairHybridMolecular::PairHybridMolecular(LAMMPS *lmp) : PairHybridOverlay(lmp) {}
/* ----------------------------------------------------------------------
modify neighbor list requests
------------------------------------------------------------------------- */
void PairHybridMolecular::init_style()
{
if (!atom->molecule_flag)
error->all(FLERR, "Pair style hybrid/molecular requires atom attribute molecule");
if (manybody_flag)
error->all(FLERR, "Pair style hybrid/molecular is not compatible with manybody potentials");
PairHybridOverlay::init_style();
// modify neighbor list requests
bool first = true;
for (auto &request : neighbor->get_pair_requests()) {
if (first) {
request->set_molskip(NeighRequest::INTRA);
first = false;
} else {
request->set_molskip(NeighRequest::INTER);
}
}
born_matrix_enable = 0;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairHybridMolecular::init_one(int i, int j)
{
// if I,J is not set explicitly:
// perform mixing only if I,I sub-style = J,J sub-style
// plus I,I and J,J need the same number of substyles
if (setflag[i][j] == 0) {
if (nmap[i][i] != nmap[j][j])
error->one(FLERR,"All pair coeffs are not set");
int num = 0;
for (int k = 0; k < nmap[i][i]; ++k) {
for (int l = 0; l < nmap[j][j]; ++l) {
if (map[i][i][k] == map[j][j][l]) {
map[i][j][num] = map[i][i][k];
++num;
nmap[i][j] = num;
}
}
}
if (nmap[i][i] != nmap[i][j])
error->one(FLERR,"All pair coeffs are not set");
}
nmap[j][i] = nmap[i][j];
// call init/mixing for all sub-styles of I,J
// set cutsq in sub-style just as Pair::init() does via call to init_one()
// set cutghost for I,J and J,I just as sub-style does
// sum tail corrections for I,J
// return max cutoff of all sub-styles assigned to I,J
// if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned
double cutmax = 0.0;
cutghost[i][j] = cutghost[j][i] = 0.0;
if (tail_flag) etail_ij = ptail_ij = 0.0;
for (int k = 0; k < nmap[i][j]; k++) {
map[j][i][k] = map[i][j][k];
double cut = styles[map[i][j][k]]->init_one(i,j);
if (styles[map[i][j][k]]->did_mix) did_mix = true;
styles[map[i][j][k]]->cutsq[i][j] = styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
if (styles[map[i][j][k]]->ghostneigh)
cutghost[i][j] = cutghost[j][i] = MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
if (tail_flag) {
etail_ij += styles[map[i][j][k]]->etail_ij;
ptail_ij += styles[map[i][j][k]]->ptail_ij;
}
cutmax = MAX(cutmax,cut);
int istyle;
for (istyle = 0; istyle < nstyles; istyle++)
if (styles[istyle] == styles[map[i][j][k]]) break;
if (styles[istyle]->trim_flag) {
if (cut > cutmax_style[istyle]) {
cutmax_style[istyle] = cut;
for (auto &request : neighbor->get_pair_requests()) {
if (styles[istyle] == request->get_requestor()) {
request->set_cutoff(cutmax_style[istyle]);
break;
}
}
}
}
}
return cutmax;
}
/* ----------------------------------------------------------------------
call sub-style to compute single interaction
error if sub-style does not support single() call
since overlay could have multiple sub-styles, sum results explicitly
------------------------------------------------------------------------- */
double PairHybridMolecular::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj, double &fforce)
{
if (nmap[itype][jtype] == 0)
error->one(FLERR,"Invoked pair single() on sub-style none");
double fone;
fforce = 0.0;
double esum = 0.0;
if (nmap[itype][jtype] == 2) {
int m = 0;
if (atom->molecule[i] != atom->molecule[j]) m = 1;
const int mystyle = map[itype][jtype][m];
if (rsq < styles[mystyle]->cutsq[itype][jtype]) {
if (styles[mystyle]->single_enable == 0)
error->one(FLERR,"Pair hybrid/molecular sub-style {} does not support single() call",
keywords[mystyle]);
if ((special_lj[mystyle] != nullptr) || (special_coul[mystyle] != nullptr))
error->one(FLERR,"Pair hybrid/molecular single() calls do not support per sub-style "
"special bond values");
esum += styles[mystyle]->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone);
fforce += fone;
}
}
if (single_extra) copy_svector(itype,jtype);
return esum;
}

41
src/pair_hybrid_molecular.h Normal file
View File

@ -0,0 +1,41 @@
/* -*- 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 PAIR_CLASS
// clang-format off
PairStyle(hybrid/molecular,PairHybridMolecular);
PairStyle(hybrid/molecular/omp,PairHybridMolecular);
// clang-format on
#else
#ifndef LMP_PAIR_HYBRID_MOLECULAR_H
#define LMP_PAIR_HYBRID_MOLECULAR_H
#include "pair_hybrid_overlay.h"
namespace LAMMPS_NS {
class PairHybridMolecular : public PairHybridOverlay {
public:
PairHybridMolecular(class LAMMPS *);
void init_style() override;
double init_one(int, int) override;
double single(int, int, int, int, double, double, double, double &) override;
};
} // namespace LAMMPS_NS
#endif
#endif

1
src/pair_hybrid_overlay.h
View File

@ -14,6 +14,7 @@
#ifdef PAIR_CLASS
// clang-format off
PairStyle(hybrid/overlay,PairHybridOverlay);
PairStyle(hybrid/overlay/omp,PairHybridOverlay);
// clang-format on
#else

1
src/pair_hybrid_scaled.h
View File

@ -14,6 +14,7 @@
#ifdef PAIR_CLASS
// clang-format off
PairStyle(hybrid/scaled,PairHybridScaled);
PairStyle(hybrid/scaled/omp,PairHybridScaled);
// clang-format on
#else

100
unittest/force-styles/tests/mol-pair-hybrid_molecular.yaml Normal file
View File

@ -0,0 +1,100 @@
---
lammps_version: 17 Apr 2024
tags:
date_generated: Sun Jun 9 11:41:13 2024
epsilon: 1e-13
skip_tests:
prerequisites: ! |
atom full
pair lj/cut
pair hybrid/molecular
pre_commands: ! ""
post_commands: ! |
pair_modify mix arithmetic
pair_modify shift yes
input_file: in.fourmol
pair_style: hybrid/molecular lj/cut 8.0 lj/cut 8.0
pair_coeff: ! |
1 1 lj/cut 1 0.02 2.5
2 2 lj/cut 1 0.005 1.0
2 4 lj/cut 1 0.005 0.5
3 3 lj/cut 1 0.02 3.2
4 4 lj/cut 1 0.015 3.1
5 5 lj/cut 1 0.015 3.1
1 1 lj/cut 2 0.02 2.5
2 2 lj/cut 2 0.005 1.0
2 4 lj/cut 2 0.005 0.5
3 3 lj/cut 2 0.02 3.2
4 4 lj/cut 2 0.015 3.1
5 5 lj/cut 2 0.015 3.1
extract: ! ""
natoms: 29
init_vdwl: 749.2470096189497
init_coul: 0
init_stress: ! |2-
2.1793857186503237e+03 2.1988957679770601e+03 4.6653994738862348e+03 -7.5956544622684351e+02 2.4751393539193327e+01 6.6652061873806679e+02
init_forces: ! |2
1 -2.3333390274530565e+01 2.6994567613591136e+02 3.3272827850621587e+02
2 1.5828554630423909e+02 1.3025008843536872e+02 -1.8629682358915150e+02
3 -1.3528903744071786e+02 -3.8704313350789641e+02 -1.4568978426110147e+02
4 -7.8711096705734178e+00 2.1350518625352004e+00 -5.5954532185292400e+00
5 -2.5176757267276137e+00 -4.0521510680612876e+00 1.2152704057983799e+01
6 -8.3190665562047559e+02 9.6394165349388811e+02 1.1509101492424440e+03
7 5.8203416066164465e+01 -3.3609013622052356e+02 -1.7179626006587687e+03
8 1.4451392646293456e+02 -1.0927476052490437e+02 3.9990594285329479e+02
9 7.9156945283109025e+01 8.5273009784086469e+01 3.5032175698457496e+02
10 5.3118875219106917e+02 -6.1040990846582008e+02 -1.8355872692632028e+02
11 -2.3530157265571860e+00 -5.9077640075588906e+00 -9.6590723956614433e+00
12 1.7527155197359406e+01 1.0633119514682473e+01 -7.9254397903886158e+00
13 8.0986409580712841e+00 -3.2098088269317300e+00 -1.4896399871387667e-01
14 -3.3852721291218524e+00 6.8636181224987947e-01 -8.7507190862837803e+00
15 -2.0454999188607306e-01 8.4846165523012136e+00 3.0131615419840623e+00
16 4.6326331471561195e+02 -3.3087730492363477e+02 -1.1893030175606582e+03
17 -4.5334322060634048e+02 3.1554297967975305e+02 1.2058423415744451e+03
18 -1.8862629870158503e-02 -3.3402022492930034e-02 3.1000492146377390e-02
19 3.1843079948447594e-04 -2.3918628211596124e-04 1.7427252652160224e-03
20 -9.9760831169755002e-04 -1.0209184785886856e-03 3.6910973051849135e-04
21 -7.1566158640374354e+01 -8.1615716383825756e+01 2.2589571940670788e+02
22 -1.0808840769631149e+02 -2.6193799449067580e+01 -1.6957912849816358e+02
23 1.7964463850759611e+02 1.0782102722442450e+02 -5.6305812731665995e+01
24 3.6591423637378952e+01 -2.1181597497621908e+02 1.1218307103182993e+02
25 -1.4851496072162067e+02 2.3907129270267117e+01 -1.2485640694398953e+02
26 1.1191134671510581e+02 1.8789783424990625e+02 1.2650143102803206e+01
27 5.1810412832328005e+01 -2.2705468907750404e+02 9.0849153441059272e+01
28 -1.8041315533250560e+02 7.7534079082878250e+01 -1.2206962452216493e+02
29 1.2861063251415737e+02 1.4952718246094852e+02 3.1216040111076957e+01
run_vdwl: 719.4532389988315
run_coul: 0
run_stress: ! |2-
2.1330157554553725e+03 2.1547730555430494e+03 4.3976512412988704e+03 -7.3873325485023713e+02 4.1743707190785464e+01 6.2788040986774649e+02
run_forces: ! |2
1 -2.0299419744961813e+01 2.6686193379336868e+02 3.2358785871037435e+02
2 1.5298617928501707e+02 1.2596516341411086e+02 -1.7961292655320207e+02
3 -1.3353630670276331e+02 -3.7923748676909099e+02 -1.4291839777232488e+02
4 -7.8374717836014449e+00 2.1276610789788282e+00 -5.5845014473593917e+00
5 -2.5014258629959465e+00 -4.0250131424457543e+00 1.2103512372172734e+01
6 -8.0681466162480240e+02 9.2165651041424792e+02 1.0270802401119468e+03
7 5.5780302775854636e+01 -3.1117544157318952e+02 -1.5746997989226002e+03
8 1.3452983973683914e+02 -1.0064660034658647e+02 3.8851792520911863e+02
9 7.6746213900459239e+01 8.2501469902247337e+01 3.3944351209160595e+02
10 5.2128033526109800e+02 -5.9920098832868109e+02 -1.8126029871233905e+02
11 -2.3573118088794365e+00 -5.8616944553482799e+00 -9.6049808813641686e+00
12 1.7503975897697526e+01 1.0626930302269722e+01 -8.0603160114673926e+00
13 8.0530313324242382e+00 -3.1756495175042607e+00 -1.4618315691984204e-01
14 -3.3416065166863160e+00 6.6492606318663194e-01 -8.6345131440736740e+00
15 -2.2253843262483208e-01 8.5025661635305205e+00 3.0369735873547175e+00
16 4.3476329769010198e+02 -3.1171099668258080e+02 -1.1135222104230593e+03
17 -4.2469864617016134e+02 2.9615424659116553e+02 1.1302578406458213e+03
18 -1.8849988250623853e-02 -3.3371648038832503e-02 3.0986306282264790e-02
19 3.0940278115793517e-04 -2.4634536779368854e-04 1.7433360016754921e-03
20 -9.8648131231171901e-04 -1.0112587092668940e-03 3.6932949186791977e-04
21 -7.0490777148272102e+01 -7.9749189729874402e+01 2.2171013458550721e+02
22 -1.0638722739944252e+02 -2.5949513934649758e+01 -1.6645597092015180e+02
23 1.7686805727889882e+02 1.0571023691370021e+02 -5.5243362166860535e+01
24 3.8206035227327128e+01 -2.1022829679057401e+02 1.1260716393332925e+02
25 -1.4918888258035886e+02 2.3762162241718102e+01 -1.2549193847418989e+02
26 1.1097064525776705e+02 1.8645512086371158e+02 1.2861565481437628e+01
27 5.0800867695850577e+01 -2.2296598219372004e+02 8.8607407764830413e+01
28 -1.7694198509380672e+02 7.6029979926844575e+01 -1.1950523558040683e+02
29 1.2614900659680345e+02 1.4694257504728049e+02 3.0893400701043564e+01
...

119
unittest/force-styles/tests/mol-pair-hybrid_molecular_lj_morse.yaml Normal file
View File

@ -0,0 +1,119 @@
---
lammps_version: 17 Apr 2024
tags:
date_generated: Sun Jun 9 11:41:13 2024
epsilon: 1e-13
skip_tests:
prerequisites: ! |
atom full
pair lj/cut
pair morse
pair hybrid/molecular
pre_commands: ! ""
post_commands: ! |
pair_modify mix arithmetic
pair_modify shift yes
input_file: in.fourmol
pair_style: hybrid/molecular lj/cut 8.0 morse 8.0
pair_coeff: ! |
1 1 lj/cut 0.02 2.5
1 2 lj/cut 0.01 1.75
1 3 lj/cut 0.02 2.85
1 4 lj/cut 0.01732050807568877293 2.8
1 5 lj/cut 0.01732050807568877293 2.8
2 2 lj/cut 0.005 1.0
2 3 lj/cut 0.01 2.1
2 4 lj/cut 0.005 0.5
2 5 lj/cut 0.00866025403784438646 2.05
3 3 lj/cut 0.02 3.2
3 4 lj/cut 0.01732050807568877293 3.15
3 5 lj/cut 0.01732050807568877293 3.15
4 4 lj/cut 0.015 3.1
4 5 lj/cut 0.015 3.1
5 5 lj/cut 0.015 3.1
1 1 morse 0.0202798941614106 2.78203488021395 2.725417159299
1 2 morse 0.0101167811264648 3.9793050302425 1.90749569018897
1 3 morse 0.0202934330695928 2.43948720203264 3.10711749999622
1 4 morse 0.0175731334238374 2.48316585521317 3.05258880102438
1 5 morse 0.0175731334238374 2.48316585521317 3.05258880102438
2 2 morse 0.00503064360487288 6.98433077606902 1.08960295117864
2 3 morse 0.0101296013842819 3.31380153807866 2.28919067558352
2 4 morse 0.00497405122588691 14.0508902925745 0.544416409093563
2 5 morse 0.00877114211614446 3.39491256196178 2.23466262511073
3 3 morse 0.0203039874239943 2.17204344301477 3.48881895084762
3 4 morse 0.0175825321440736 2.20660439192238 3.43428999287994
3 5 morse 0.0175825321440736 2.20660439192238 3.43428999287994
4 4 morse 0.0152259201379927 2.24227873774009 3.37976131582396
4 5 morse 0.0152259201379927 2.24227873774009 3.37976131582396
5 5 morse 0.0152259201379927 2.24227873774009 3.37976131582396
extract: ! ""
natoms: 29
init_vdwl: 642.2857035487534
init_coul: 0
init_stress: ! |2-
1.5446101412530770e+03 1.7762864877230827e+03 4.3406629723991382e+03 -2.4168966825090698e+02 -4.2399281241000449e+02 1.0313398732648857e+03
init_forces: ! |2
1 -2.3337658999820331e+01 2.6994849929388920e+02 3.3272731927204762e+02
2 1.5828475525620013e+02 1.3025198646645657e+02 -1.8629727921784635e+02
3 -1.3530883593401191e+02 -3.8702990176780906e+02 -1.4568955729804821e+02
4 -7.8720974048138279e+00 2.1357285334450031e+00 -5.5956697614574145e+00
5 -2.5179872246085657e+00 -4.0521698308193113e+00 1.2152426584580066e+01
6 -1.0627539753780246e+02 3.7096636535889553e+02 1.6577829288882351e+03
7 6.2545311732649182e+01 -3.3857732069504033e+02 -1.7085638850072914e+03
8 -5.8544805218664203e+02 4.8619364551174579e+02 -1.1637016110669298e+02
9 7.9157108381891348e+01 8.5271268694585373e+01 3.5031950365280102e+02
10 5.3119792393751504e+02 -6.1042010285634149e+02 -1.8355328074325874e+02
11 -2.3525975498215241e+00 -5.9087312467830966e+00 -9.6592450379633981e+00
12 1.7514418718210869e+01 1.0633490551535797e+01 -7.9269392455530774e+00
13 8.0989976759274995e+00 -3.2092537736156967e+00 -1.4857368826935238e-01
14 -3.3831155658252769e+00 6.8583230320825528e-01 -8.7521044176224088e+00
15 -2.0078689909299524e-01 8.4842702953333049e+00 3.0127255727602438e+00
16 4.6326713730462569e+02 -3.3088113342819690e+02 -1.1892994362776039e+03
17 -4.5333858978430715e+02 3.1553830929611308e+02 1.2058459850635415e+03
18 -2.0452760511035117e-02 -3.1643654311568722e-02 2.7038063904994043e-02
19 6.7683676563436612e-04 5.5257878746558435e-04 3.4368881128008698e-04
20 -7.2687428125721442e-04 -6.4369255876917293e-04 -7.3055429537143658e-05
21 -7.1559739956792001e+01 -8.1623086552077027e+01 2.2588907386171752e+02
22 -1.0808827580116072e+02 -2.6193999003846525e+01 -1.6957935616989226e+02
23 1.7964485683675787e+02 1.0782085824265097e+02 -5.6306005220676077e+01
24 3.6592705490564995e+01 -2.1181065162415683e+02 1.1218892889291895e+02
25 -1.4851386780240065e+02 2.3908077402265249e+01 -1.2485428291933182e+02
26 1.1191162087301427e+02 1.8789805761651786e+02 1.2650515090480996e+01
27 5.1805674347396028e+01 -2.2705907095144718e+02 9.0852216146745278e+01
28 -1.8041330407833553e+02 7.7533936594765862e+01 -1.2206953762508090e+02
29 1.2861029896870866e+02 1.4952683033680830e+02 3.1216382013473979e+01
run_vdwl: 612.532027437114
run_coul: 0
run_stress: ! |2-
1.5079612443981077e+03 1.7300922919292559e+03 4.0645404509935297e+03 -2.2501656863720930e+02 -3.8440649220519191e+02 9.6872117407090468e+02
run_forces: ! |2
1 -2.0303664163328520e+01 2.6686344032747223e+02 3.2358595259116754e+02
2 1.5298537986232245e+02 1.2596704781215615e+02 -1.7961336269199933e+02
3 -1.3340899705451542e+02 -3.7937983996408559e+02 -1.4304731523709151e+02
4 -7.8384621494638589e+00 2.1283385602122546e+00 -5.5847156581558837e+00
5 -2.5017213737297879e+00 -4.0250398342732856e+00 1.2103221315331954e+01
6 -1.0069389067343234e+02 3.4070704568948952e+02 1.4995989516904360e+03
7 5.6637967265264535e+01 -3.0800622792471370e+02 -1.5506690902144287e+03
8 -5.6211472712331181e+02 4.6807109251305587e+02 -1.0712697415342610e+02
9 7.6069407842570087e+01 8.0619386602314336e+01 3.3418245682682038e+02
10 5.1153094792744940e+02 -5.8811350830044319e+02 -1.7678443640175220e+02
11 -2.3568347685072353e+00 -5.8621926046697821e+00 -9.6044928327397194e+00
12 1.7464493846680174e+01 1.0623724119993449e+01 -8.0428896399820058e+00
13 8.0533880511151636e+00 -3.1750913604087598e+00 -1.4579751155151741e-01
14 -3.3394709175794697e+00 6.6440083564233665e-01 -8.6358918170724337e+00
15 -2.1889102849274136e-01 8.5021207266300767e+00 3.0365920030191385e+00
16 4.3476429812120438e+02 -3.1171100998816553e+02 -1.1135203700101538e+03
17 -4.2469867717957146e+02 2.9615706504039451e+02 1.1302529422580635e+03
18 -2.0438484585882859e-02 -3.1629705308846789e-02 2.7018991453856001e-02
19 6.6852492737214141e-04 5.4625495751892833e-04 3.4152455215868423e-04
20 -7.1620896561888831e-04 -6.3500663583982407e-04 -7.1700429767679216e-05
21 -7.0484355349417541e+01 -7.9756579011007801e+01 2.2170349595393003e+02
22 -1.0638710558378827e+02 -2.5949713015442118e+01 -1.6645620880910886e+02
23 1.7686828769914484e+02 1.0571007956441812e+02 -5.5243556600722584e+01
24 3.8207315780355771e+01 -2.1022299522635379e+02 1.1261302962088173e+02
25 -1.4918780612527124e+02 2.3763111781582175e+01 -1.2548982817638399e+02
26 1.1097093258138914e+02 1.8645536456663984e+02 1.2861936841695268e+01
27 5.0796145560529055e+01 -2.2297034959629724e+02 8.8610475497234873e+01
28 -1.7694213444502824e+02 7.6029840540845129e+01 -1.1950515191029582e+02
29 1.2614865956603674e+02 1.4694220660200173e+02 3.0893738250707873e+01
...